Standard Errors From the Observed Information in Joint Models

Why Standard Errors Need the Whole Joint Model

A joint model estimates longitudinal and event-time parameters together, often while integrating over latent subject-specific random effects. An estimate can be numerically available without its uncertainty being clear: standard errors depend not only on the curvature of the likelihood in one parameter, but also on how that parameter trades off against every other estimated parameter. The observed information matrix captures this local curvature and dependence.

Let $D$ denote the observed longitudinal and time-to-event data, and let $\theta$ be the vector of all model parameters, in a stated parameterization. The integrated likelihood is obtained by integrating the subject-specific random effects out of the joint likelihood. For independent subjects, it can be written as a product of subject-level integrated contributions. Write $\ell(\theta)=\log L(\theta;D)$ for its log-likelihood and $\hat{\theta}$ for an interior maximum-likelihood estimate.

Definition: The observed information matrix at $\hat{\theta}$ is the negative Hessian of the integrated log-likelihood, $$ J(\hat{\theta})=-\left.\frac{\partial^2\ell(\theta)}{\partial\theta\,\partial\theta^\mathsf{T}}\right|_{\theta=\hat{\theta}}. $$ If $J(\hat{\theta})$ is nonsingular and positive definite, the usual approximate covariance matrix of $\hat{\theta}$ is $J(\hat{\theta})^{-1}$. The approximate standard error of component $\hat{\theta}_j$ is $\sqrt{[J(\hat{\theta})^{-1}]_{jj}}$.

The derivatives are of the likelihood after the random effects have been integrated out, not of a likelihood that treats their unknown realized values as additional fixed parameters. The earlier material on observed information and Louis’ identity gives a route to the observed information in latent-variable settings; direct differentiation of the integrated log-likelihood is another route. Whichever route is used, the matrix must correspond to the likelihood and parameter vector used to obtain the reported estimates.

Why the Inverse Determines the Standard Errors

The approximation has a local likelihood justification. At an interior maximum, the score is zero. For a small displacement $h=\theta-\hat{\theta}$, a second-order Taylor expansion gives

\ell(\hat{\theta}+h) \approx \ell(\hat{\theta})+ \left.\frac{\partial\ell}{\partial\theta}\right|_{\hat{\theta}}^\mathsf{T}h +\frac{1}{2}h^\mathsf{T} \left.\frac{\partial^2\ell}{\partial\theta\,\partial\theta^\mathsf{T}}\right|_{\hat{\theta}}h = \ell(\hat{\theta})-\frac{1}{2}h^\mathsf{T}J(\hat{\theta})h.

When the likelihood is locally well approximated by this quadratic form, its shape around the maximum is like a multivariate normal density with covariance $J(\hat{\theta})^{-1}$. This gives the standard errors from the diagonal entries of the inverse. The off-diagonal entries are approximate covariances between estimates. They matter: correlated parameter estimates can have different marginal variances from what one would infer by considering each diagonal information entry alone.

This is an approximation, not an exact finite-sample distribution. Its usual justification relies on regularity conditions such as a correctly specified, identifiable model; a sufficiently informative sample; differentiability of the integrated log-likelihood; and a true parameter not on a problematic boundary. In joint models, the likelihood and its derivatives must also be evaluated accurately enough that integration or optimization error does not dominate the curvature calculation. If censoring is not explicitly modeled, the likelihood must be appropriate under the assumed conditionally independent censoring mechanism.

Worked Example: Inverting a Two-Parameter Information Matrix

Suppose a hypothetical joint-model fit estimates two dimensionless coefficients: $\beta$, a longitudinal fixed effect for a standardized predictor, and $\alpha$, an association coefficient for the standardized longitudinal marker in the event model. At the maximum, the observed information in the $(\beta,\alpha)$ parameterization is

J(\hat{\theta})= \begin{pmatrix} 100 & 20\\ 20 & 25 \end{pmatrix}, \qquad \hat{\beta}=0.70,\quad \hat{\alpha}=0.35.

The determinant is $100(25)-20(20)=2{,}100$, so the inverse is

J(\hat{\theta})^{-1} = \frac{1}{2{,}100} \begin{pmatrix} 25 & -20\\ -20 & 100 \end{pmatrix} = \begin{pmatrix} 0.01190 & -0.009524\\ -0.009524 & 0.04762 \end{pmatrix}.

Thus $SE(\hat{\beta})=\sqrt{0.01190}=0.1091$, and $SE(\hat{\alpha})=\sqrt{0.04762}=0.2182$. For an approximate 95% Wald interval, use the estimate plus or minus $1.96$ standard errors. The intervals are $0.70\pm1.96(0.1091)=[0.4862,0.9138]$ for $\beta$, and $0.35\pm1.96(0.2182)=[-0.0777,0.7777]$ for $\alpha$, with endpoints rounded. These intervals describe uncertainty in the dimensionless coefficients on their stated scales; the interval for $\alpha$ includes zero.

The covariance is $-0.009524$, and the estimated correlation is $-0.009524/\sqrt{(0.01190)(0.04762)}=-0.400$. This dependence is why the standard errors come from the inverse matrix rather than from reciprocals of the diagonal information entries.

Parameter Scales and Transformations

An information matrix is tied to the coordinates in which it is calculated. A standard error for a log standard deviation is not a standard error for the standard deviation itself. Likewise, a coefficient on a log-hazard scale is not a standard error for a hazard ratio. To report uncertainty on another scale, transform the covariance appropriately.

For a scalar quantity $g(\theta)$, the delta method uses a first-order Taylor expansion around $\hat{\theta}$. If $V=J(\hat{\theta})^{-1}$, then

\widehat{\operatorname{Var}}\{g(\hat{\theta})\} \approx \nabla g(\hat{\theta})^\mathsf{T}V\nabla g(\hat{\theta}).

For a vector-valued transformation, use its Jacobian $G$ to obtain the approximate covariance $GVG^\mathsf{T}$. The transformation must match the reported quantity: exponentiating a coefficient changes both its interpretation and its uncertainty scale.

Worked Example: Standard Error for a Standard Deviation

Suppose a hypothetical model estimates a random-effects standard deviation using $q=\log(\sigma)$, and gives $\hat{q}=-0.2231$ with $SE(\hat{q})=0.1200$. The estimated standard deviation is $\hat{\sigma}=\exp(-0.2231)\approx0.8000$. To express its approximate standard error on the natural scale, use $g(q)=\exp(q)$, whose derivative is $g'(q)=\exp(q)$. Therefore,

SE(\hat{\sigma}) \approx |g'(\hat{q})|SE(\hat{q}) = 0.8000(0.1200) = 0.09600.

The approximate standard error of the estimated random-effects standard deviation is $0.09600$ in the units of the longitudinal outcome. It is not $0.1200$: that value measures uncertainty in the log standard deviation, which is dimensionless and on a different scale. This delta-method result is local; it need not describe the full uncertainty well when the estimate is near a boundary or the distribution is strongly asymmetric.

Diagnostics and Failure Modes

A numerical inverse is not, by itself, evidence that the standard errors are trustworthy. Before reporting them, check the fit and the curvature calculation. The following checks are especially relevant when the likelihood involves random-effects integration.

Check the maximum. The score should be close to zero at the reported estimate. A noticeably nonzero score suggests optimization has not located a stationary point, so the local quadratic approximation may be misleading.
Check curvature. The observed information should be symmetric up to numerical tolerance and positive definite at a well-behaved interior maximum. A negative eigenvalue can indicate a saddle point, a boundary complication, or inaccurate derivatives. A zero or very small eigenvalue signals a flat or poorly identified direction.
Check numerical stability. Recalculate derivatives and standard errors with tighter numerical tolerances or a more accurate integration calculation. Large changes suggest that numerical error is affecting the Hessian or its inverse.
Check scale and parameterization. Parameters measured in very different units can make a matrix poorly conditioned numerically. Rescaling parameters can help computation, but any reported standard errors must be returned to the intended parameter scale.
Check the full covariance. Large estimated correlations or very large standard errors can reveal weak separation between parameters, including longitudinal and survival components. Such results call for investigation of identifiability and model fit rather than automatic acceptance of a software output.

Worked Example: Weak Identification Despite an Invertible Matrix

Consider a hypothetical pair of dimensionless parameters, such as a longitudinal coefficient and a survival association coefficient. Their observed information is

J= \begin{pmatrix} 4.0 & 3.9\\ 3.9 & 4.0 \end{pmatrix}.

Its determinant is $4.0(4.0)-3.9(3.9)=0.79$, so

J^{-1} = \frac{1}{0.79} \begin{pmatrix} 4.0 & -3.9\\ -3.9 & 4.0 \end{pmatrix} = \begin{pmatrix} 5.063 & -4.937\\ -4.937 & 5.063 \end{pmatrix}.

Each standard error is $\sqrt{5.063}=2.250$, and the covariance-implied correlation is $-4.937/5.063=-0.975$. The eigenvalues of $J$ are $0.1$ and $7.9$: the small eigenvalue indicates a nearly flat direction even though the matrix is invertible. For hypothetical estimates $\hat{\theta}_1=0.20$ and $\hat{\theta}_2=0.30$, the Wald intervals would be approximately $[-4.21,4.61]$ and $[-4.11,4.71]$. These very wide intervals warn that the data provide limited precision for the separate parameters. Reporting the inverse without investigating the weak direction would conceal the key feature of this fit.

How This Approximation Compares With Alternatives

The observed-information calculation uses the negative Hessian at the fitted parameter value. Expected information instead uses the expectation of the negative Hessian under the model. Under regular conditions and a correctly specified model, these approaches are asymptotically related, but their finite-sample values can differ. Expected information may be more stable in some settings, while observed information reflects the curvature at the data actually observed. Neither choice corrects an inaccurate integrated likelihood.

A sandwich covariance estimator is another alternative when the likelihood may be misspecified but independent subject-level contributions are available. It combines an estimate of likelihood curvature with the empirical variability of subject-level scores. It targets a different robustness goal; it does not repair poor identification, an unsuitable censoring assumption, or numerical integration error. Resampling methods can also estimate variability, but they require repeated fitting and can be unstable when some resampled fits fail or are weakly identified. For a well-specified, well-behaved joint model, inverse observed information is a direct and commonly used local approximation.

The approximation has limitations. It can be unreliable with small samples, strong likelihood asymmetry, parameters near boundaries, or weak identification. For a variance component at zero, for example, the usual interior, symmetric normal approximation may not apply. A Wald interval can also extend beyond a parameter’s possible range. In such cases, do not interpret a standard error as a complete description of uncertainty; consider whether a method designed for the particular boundary or asymmetry is appropriate.

Key takeaway: Compute the inverse of the observed information for the full parameter vector, then take square roots of its diagonal entries. Interpret each standard error on the parameterization used, transform it when necessary, and assess whether the likelihood’s local curvature is both numerically reliable and sufficiently informative.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

If the first diagonal entry of $J^{-1}$ is $0.0400$, what is the approximate standard error of the first parameter?
Why is the inverse of a diagonal information entry generally not the correct variance when parameters are correlated?
A log-scale estimate is $\hat{q}=0$ with $SE(\hat{q})=0.15$, and $\sigma=\exp(q)$. Use the delta method to find the approximate standard error of $\hat{\sigma}$.
What does a very small eigenvalue of the observed information suggest, and why should it prompt further checks?

Solutions. (1) The standard error is $\sqrt{0.0400}=0.2000$. (2) The marginal variances are the diagonal entries of the full inverse matrix; off-diagonal information changes the inverse when parameters are correlated. (3) Since $g'(q)=\exp(q)$, $SE(\hat{\sigma})\approx\exp(0)(0.15)=0.1500$, on the natural scale of $\sigma$. (4) It indicates a nearly flat likelihood direction and possible weak identification. The resulting inverse can be unstable, so check the fit, numerical accuracy, and parameter dependence before interpreting the standard errors.