29
$\begingroup$

Given an output from optim with a hessian matrix, how to calculate parameter confidence intervals using the hessian matrix?

fit<-optim(..., hessian=T)

hessian<-fit$hessian

I am mostly interested in the context of maximum likelihood analysis, but am curious to know if the method can be expanded beyond.

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ This question is too vague. What kind of confidence intervals? What sort of models are you interested on? Consider taking a look at the literature before posting 3 questions in a row. $\endgroup$
    – user10525
    Apr 24, 2012 at 14:35
  • $\begingroup$ The method should be independent of the interval and model. $\endgroup$
    – Etienne Low-Décarie
    Apr 24, 2012 at 15:25
  • $\begingroup$ What function do you optimize ? $\endgroup$
    – Stéphane Laurent
    Apr 24, 2012 at 17:36
  • $\begingroup$ I was told I should be able to do this independent of the model I use. $\endgroup$
    – Etienne Low-Décarie
    Apr 25, 2012 at 20:38
  • $\begingroup$ Corey Chivers provided the answer. $\endgroup$
    – Etienne Low-Décarie
    Apr 25, 2012 at 21:12

1 Answer 1

Reset to default
41
$\begingroup$

If you are maximising a likelihood then the covariance matrix of the estimates is (asymptotically) the inverse of the negative of the Hessian. The standard errors are the square roots of the diagonal elements of the covariance (from elsewhere on the web!, from Prof. Thomas Lumley and Spencer Graves, Eng.).

For a 95% confidence interval

fit<-optim(pars,li_func,control=list("fnscale"=-1),hessian=TRUE,...)
fisher_info<-solve(-fit$hessian)
prop_sigma<-sqrt(diag(fisher_info))
prop_sigma<-diag(prop_sigma)
upper<-fit$par+1.96*prop_sigma
lower<-fit$par-1.96*prop_sigma
interval<-data.frame(value=fit$par, upper=upper, lower=lower)

Note that:

  • If you are maximizing the log(likelihood), then the NEGATIVE of the hessian is the "observed information" (such as here).
  • If you MINIMIZE a "deviance" = (-2)*log(likelihood), then the HALF of the hessian is the observed information.
  • In the unlikely event that you are maximizing the likelihood itself, you need to divide the negative of the hessian by the likelihood to get the observed information.

See this for further limitations due to optimization routine used.

Improve this answer
$\endgroup$
12
  • 3
    $\begingroup$ (+1) The inverse of the negative hessian is an estimator of the asymptotic covariance matrix - I know this appears in your code but I think it's important to point out. $\endgroup$
    – Macro
    Apr 26, 2012 at 12:20
  • 1
    $\begingroup$ Excellent answer, should the upper and lower bounds read upper<-fit$par+1.96*(prop_sigma/sqrt(n)) lower<-fit$par-1.96*(prop_sigma/sqrt(n)) ? Thanks $\endgroup$
    – forecaster
    Jun 16, 2015 at 0:03
  • 5
    $\begingroup$ Why not delete line 4? $\endgroup$
    – Jason
    Sep 12, 2016 at 19:34
  • 3
    $\begingroup$ Is including the line prop_sigma<-diag(prop_sigma) an error? $\endgroup$
    – Mark Miller
    Jan 21, 2017 at 2:58
  • 1
    $\begingroup$ FWIW numDeriv::hessian() might be worth a look. stat.ethz.ch/pipermail/r-help/2009-August/401566.html $\endgroup$
    – swihart
    Apr 1, 2020 at 14:50

Your Answer

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

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