2
$\begingroup$

I have run an MLE estimation using R to estimate parameters $\theta_1, \theta_2, \theta_3$, which are:

 0.0002022146  0.0222026625 -6.1910421067

the variance matrix is:

0.00000004574446 0.000001708064 0.0000008477115
0.00000170806367 0.000115362174 0.0000572542168
0.00000084771147 0.000057254217 0.0024674380942

Now, I need to simulate standard deviations for parameters $\kappa=\frac{\theta_1}{\theta_2}$. I have tried:

library(MASS)
df <- mvrnorm(n = 1000000, mu = fit@coef, Sigma = fit@vcov)
kappa<-df[,1]/df[,2]
sdkappa <- sd(kappa)

Is this the correct way? I doubt it, since the mean(kappa) is quite different from $\frac{\hat{\theta_1}}{\hat{\theta_2}}$. What is the proper way of doing this?

$\endgroup$

1 Answer 1

2
$\begingroup$

Briefly: your approach is reasonable. It assumes that the sampling distribution of the MLE is approximately multivariate Normal, and it only makes sense (more generally the ratio of Normals only makes sense) if the mean of the denominator is far from zero ($\textrm{CV}_{\textrm{denom}}$ is large): see this blog post by Andrew Gelman. In your case $\textrm{CV}_2$ is $0.022/\sqrt{0.00011} \approx 2.09$, which is probably OK (I would really worry if CV were $<1$).

  • I don't know the formal name for this approach: in Ecological Models and Data (2008) ($\S$7.5.3) I call the confidence intervals derived in this way "population prediction intervals".
  • More approximately you can use the delta method (search here and elsewhere).
  • In general $f(\bar x, \bar y)$ is not the same as $\overline{f(x,y)}$: this is called Jensen's Inequality. (There is a delta-method-like approximation for estimating the magnitude of JI in particular cases.)
  • if the distribution of $\theta_1/\theta_2$ is markedly non-Normal, the standard deviation may not be as meaningful as you think (try computing the (0.025,0.975) quantiles of your sampled value and comparing it to the confidence intervals based on $\pm 1.96 \sigma$). (It would be useful to know why you want to calculate standard deviations of $\theta_1/\theta_2$, i.e. how do you plan to use this quantity going forward?)

Other ways of doing this that are less approximate include:

  • re-parameterizing the model to use $\theta_r= \theta_1/\theta_2$ as a primary parameter (then finding profile confidence intervals for $\theta_r$)
  • bootstrapping confidence intervals for $\theta_r$
  • Bayesian approaches (e.g. MCMC) to find the posterior marginal credible intervals of $\theta_r$ based on samples of $\{\theta_1, \theta_2\}$.
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.