Here we describe the details of the calculations for the Kenward-Roger degrees of freedom and the adjusted covariance matrix of the coefficients.
The model definition is the same as what we have in Details of the model fitting in
mmrm. We are using the same notations.
For each subject i we observe a vector Yi = (yi1, …, yimi)⊤ ∈ ℝmi and given a design matrix Xi ∈ ℝmi × p and a corresponding coefficient vector β ∈ ℝp we assume that the observations are multivariate normal distributed: Yi ∼ N(Xiβ, Σi) where the covariance matrix Σi ∈ ℝmi × mi is derived by subsetting the overall covariance matrix Σ ∈ ℝm × m appropriately by Σi = Gi−1/2Si⊤ΣSiGi−1/2 where the subsetting matrix Si ∈ {0, 1}m × mi contains in each of its mi columns contains a single 1 indicating which overall time point is matching tih. Gi ∈ ℝ > 0mi × mi is the diagonal weight matrix.
Conditional on the design matrices Xi, the coefficient vector β and the covariance matrix Σ we assume that the observations are independent between the subjects.
We can write the linear model for all subjects together as Y = Xβ + ϵ where Y ∈ ℝN combines all subject specific observations vectors Yi such that we have in total $N = \sum_{i = 1}^{n}{m_i}$ observations, X ∈ ℝN × p combines all subject specific design matrices and ϵ ∈ ℝN has a multivariate normal distribution ϵ ∼ N(0, Ω) where Ω ∈ ℝN × N is block-diagonal containing the subject specific Σi covariance matrices on the diagonal and 0 in the remaining entries.
The mathematical derivation of the Kenward-Roger method is based on
the Taylor expansion of the obtained covariance matrix of β̂ 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, Ω can be
represented as a function of covariance matrix parameters θ = (θ1, …, θk)⊤,
i.e. Ω(θ). Here after
model fitting with mmrm, we obtain the estimate β̂ = Φ(θ̂)X⊤Ω(θ̂)−1Y,
where Φ(θ) = {X⊤Ω(θ)−1X}−1
is the asymptotic covariance matrix of β̂. However, Kackar and Harville (1984) suggests that
although the β̂ is unbiased for
β, the covariance matrix,
Φ̂ = {X⊤Ω̂X}−1
can be biased. They showed that the variability of β̂ can be partitioned into two
components,
ΦA = Φ + Λ
where Φ is the variance-covariance matrix of the asymptotic distribution of β̂ as n → ∞ as defined above, and Λ represents the amount to which the asymptotic variance-covariance matrix underestimates ΦA.
Based on a Taylor series expansion around θ, Λ 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 inverse of the Hessian matrix of the log-likelihood function of θ evaluated at the estimate θ̂, i.e. the observed Fisher Information matrix as a consistent estimator of the variance-covariance matrix of θ̂.
Again, based on a Taylor series expansion about θ, 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 θ̂, $$ 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 $$
substituting Φ and Λ back to the Φ̂A, we have
$$ \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 Ω(θ̂) replaces Ω(θ) in the right-hand side.
Please note that, if we ignore Rhj, the second-order derivatives, we will get a different estimate of adjusted covariance matrix, and we call this the linear Kenward-Roger approximation.
In mmrm models, Ω is a block-diagonal matrix, hence we can calculate P, Q and R for each subject 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} $$
The derivative of the overall covariance matrix Σ 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^T}}{\partial{\theta_j}} + \frac{\partial{L}}{\partial{\theta_j}}\frac{\partial{L^\top}}{\partial{\theta_h}} $$
The derivatives of Σ−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} $$
If a subject do not have all visits, the corresponding covariance matrix can be represented as Σi = Si⊤ΣSi
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 Σi−1, $\frac{\partial\Sigma_i^{-1}}{\partial{\theta_h}}$ can be calculated through Σi−1 and $\frac{\partial{\Sigma_i}}{\partial{\theta_h}}$ using the above.
When fitting 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, the number of parameters is increased by B times. With the fact that for each
group, the corresponding θ
will not affect other parts, we will have block-diagonal P, Q and R matrices, where the blocks are
given by:
$$ P_h = \begin{pmatrix} P_{h, 1} & \dots & P_{h, B} \\ \end{pmatrix} $$
$$ Q_{hj} = \begin{pmatrix} Q_{hj, 1} & 0 & \dots & \dots & 0 \\ 0 & Q_{hj, 2} & 0 & \dots & 0\\ \vdots & & \ddots & & \vdots \\ \vdots & & & \ddots & \vdots \\ 0 & \dots & \dots & 0 & Q_{hj, B} \end{pmatrix} $$
$$ R_{hj} = \begin{pmatrix} R_{hj, 1} & 0 & \dots & \dots & 0 \\ 0 & R_{hj, 2} & 0 & \dots & 0\\ \vdots & & \ddots & & \vdots \\ \vdots & & & \ddots & \vdots \\ 0 & \dots & \dots & 0 & R_{hj, B} \end{pmatrix} $$
Use Pj, b to denote the block diagonal part for group b, we have $$ P_{h,b} = \sum_{g(i) = b}{P_{ih}} = \sum_{g(i) = b}{X_i^\top \frac{\partial{\Sigma_i^{-1}}}{\partial \theta_h} X_i} $$
$$ Q_{hj,b} = \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,b} = \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} $$
Similarly for R.
Under weights mmrm model, the covariance matrix for subject i, can be represented as
Σi = Gi−1/2Si⊤ΣSiGi−1/2
Where Gi is a diagonal matrix of the weights. Then, when deriving P, Q and R, there are no mathematical differences as they are constant, and having Gi in addition to Si does not change the algorithms and we can simply multiply the formulas with Gi−1/2, similarly as above for the subsetting matrix.
Suppose we are testing the linear combination of β, Cβ with C ∈ ℝc × p, we can use the following F-statistic $$ F = \frac{1}{c} (\hat\beta - \beta)^\top C (C^\top \hat\Phi_A C)^{-1} C^\top (\hat\beta - \beta) $$ and F* = λF follows exact Fc, m distribution.
When we have only one coefficient to test, then C is a column vector containing a single 1 inside. We can still follow the same calculations as for the multi-dimensional case. This recovers the degrees of freedom results of the Satterthwaite method.
λ and m can be calculated through
M = C(C⊤ΦC)−1C⊤
$$ A_1 = \sum_{h=1}^k{\sum_{j=1}^k{W_{hj} tr(M \Phi P_h \Phi) tr(M \Phi P_j \Phi)}} $$
$$ A_2 = \sum_{h=1}^k{\sum_{j=1}^k{W_{hj} tr(M \Phi P_h \Phi M \Phi P_j \Phi)}} $$
$$ 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 B}{(1-c_2 B)^2(1-c_3 B)}\right\}}$$
$$\rho = \frac{V^{*}}{2(E^*)^2}$$
$$m = 4 + \frac{c+2}{c\rho - 1}$$ $$\lambda = \frac{m}{E^*(m-2)}$$
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 elements of the Cholesky factor of Σ (see parameterization), the second-order derivatives of Σ over our parameters, are non-zero matrices. However, if we use the elements of Σ as our parameters, then the second-order derivatives are zero matrices. However, the differences can be small and will not affect the inference. If you would like to match SAS results for the unstructured covariance model, you can use the linear Kenward-Roger approximation.
mmrmIn package mmrm, we have implemented Kenward-Roger
calculations based on the previous sections. Specially, for the
first-order and second-order derivatives, we use automatic
differentiation to obtain the results easily for non-spatial covariance
structure. For spatial covariance structure, we derive the exact
results.
For spatial exponential covariance structure, we have
θ = (θ1, θ2) σ = eθ1 $$\rho = \frac{e^{\theta_2}}{1 + e^{\theta_2}}$$
Σij = σρdij where dij 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) $$