Matrix repair

When fitting a statistical model we want two things as a minimum:

1. The parameter estimates, e.g. the maximum-likelihood estimates (MLEs), and
2. The estimated variance-covariance matrix, \(\hat V\), for those estimates.

We can get both from the log-likelihood: the MLEs maximise the value of the log-likelihood function, and an approximation for the covariance matrix comes from inverting the negative information matrix, \(\mathcal J\), i.e. the matrix of second partial derivatives evaluated at the MLEs. However, the limitations of computer arithmetic can sometimes get in the way, as shown in the information matrix in Figure 1:

Figure 1. Information matrix, \(\mathcal J\), for a five-parameter model. Only the leading diagonal is shown for clarity.  Source: own calculations.

\[\mathcal J=\left(\begin{matrix}-11584.8 & \circ & \circ & \circ & \circ \\
\circ & -10062.7 & \circ & \circ & \circ \\
\circ & \circ & -167.653 & \circ & \circ \\
\circ & \circ & \circ & -7275.98 & \circ \\
\circ & \circ & \circ & \circ & -0.00291038\end{matrix}\right)\]

The issue with \(\mathcal J\) in Figure 1 lies in the bottom right: the pure second partial derivative \(\mathcal J_{5,5}\) is tiny relative to the pure second partial derivatives for the other parameters. This causes a problem when inverting to estimate the covariance matrix, as shown in Figure 2:

Figure 2. \(\hat V=-\mathcal J^{-1}\) for the five-parameter model. Only the leading diagonal is shown for clarity.  Source: own calculations.

\[-\mathcal J^{-1}=\left(\begin{matrix}-0.000201 & \circ & \circ & \circ & \circ \\
\circ & -0.000181 & \circ & \circ & \circ \\
\circ & \circ & -00642 & \circ & \circ \\
\circ & \circ & \circ & -0.000136 & \circ \\
\circ & \circ & \circ & \circ & -3980.715\end{matrix}\right)\]

The leading diagonal of \(\hat V\) is supposed to be the Cramer-Rao lower bound for the variances of the parameter estimates, i.e. they should all be positive values.  However, the unreliable value for \(\mathcal J_{5,5}\) has rippled through the entire matrix inversion and ruined the estimates for everything else.  Can anything be done to recover the situation?

One avenue of exploration is to try alternative values for \(\mathcal J_{5,5}\) and see what the impact is on the leading diagonal of \(\hat V\).  Table 1 shows the estimated variances for a series of alternative values of \(\mathcal J_{5,5}\):

Table 1.  Sensitivity of variance estimates to alternative values for \(\mathcal J_{5,5}\).  Source: own calculations using the R script on the right.

Table 1 shows that the variance estimates for the first four parameters are robust to a wide range of alternative values of \(\mathcal J_{5,5}\). Inspection of the resulting coefficients of variation (not shown) suggests that these are indeed sensible variance estimates. Thus, it is possible to recover something from an information matrix even where one of the values is unreliable.  The R script on the right gives the full details of this example; feel free to adapt it when trying to repair your own matrices!