Say you have $y = f(\mathbf{X, \theta})$ and you want $\partial y_i/\partial \mathbf{\theta}$ at the optimized value of $\theta$. Further assume that it is feasible -- but a hassle -- to calculate these derivatives analytically, and simple to compute a numerical Jacobian, using for example jacobian in the numDeriv package in R.
Under what conditions is the analytical approach preferable? Given that these are just first derivatives -- slope coefs -- it seems like the numerical approach should always work just fine, and avoid any risk of pencil-and-paper mistakes. Is this correct?