Here we describe the details of the calculations for the Kenward-Roger degrees of freedom and the adjusted covariance matrix of the coefficients.
Model definition
The model definition is the same as what we have in Details of the model fitting in
mmrm. We are using the same notations.
Linear model
For each subject \(i\) we observe a vector \[ Y_i = (y_{i1}, \dotsc, y_{im_i})^\top \in \mathbb{R}^{m_i} \] and given a design matrix \[ X_i \in \mathbb{R}^{m_i \times p} \] and a corresponding coefficient vector \(\beta \in \mathbb{R}^{p}\) we assume that the observations are multivariate normal distributed: \[ Y_i \sim N(X_i\beta, \Sigma_i) \] where the covariance matrix \(\Sigma_i \in \mathbb{R}^{m_i \times m_i}\) is derived by subsetting the overall covariance matrix \(\Sigma \in \mathbb{R}^{m \times m}\) appropriately by \[ \Sigma_i = G_i^{-1/2} S_i^\top \Sigma S_i G_i^{-1/2} \] where the subsetting matrix \(S_i \in \{0, 1\}^{m \times m_i}\) contains in each of its \(m_i\) columns a single 1 indicating which overall time point is matching the \(h\)-th observed time point \(t_{ih}\) of subject \(i\). \(G_i \in \mathbb{R}^{m_i \times m_i}\) is a diagonal matrix with strictly positive weights on its diagonal.
Conditional on the design matrices \(X_i\), the coefficient vector \(\beta\) and the covariance matrix \(\Sigma\) we assume that the observations are independent between the subjects.
We can write the linear model for all subjects together as \[ Y = X\beta + \epsilon \] where \(Y \in \mathbb{R}^N\) combines all subject specific observations vectors \(Y_i\) such that we have in total \(N = \sum_{i = 1}^{n}{m_i}\) observations, \(X \in \mathbb{R}^{N \times p}\) combines all subject specific design matrices and \(\epsilon \in \mathbb{R}^N\) has a multivariate normal distribution \[ \epsilon \sim N(0, \Omega) \] where \(\Omega \in \mathbb{R}^{N \times N}\) is block-diagonal containing the subject specific \(\Sigma_i\) covariance matrices on the diagonal and 0 in the remaining entries.
Mathematical Details of Kenward-Roger method
The mathematical derivation of the Kenward-Roger method is based on
the Taylor expansion of the obtained covariance matrix of \(\hat\beta\) to get a more accurate estimate
for it. All these derivations are based on the restricted maximum
likelihood. Following the same notation, the
covariance matrix, \(\Omega\) can be
represented as a function of covariance matrix parameters \(\theta = (\theta_1, \dotsc,
\theta_k)^\top\), i.e. \(\Omega(\theta)\). Here after model fitting
with mmrm, we obtain the estimate \(\hat\beta =
\Phi(\hat\theta)X^\top\Omega(\hat\theta)^{-1}Y\), where \(\Phi(\theta) = \left\{X^\top \Omega(\theta)^{-1}
X\right\} ^{-1}\) is the asymptotic covariance matrix of \(\hat\beta\). However, Kackar and Harville (1984) suggests that
although the \(\hat\beta\) is unbiased
for \(\beta\), the covariance matrix,
\(\hat\Phi = \left\{X^\top \hat\Omega^{-1}
X\right\}^{-1}\) can be biased. They showed that the variability
of \(\hat\beta\) can be partitioned
into two components,
\[ \Phi_A = \Phi + \Lambda \]
where \(\Phi\) is the variance-covariance matrix of the asymptotic distribution of \(\hat\beta\) as \(n\rightarrow \infty\) as defined above, and \(\Lambda\) represents the amount to which the asymptotic variance-covariance matrix underestimates \(\Phi_A\).
Based on a Taylor series expansion around \(\theta\), \(\Lambda\) can be approximated by
\[ \Lambda \simeq \Phi \left\{\sum_{h=1}^k{\sum_{j=1}^k{W_{hj}(Q_{hj} - P_h \Phi P_j)} }\right\} \Phi \] where \[ P_h = X^\top \frac{\partial{\Omega^{-1}}}{\partial \theta_h} X \] \[ Q_{hj} = X^\top \frac{\partial{\Omega^{-1}}}{\partial \theta_h} \Omega \frac{\partial{\Omega^{-1}}}{\partial \theta_j} X \]
\(W\) is the variance-covariance matrix of \(\hat\theta\). It is estimated by the inverse of the observed information, i.e. the inverse of the Hessian matrix of the negative restricted log-likelihood function of \(\theta\), evaluated at the estimate \(\hat\theta\).
Again, based on a Taylor series expansion about \(\theta\), Kenward and Roger (1997) show that \[ \hat\Phi \simeq \Phi + \sum_{h=1}^k{(\hat\theta_h - \theta_h)\frac{\partial{\Phi}}{\partial{\theta_h}}} + \frac{1}{2} \sum_{h=1}^k{\sum_{j=1}^k{(\hat\theta_h - \theta_h)(\hat\theta_j - \theta_j)\frac{\partial^2{\Phi}}{\partial{\theta_h}\partial{\theta_j}}}} \] Ignoring the possible bias in \(\hat\theta\), \[ E(\hat\Phi) \simeq \Phi + \frac{1}{2} \sum_{h=1}^k{\sum_{j=1}^k{W_{hj}\frac{\partial^2{\Phi}}{\partial{\theta_h}\partial{\theta_j}}}} \] Using previously defined notations, this can be further written as \[ \frac{\partial^2{\Phi}}{\partial{\theta_h}\partial{\theta_j}} = \Phi (P_h \Phi P_j + P_j \Phi P_h - Q_{hj} - Q_{jh} + R_{hj}) \Phi \] where \[ R_{hj} = X^\top\Omega^{-1}\frac{\partial^2\Omega}{\partial{\theta_h}\partial{\theta_j}} \Omega^{-1} X \]
Combining this with \(\Phi_A = \Phi + \Lambda\), we obtain the adjusted estimate
\[ \hat\Phi_A = \hat\Phi + 2\hat\Phi \left\{\sum_{h=1}^k{\sum_{j=1}^k{W_{hj}(Q_{hj} - P_h \hat\Phi P_j - \frac{1}{4}R_{hj})} }\right\} \hat\Phi \]
where \(\Omega(\hat\theta)\) replaces \(\Omega(\theta)\) in the right-hand side. In the remainder of this vignette, all of \(\Phi\), \(\Phi_A\), \(P_h\), \(Q_{hj}\), \(R_{hj}\) and \(W\) are evaluated at \(\hat\theta\), and we drop the hats for brevity.
Please note that, if we ignore \(R_{hj}\), the second-order derivatives, we will get a different estimate of adjusted covariance matrix, and we call this the linear Kenward-Roger approximation.
Special Considerations for mmrm models
In mmrm models, \(\Omega\) is a block-diagonal matrix, hence we can calculate \(P\), \(Q\) and \(R\) for each of the \(n\) subjects and add them up.
\[ P_h = \sum_{i=1}^{n}{P_{ih}} = \sum_{i=1}^{n}{X_i^\top \frac{\partial{\Sigma_i^{-1}}}{\partial \theta_h} X_i} \]
\[ Q_{hj} = \sum_{i=1}^{n}{Q_{ihj}} = \sum_{i=1}^{n}{X_i^\top \frac{\partial{\Sigma_i^{-1}}}{\partial \theta_h} \Sigma_i \frac{\partial{\Sigma_i^{-1}}}{\partial \theta_j} X_i} \]
\[ R_{hj} = \sum_{i=1}^{n}{R_{ihj}} = \sum_{i=1}^{n}{X_i^\top\Sigma_i^{-1}\frac{\partial^2\Sigma_i}{\partial{\theta_h}\partial{\theta_j}} \Sigma_i^{-1} X_i} \]
Derivative of the overall covariance matrix \(\Sigma\)
The derivative of the overall covariance matrix \(\Sigma\) with respect to the variance parameters can be calculated through the derivatives of the Cholesky factor, and hence obtained through automatic differentiation, following matrix identities calculations. \[ \frac{\partial{\Sigma}}{\partial{\theta_h}} = \frac{\partial{LL^\top}}{\partial{\theta_h}} = \frac{\partial{L}}{\partial{\theta_h}}L^\top + L\frac{\partial{L^\top}}{\partial{\theta_h}} \]
\[ \frac{\partial^2{\Sigma}}{\partial{\theta_h}\partial{\theta_j}} = \frac{\partial^2{L}}{\partial{\theta_h}\partial{\theta_j}}L^\top + L\frac{\partial^2{L^\top}}{\partial{\theta_h}\partial{\theta_j}} + \frac{\partial{L}}{\partial{\theta_h}}\frac{\partial{L^\top}}{\partial{\theta_j}} + \frac{\partial{L}}{\partial{\theta_j}}\frac{\partial{L^\top}}{\partial{\theta_h}} \]
Derivative of the \(\Sigma^{-1}\)
The derivatives of \(\Sigma^{-1}\) can be calculated through
\[ \frac{\partial{\Sigma\Sigma^{-1}}}{\partial{\theta_h}}\\ = \frac{\partial{\Sigma}}{\partial{\theta_h}}\Sigma^{-1} + \Sigma\frac{\partial{\Sigma^{-1}}}{\partial{\theta_h}} \\ = 0 \] \[ \frac{\partial{\Sigma^{-1}}}{\partial{\theta_h}} = - \Sigma^{-1} \frac{\partial{\Sigma}}{\partial{\theta_h}}\Sigma^{-1} \]
Subjects with missed visits
For unit weights, if a subject does not have all visits, the corresponding covariance matrix can be represented as \[ \Sigma_i = S_i^\top \Sigma S_i \]
and the derivatives can be obtained through
\[ \frac{\partial{\Sigma_i}}{\partial{\theta_h}} = S_i^\top \frac{\partial{\Sigma}}{\partial{\theta_h}} S_i \]
\[ \frac{\partial^2{\Sigma_i}}{\partial{\theta_h}\partial{\theta_j}} = S_i^\top \frac{\partial^2{\Sigma}}{\partial{\theta_h}\partial{\theta_j}} S_i \]
The derivative of the \(\Sigma_i^{-1}\), \(\frac{\partial\Sigma_i^{-1}}{\partial{\theta_h}}\) can be calculated through \(\Sigma_i^{-1}\) and \(\frac{\partial{\Sigma_i}}{\partial{\theta_h}}\) using the above.
Scenario under group specific covariance estimates
For unit weights in grouped mmrm models, the covariance
matrix for subject \(i\) of group \(g(i)\) can be written as \[
\Sigma_i = S_i^\top \Sigma_{g(i)} S_i.
\] Assume there are \(B\)
groups. Each group \(b\) has its own
covariance matrix \(\Sigma_b\) with its
own parameters, so the number of covariance parameters is multiplied by
\(B\). Let \(\mathcal G_b\) denote the indices of the
covariance parameters of group \(b\). A
parameter \(\theta_h\) with \(h \in \mathcal G_b\) only affects \(\Sigma_b\), hence \(\partial\Sigma_i/\partial\theta_h = 0\) and
\(\partial\Sigma_i^{-1}/\partial\theta_h =
0\) for all subjects \(i\) with
\(g(i) \neq b\). Therefore, for \(h \in \mathcal G_b\), only the subjects of
group \(b\) contribute: \[
P_h = \sum_{g(i) = b}{P_{ih}} = \sum_{g(i) = b}{X_i^\top
\frac{\partial{\Sigma_i^{-1}}}{\partial \theta_h} X_i}.
\] Each \(P_h\) is a \(p \times p\) matrix as before. For \(Q_{hj}\) and \(R_{hj}\) both parameters need to affect the
same subject, so they are zero matrices if \(\theta_h\) and \(\theta_j\) belong to different groups. For
\(h, j \in \mathcal G_b\) we have \[
Q_{hj} = \sum_{g(i) = b}{Q_{ihj}} = \sum_{g(i) = b}{X_i^\top
\frac{\partial{\Sigma_i^{-1}}}{\partial \theta_h} \Sigma_i
\frac{\partial{\Sigma_i^{-1}}}{\partial \theta_j} X_i}
\] \[
R_{hj} = \sum_{g(i) = b}{R_{ihj}} = \sum_{g(i) =
b}{X_i^\top\Sigma_i^{-1}\frac{\partial^2\Sigma_i}{\partial{\theta_h}\partial{\theta_j}}
\Sigma_i^{-1} X_i}.
\] Hence only the \(Q_{hj}\) and
\(R_{hj}\) within each group need to be
calculated and stored. Note that this does not hold for the products
\(P_h \Phi P_j\), which are in general
not zero for parameters from different groups. Therefore the cross-group
terms \(W_{hj} P_h \Phi P_j\) still
contribute to the adjusted covariance matrix \(\Phi_A\).
Scenario under weighted mmrm
Under a weighted mmrm model, the covariance matrix for subject \(i\) can be represented as
\[ \Sigma_i = G_i^{-1/2} S_i^\top \Sigma S_i G_i^{-1/2} \]
where \(G_i\) is a diagonal matrix of fixed weights. Write \(\bar\Sigma_i = S_i^\top\Sigma S_i\) (or use \(\Sigma_{g(i)}\) for grouped models). Since the weights do not depend on \(\theta\), the covariance derivatives are \[ \frac{\partial\Sigma_i}{\partial\theta_h} = G_i^{-1/2}\frac{\partial\bar\Sigma_i}{\partial\theta_h}G_i^{-1/2}, \qquad \frac{\partial^2\Sigma_i}{\partial\theta_h\partial\theta_j} = G_i^{-1/2}\frac{\partial^2\bar\Sigma_i}{\partial\theta_h\partial\theta_j}G_i^{-1/2}. \] The inverse covariance instead satisfies \(\Sigma_i^{-1} = G_i^{1/2}\bar\Sigma_i^{-1}G_i^{1/2}\), so \[ \frac{\partial\Sigma_i^{-1}}{\partial\theta_h} = G_i^{1/2}\frac{\partial\bar\Sigma_i^{-1}}{\partial\theta_h}G_i^{1/2}. \] These weighted matrices are inserted into the same definitions of \(P_h\), \(Q_{hj}\) and \(R_{hj}\). See subject contributions with weights below for how the weights then enter the subject contributions to \(P_h\) and \(Q_{hj}\).
Contracted linear covariance adjustment
For vcov = "Kenward-Roger-Linear" the \(R_{hj}\) terms are omitted, and the
adjusted covariance matrix above only requires two \(p\times p\) sums:
\[ S_Q = \sum_{h=1}^k\sum_{j=1}^k W_{hj}Q_{hj}, \qquad S_P = \sum_{h=1}^k\sum_{j=1}^k W_{hj}P_h\Phi P_j, \qquad \Phi_A = \Phi + 2\Phi(S_Q-S_P)\Phi. \]
This is the linear Kenward-Roger formula with the outer
multiplications by \(\Phi\) moved
outside the sums. Instead of storing all \(Q_{hj}\) matrices,
h_get_kr_comp() returns the \(P_h\) matrices together with \(S_Q\) (so kr_comp$Q and
kr_comp$R are NULL and
kr_comp$S_Q is a \(p\times
p\) matrix), and h_var_adj_contracted() calculates
\(S_P\) and \(\Phi_A\). The result is the same as the
pairwise calculation up to floating-point error. Full Kenward-Roger
continues to use the pairwise \(Q_{hj}\) and \(R_{hj}\).
Subject contributions with weights
As in the weighted model, write the subject covariance matrix as \(\Sigma_i = G_i^{-1/2}\bar\Sigma_iG_i^{-1/2}\), where \(\bar\Sigma_i = S_i^\top\Sigma_{g(i)}S_i\) is the unweighted covariance matrix of the observed visits of subject \(i\) in group \(g(i)\). For spatial covariance structures, \(\bar\Sigma_i\) is instead calculated from the distances between the subject’s observations. Let
\[ X_i^\ast = G_i^{1/2}X_i, \qquad D_{ih} = \frac{\partial\bar\Sigma_i^{-1}}{\partial\theta_h} = -\bar\Sigma_i^{-1}\frac{\partial\bar\Sigma_i}{\partial\theta_h}\bar\Sigma_i^{-1}. \]
The weights do not depend on \(\theta\), so \(\partial\Sigma_i^{-1}/\partial\theta_h = G_i^{1/2}D_{ih}G_i^{1/2}\). Inserting this into the subject contributions cancels the inner weight factors in \(Q_{ihj}\), and the weights enter exactly once through \(X_i^\ast\):
\[ P_{ih} = (X_i^\ast)^\top D_{ih}X_i^\ast, \qquad Q_{ihj} = (X_i^\ast)^\top D_{ih}\bar\Sigma_i D_{ij}X_i^\ast. \]
Contraction of \(Q\)
As in the grouped model, \(\mathcal G_b\) denotes the indices of the covariance parameters of group \(b\), so that \(\theta_h\) with \(h \notin \mathcal G_{g(i)}\) does not affect \(\bar\Sigma_i\) and \(D_{ih} = 0\). Hence only the diagonal blocks \(W^{(b)} = (W_{hj})_{h,j\in\mathcal G_b}\) enter \(S_Q\):
\[ S_Q = \sum_{i=1}^{n}\sum_{h,j\in\mathcal G_{g(i)}} W_{hj}\, (X_i^\ast)^\top D_{ih}\bar\Sigma_i D_{ij}X_i^\ast. \]
Each block is factored as follows, with the rows of \(U_b\) indexed by the parameter labels in \(\mathcal G_b\):
\[ W^{(b)} = U_b\operatorname{diag}(s_b)U_b^\top, \qquad \text{i.e.} \quad W_{hj} = \sum_{\ell}s_{b\ell}(U_b)_{h\ell}(U_b)_{j\ell} \quad\text{for } h,j\in\mathcal G_b. \]
If \(W^{(b)}\) is positive definite, \(U_b\) is its lower Cholesky factor and all \(s_{b\ell} = 1\). Otherwise, with the eigendecomposition \(W^{(b)} = V\operatorname{diag}(\omega)V^\top\), the columns are \((U_b)_{\cdot\ell} = \sqrt{|\omega_\ell|}\,V_{\cdot\ell}\) and \(s_{b\ell} = \operatorname{sign}(\omega_\ell)\), with \(s_{b\ell} = 1\) for \(\omega_\ell = 0\). Every term is retained, including negative eigenvalues, so this is an exact rearrangement and not a low-rank approximation.
Substituting the factorization for subject \(i\) in group \(b = g(i)\) and regrouping the double sum gives
\[ \sum_{h,j\in\mathcal G_b} W_{hj}D_{ih}\bar\Sigma_iD_{ij} = \sum_\ell s_{b\ell} \Big(\sum_{h\in\mathcal G_b}(U_b)_{h\ell}D_{ih}\Big)\bar\Sigma_i \Big(\sum_{j\in\mathcal G_b}(U_b)_{j\ell}D_{ij}\Big) = \sum_\ell s_{b\ell}\,\Delta_{i\ell}\bar\Sigma_i\Delta_{i\ell}, \]
where the derivative combinations
\[ \Delta_{i\ell} = \sum_{h\in\mathcal G_{g(i)}}(U_{g(i)})_{h\ell}D_{ih} \]
are symmetric because each \(D_{ih}\) is symmetric. Therefore
\[ S_Q = \sum_{i=1}^{n}\sum_\ell s_{g(i)\ell}\, (\Delta_{i\ell}X_i^\ast)^\top\,\bar\Sigma_i\,(\Delta_{i\ell}X_i^\ast). \]
For a group with \(|\mathcal G_b|\) covariance parameters, each subject now needs \(|\mathcal G_b|\) such products instead of \(|\mathcal G_b|^2\). For non-spatial covariance structures, \(\bar\Sigma_i\), \(D_{ih}\) and \(\Delta_{i\ell}\) only depend on the group and the observed visits, so they are calculated once per group and visit pattern. For spatial covariance structures they are calculated for each subject, because the distances can differ between subjects.
\(W\) is the inverse of the Hessian
matrix and therefore symmetric up to rounding errors. The factorization
uses the symmetric part \((W^{(b)} +
W^{(b)\top})/2\), which yields the symmetric part of the pairwise
sum because \(Q_{jh} = Q_{hj}^\top\).
If \(W\) contains non-finite values,
which mmrm() reports as a convergence problem, \(S_Q\) and hence \(\Phi_A\) are non-finite, as with the
pairwise calculation.
Contraction of \(P\)
In contrast to \(Q_{hj}\), the product \(P_h\Phi P_j\) is in general not zero for parameters \(h\) and \(j\) from different groups. In addition, \(W\) is in general not block diagonal, because the restricted likelihood couples the groups through the common coefficients \(\beta\). Therefore \(S_P\) uses the full \(W\), including the cross-group entries:
\[ S_P = \sum_{h=1}^k P_h\Phi\overline P_h, \qquad \overline P_h = \sum_{j=1}^k W_{hj}P_j. \]
All \(\overline P_h\) are obtained with one matrix product: if the \(h\)th column of the \(p^2\times k\) matrix \(\mathcal P\) is \(\operatorname{vec}(P_h)\), then the \(h\)th column of \(\mathcal P W^\top\) is \(\operatorname{vec}(\overline P_h)\). The \(P_h\) and \(W\) are also used unchanged for the degrees of freedom.
Inference
Suppose we are testing \(C\beta = C\beta_0\), where \(\beta_0\) specifies the null hypothesis and \(C \in \mathbb{R}^{c\times p}\) has full row rank. We use the following F-statistic \[ F = \frac{1}{c} (C(\hat\beta - \beta_0))^\top (C \Phi_A C^\top)^{-1} C(\hat\beta - \beta_0) \] and \[ F^* = \lambda F \] is approximated by an \(F_{c,\nu}\) distribution.
\(\lambda\) and \(\nu\) can be calculated through
\[ M = C^\top (C \Phi C^\top)^{-1} C \]
\[ A_1 = \sum_{h=1}^k{\sum_{j=1}^k{W_{hj} \operatorname{tr}(M \Phi P_h \Phi) \operatorname{tr}(M \Phi P_j \Phi)}} \]
\[ A_2 = \sum_{h=1}^k{\sum_{j=1}^k{W_{hj} \operatorname{tr}(M \Phi P_h \Phi M \Phi P_j \Phi)}} \]
\[ \mathcal B = \frac{1}{2c}(A_1 + 6A_2) \]
\[ g = \frac{(c+1)A_1 - (c+4)A_2}{(c+2)A_2} \]
\[ c_1 = \frac{g}{3c+2(1-g)} \]
\[ c_2 = \frac{c-g}{3c+2(1-g)} \]
\[ c_3 = \frac{c+2-g}{3c+2(1-g)} \] \[E^*={\left\{1-\frac{A_2}{c}\right\}}^{-1}\] \[V^*=\frac{2}{c}{\left\{\frac{1+c_1 \mathcal B}{(1-c_2 \mathcal B)^2(1-c_3 \mathcal B)}\right\}}\]
\[\eta = \frac{V^{*}}{2(E^*)^2}\]
\[\nu = 4 + \frac{c+2}{c\eta - 1}\] \[\lambda = \frac{\nu}{E^*(\nu-2)}\]
Degrees of freedom in contrast space
h_kr_df() evaluates the same \(A_1\) and \(A_2\) in the smaller contrast space, for
both full and linear Kenward-Roger and all covariance structures. Let
\(H_C = C\Phi C^\top = TT^\top\), where
\(T\) is the lower triangular Cholesky
factor, and set
\[ Z = T^{-1}C\Phi, \qquad K_h = ZP_hZ^\top. \]
Triangular solves obtain \(Z\) without explicitly inverting \(H_C\). Hypotheses whose \(H_C\) is numerically singular or not positive definite are rejected. Each \(K_h\) is a symmetric \(c\times c\) matrix. Cyclic invariance of trace gives
\[ A_1 = t^\top Wt, \qquad t_h = \operatorname{tr}(K_h), \qquad A_2 = \sum_{h,j} W_{hj}\langle K_h,K_j\rangle_F, \]
where \(\langle\cdot,\cdot\rangle_F\) is the Frobenius inner product. With \(\mathcal{K}\) containing \(\operatorname{vec}(K_h)\) as its \(h\)th column, the implementation calculates \(A_2 = \operatorname{sum}((\mathcal{K}W)\odot\mathcal{K})\). This replaces the pairwise products of \(p\times p\) matrices by a dense contraction with \(c^2\) rows. The remaining scalar formulas above are unchanged. The full covariance \(W\) is retained, including cross-group entries in grouped models.
One-dimensional shortcut
For a single row contrast \(C = l^\top\) with \(l \in \mathbb{R}^p\), define
\[ a_h = \frac{l^\top\Phi P_h\Phi l}{l^\top\Phi l}. \]
Then \(A_1 = A_2 = a^\top Wa\), and the formulas simplify to
\[ \nu = \frac{2}{a^\top Wa}, \qquad \lambda = 1. \]
The implementation uses this shortcut directly, avoiding Cholesky normalization and the general moment formulas. It applies to any scalar contrast with positive variance \(l^\top\Phi l\), including a combination of coefficients; other contrasts are rejected. If \(a^\top Wa = 0\), the limit \(\nu = \infty\) is returned. Since \(\partial\Phi/\partial\theta_h = -\Phi P_h\Phi\), this is exactly the scalar Satterthwaite degrees of freedom based on the unadjusted variance \(l^\top\Phi l\). The standard error and test statistic still use \(\Phi_A\) for both full and linear KR. Consequently, their tests and confidence intervals can differ from Satterthwaite even when their scalar degrees of freedom agree. Multi-dimensional KR retains its own degrees of freedom and F scaling.
Parameterization methods and Kenward-Roger
While the Kenward-Roger adjusted covariance matrix is adopting a Taylor series to approximate the true value, the choices of parameterization can change the result. In a simple example of unstructured covariance structure, in our current approach, where the parameters are log diagonal entries and scaled off-diagonal entries of the Cholesky factor of \(\Sigma\) (see parameterization), the second-order derivatives of \(\Sigma\) over our parameters, are non-zero matrices. However, if we use the elements of \(\Sigma\) as our parameters, then the second-order derivatives are zero matrices. However, the differences are usually small. If you would like to match SAS results for the unstructured covariance model, you can use the linear Kenward-Roger approximation.
Implementations in mmrm
In package mmrm, we have implemented Kenward-Roger
calculations based on the previous sections. For non-spatial covariance
structures, the first-order and second-order derivatives are obtained by
automatic differentiation of the Cholesky factor. For spatial covariance
structures, we derive the exact results below and evaluate them
analytically on demand.
For vcov = "Kenward-Roger-Linear", only the first
derivatives of the subject covariance matrices and their inverses are
calculated and cached. For non-spatial covariance structures, automatic
differentiation is applied once to the Cholesky factor, and the nested
differentiation for the second derivatives is skipped. Missing-visit
patterns cache only the required first derivatives, and for spatial
covariance structures the second derivatives are not evaluated either.
Full Kenward-Roger still computes the second derivatives and the
parameter-pair \(Q_{hj}\) and \(R_{hj}\) components.
Spatial Exponential Derivatives
For spatial exponential covariance, \(\sigma\) denotes a variance (the diagonal of \(\Sigma\)), and \(\rho\) is the correlation at unit distance. For one covariance group, we have
\[\theta = (\theta_1,\theta_2)^\top\] \[\sigma = e^{\theta_1}\] \[\rho = \frac{e^{\theta_2}}{1 + e^{\theta_2}}\]
\[\Sigma_{ij} = \sigma \rho^{d_{ij}}\] where \(d_{ij}\) is the distance between time point \(i\) and time point \(j\).
So the first-order derivatives can be written as:
\[ \frac{\partial{\Sigma_{ij}}}{\partial\theta_1} = \frac{\partial\sigma}{\partial\theta_1} \rho^{d_{ij}}\\ = e^{\theta_1}\rho^{d_{ij}} \\ = \Sigma_{ij} \]
\[ \frac{\partial{\Sigma_{ij}}}{\partial\theta_2} = \sigma\frac{\partial{\rho^{d_{ij}}}}{\partial\theta_2} \\ = \sigma\rho^{d_{ij}-1}{d_{ij}}\frac{\partial\rho}{\partial\theta_2}\\ = \sigma\rho^{d_{ij}-1}{d_{ij}}\rho(1-\rho) \\ = \sigma \rho^{d_{ij}} {d_{ij}} (1-\rho) \]
Second-order derivatives can be written as:
\[ \frac{\partial^2{\Sigma_{ij}}}{\partial\theta_1\partial\theta_1}\\ = \frac{\partial\Sigma_{ij}}{\partial\theta_1}\\ = \Sigma_{ij} \]
\[ \frac{\partial^2{\Sigma_{ij}}}{\partial\theta_1\partial\theta_2} = \frac{\partial^2{\Sigma_{ij}}}{\partial\theta_2\partial\theta_1} \\ = \frac{\partial\Sigma_{ij}}{\partial\theta_2}\\ = \sigma\rho^{d_{ij}-1}{d_{ij}}\rho(1-\rho)\\ = \sigma\rho^{d_{ij}}{d_{ij}}(1-\rho) \]
\[ \frac{\partial^2{\Sigma_{ij}}}{\partial\theta_2\partial\theta_2}\\ = \frac{\partial{\sigma\rho^{d_{ij}}{d_{ij}}(1-\rho)}}{\partial\theta_2}\\ = \sigma\rho^{d_{ij}}{d_{ij}}(1-\rho)(d_{ij} (1-\rho) - \rho) \]
Spatial Gaussian Derivatives
For spatial Gaussian covariance structure, we use the same parameterization as spatial exponential, with the squared distance in the exponent:
\[\theta = (\theta_1,\theta_2)^\top\] \[\sigma = e^{\theta_1}\] \[\rho = \frac{e^{\theta_2}}{1 + e^{\theta_2}}\]
\[\Sigma_{ij} = \sigma \rho^{d_{ij}^2}\] where \(d_{ij}\) is the distance between time point \(i\) and time point \(j\).
The derivations parallel the spatial exponential case, with \(d_{ij}\) replaced by \(d_{ij}^2\) where it appears as an exponent of \(\rho\) (and as a multiplicative factor brought down by the chain rule). The first-order derivatives are:
\[ \frac{\partial{\Sigma_{ij}}}{\partial\theta_1} = \frac{\partial\sigma}{\partial\theta_1} \rho^{d_{ij}^2}\\ = e^{\theta_1}\rho^{d_{ij}^2} \\ = \Sigma_{ij} \]
\[ \frac{\partial{\Sigma_{ij}}}{\partial\theta_2} = \sigma\frac{\partial{\rho^{d_{ij}^2}}}{\partial\theta_2} \\ = \sigma\rho^{d_{ij}^2-1}{d_{ij}^2}\frac{\partial\rho}{\partial\theta_2}\\ = \sigma\rho^{d_{ij}^2-1}{d_{ij}^2}\rho(1-\rho) \\ = \sigma \rho^{d_{ij}^2} {d_{ij}^2} (1-\rho) \]
Second-order derivatives:
\[ \frac{\partial^2{\Sigma_{ij}}}{\partial\theta_1\partial\theta_1}\\ = \frac{\partial\Sigma_{ij}}{\partial\theta_1}\\ = \Sigma_{ij} \]
\[ \frac{\partial^2{\Sigma_{ij}}}{\partial\theta_1\partial\theta_2} = \frac{\partial^2{\Sigma_{ij}}}{\partial\theta_2\partial\theta_1} \\ = \frac{\partial\Sigma_{ij}}{\partial\theta_2}\\ = \sigma\rho^{d_{ij}^2}{d_{ij}^2}(1-\rho) \]
\[ \frac{\partial^2{\Sigma_{ij}}}{\partial\theta_2\partial\theta_2}\\ = \frac{\partial{\sigma\rho^{d_{ij}^2}{d_{ij}^2}(1-\rho)}}{\partial\theta_2}\\ = \sigma\rho^{d_{ij}^2}{d_{ij}^2}(1-\rho)(d_{ij}^2 (1-\rho) - \rho) \]
As discussed in Parameterization methods and Kenward-Roger, the inverse-logit reparameterization of \(\rho\) produces non-zero second-order derivatives where SAS’s natural \((0, 1)\)-scaled parameterization would give zero. The likelihood, \(\beta\) estimates and KR-Linear adjusted standard errors all match SAS to numerical precision; small (typically sub-percent) differences remain on the default Kenward-Roger standard error. To reproduce SAS exactly, use the linear Kenward-Roger approximation.