# How to get asymptotic covariance matrix when observed information matrix is singular

I'm fitting different models by Maximum Likelihood. To do this I'm using a stochastic version of Newton-Raphson algorithm, where both the gradient and the Hessian of the likelihood are estimated at each step.

In most cases I reach convergence, but then I often encounter the following problem: the estimate Hessian $\hat H$ at convergence is negative definite. This is a problem because I can't invert it to get an estimate of the observed Fisher information.

This happens because one or more parameters are weakly identified. On the other hand other parameters seem to be well identified and their corresponding entries in the estimated Hessian seem to make sense. Identifiability is difficult to assess beforehand for the dynamical models I'm trying to fit.

What I would like is an approach that (starting from the $\hat H$) points out what set of parameters are not identifiable and that gives variances for the remaining parameters.

I tried to do this by tilting the smallest eigenvalue of $\hat H$ in order to get a better conditioning number. The results seem arbitrary: depending on the conditioning number I want to achieve I get different variances for the parameters.

EDIT: This is a typical example of an Hessian I can end up with:

H <- matrix(c( 67.23586, 10.477815, 138.696877,
10.47782, -3.238982,   2.592774,
138.69688,  2.592774, 473.161347 ), 3, 3, byrow = TRUE)


You see that I have a negative entry in the main diagonal: the second parameter is weakly identifiable, while the other well identified. The other situation I often encounter is that of sub-group of parameters that are redundant - highly correlated.

• Can you provide a reproducible example? I think you could probably use singular value decomposition to decompose the Hessian into singular and non-singular terms ... (or equivalently [???], compute a generalized inverse, as in the ginv function of the MASS package in R ...) – Ben Bolker Nov 6 '13 at 23:14
• That sound like something I might be very interested in! I added a simple example to the question. – Matteo Fasiolo Nov 6 '13 at 23:50
• You haven't said how you're estimating the Hessian. If you're using a quasi-Newton method (such as BFGS) you should be aware that the matrices produced by the quasi-Newton methods are often poor estimates of the actual Hessian and they shouldn't be relied upon for doing things like computing the Fisher information. – Brian Borchers Nov 7 '13 at 3:49
• I'm basically approximating the Hessian and Gradient by using local quadratic regression, where the regression coefficients are estimates of the first and second derivatives of the likelihood. I can always increase the sample size to get better estimates, but if a parameter is not identifiable the likelihood will still be flat in that direction. – Matteo Fasiolo Nov 7 '13 at 10:04
• This is a suggestion only. If the parameter is weakly identified, i.e. log-likelihood does not change when the parameter is changed, then the partial derivative of log-likelihood of said parameter should be close to zero for any values of that parameter. This way you could identify bad parameters before the optimisation. – mpiktas Nov 11 '13 at 13:19

## 2 Answers

It is less general than your problem, but I think this article may help: Ridge Estimators in Logistic Regression. Basically, instead of maximizing the (log-)likelihood $L$, you maximize $$L_\lambda(\beta) = L(\beta) - \lambda \|\beta\|_2^2 .$$ This regularization increases all the eigenvalues of the Hessian matrix, so it is not so far away from what you are doing right now.

I post here my solution to this problem even though it more of a quick fix than real solution. It seems to work quite well when there are one or more uncorrelated parameters that are not identifiable, not sure how well it does when there is a set of highly correlated parameters.

Start with an initial Hessian $H_0$ and index $k = 0$ and than follow this procedure:

1. At iteration $k$ verify whether $H_k$ is positive definite "enough" by checking whether the ratio between smallest and biggest eigenvalues is above a certain tolerance. If it is then stop otherwise go to step to 2.

2. Increase the small eigenvalues so that the matrix is PD according to the above criterion.

3. Invert the modified Hessian in order to get the corresponding covariance matrix $\Sigma_k$.

4. Find the parameter with the highest variance (the largest diagonal entry of $\Sigma_k$) and eliminate the corresponding row and column in $H_k$ in order to obtain $H_{k+1}$. Go back to step 1.

In step 4 the idea is that the parameter with the highest variance is the least identifiable. Given that the eigenvectors corresponding to smallest eigenvalues lay mostly in its direction, if we get rid of it the conditioning number of the Hessian should improve.

In the simple case of the above $3 \times 3$ matrix we get:

CH <- cleanHessian(H)

CH
$hessian [,1] [,2] [1,] 67.23586 138.6969 [2,] 138.69688 473.1613$badParam
 2


The function output are a Hessian where the entries corresponding to the unidentifiable parameters have been removed, and a vector containing the indexes of these parameters. In this case the second parameter was not identified. I can now invert the reduced Hessian to get the asymptotic covariance matrix:

solve(CH$hessian) [,1] [,2] [1,] 0.03762240 -0.011028182 [2,] -0.01102818 0.005346114  According to my experience with the model I'm fitting here these covariances are credible. Here is the R code for the procedure:  cleanHessian <- function(hessian, tol = 1e-8) { badParam <- numeric() repeat{ eDec <- eigen(hessian) # Find small eigenvalues lowEigen <- which(eDec$values < eDec$values * tol) # If there are none, I consider the hessian to be PD and I exit if( length(lowEigen) == 0 ) break # Increase the small eigenvalues to a threshold eDec$values[lowEigen] <- eDec$values * tol # Invert the modified hessian to get the covariance COV <- eDec$vectors%*%(t(eDec$vectors) / eDec$value)

# I identify the parameter with the highest variance, I remove the corresponding
# elements from the Hessian and I store its index
bad <- which.max( diag(COV) )
hessian <- hessian[-bad, -bad]
offset <- sum( badParam <= bad )

badParam <- c(badParam, bad + offset)

}

return( list("hessian" = hessian, "badParam" = badParam) )
}