--- title: "Kenward-Roger" package: mmrm bibliography: '`r system.file("REFERENCES.bib", package = "mmrm")`' csl: '`r system.file("jss.csl", package = "mmrm")`' output: rmarkdown::html_vignette: toc: true vignette: | %\VignetteIndexEntry{Kenward-Roger} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` 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`](algorithm.html). 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](algorithm.html#covariance-matrix-model), 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, @kackar1984 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$, @kenward1997 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](https://en.wikipedia.org/wiki/Matrix_calculus#Identities_in_differential_form). \[ \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](#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](#scenario-under-weighted-mmrm), 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](#special-considerations-for-mmrm-models) 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](#scenario-under-group-specific-covariance-estimates), $\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). ## 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](https://en.wikipedia.org/wiki/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})$, where $\cdot \odot \cdot$ is the [Hadamard product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices)). 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](satterthwaite.html#one-dimensional-contrast) 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](algorithm.html#covariance-matrix-model)), 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. ## Implementation 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](#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. The closest results to SAS are again obtain when using the linear Kenward-Roger method in `mmrm`. # References