--- title: "Satterthwaite" 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{Satterthwaite} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console markdown: wrap: 72 --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` Here we describe the details of the Satterthwaite degrees of freedom calculations. ## Satterthwaite degrees of freedom for asymptotic covariance In @Christensen2018 the Satterthwaite degrees of freedom approximation based on normal models is well detailed and the computational approach for models fitted with the `lme4` package is explained. We follow the algorithm and explain the implementation in this `mmrm` package. The model definition is the same as in [Details of the model fitting in `mmrm`](algorithm.html). We are also using the same notation as in the [Details of the Kenward-Roger calculations](kenward.html). In particular, we assume we have a full row rank contrast matrix $C \in \mathbb{R}^{c\times p}$ with which we want to test the linear hypothesis $C\beta = 0$. Further, $W(\hat\theta)$ is the covariance estimate of $\hat\theta$: the inverse observed information, obtained by inverting the Hessian of the negative log-likelihood evaluated at $\hat\theta$ (the restricted likelihood for REML fits). $\Phi(\theta) = \left\{X^\top \Omega(\theta)^{-1} X\right\} ^{-1}$ is the asymptotic covariance matrix of $\hat\beta$, and $k$ is the number of covariance parameters in $\theta$. ### One-dimensional contrast We start with the case of a one-dimensional contrast, i.e. $c = 1$. The Satterthwaite adjusted degrees of freedom for the corresponding t-test are then defined as: \[ \hat\nu(\hat\theta) = \frac{2f(\hat\theta)^2}{f{'}(\hat\theta)^\top W(\hat\theta) f{'}(\hat\theta)} \] where $f(\hat\theta) = C \Phi(\hat\theta) C^\top$ is the scalar in the numerator and we can identify it as the variance estimate for the estimated scalar contrast $C\hat\beta$. The computational challenge is essentially to evaluate the denominator in the expression for $\hat\nu(\hat\theta)$, which amounts to computing the $k$-dimensional gradient $f{'}(\hat\theta)$ of $f(\theta)$ (for the given contrast matrix $C$) at the estimate $\hat\theta$. We already have the variance-covariance matrix $W(\hat\theta)$ of the variance parameter vector $\theta$ from the model fitting. #### Jacobian approach However, if we proceeded in a naive way here, we would need to recompute the denominator again for every chosen $C$. This would be slow, e.g. when changing $C$ every time we want to test a single coefficient within $\beta$. It is better to instead evaluate the gradient of the matrix valued function $\Phi(\theta)$, which is therefore the Jacobian, with regards to $\theta$, $\mathcal{J}(\theta) = \nabla_\theta \Phi(\theta)$. Imagine $\mathcal{J}(\theta)$ as the 3-dimensional array with $k$ faces of size $p\times p$. Left and right multiplying each face by $C$ and $C^\top$ respectively leads to the $k$-dimensional gradient $f'(\theta) = C \mathcal{J}(\theta) C^\top$. Therefore for each new contrast $C$ we just need to perform simple matrix multiplications, which is fast (see `h_gradient()` where this is implemented). Thus, having computed the estimated Jacobian $\mathcal{J}(\hat\theta)$, it is only a matter of putting the different quantities together to compute the estimate of the denominator degrees of freedom, $\hat\nu(\hat\theta)$. #### Jacobian calculation Currently, we evaluate the gradient of $\Phi(\theta)$ through function `h_jac_list()`. It uses automatic differentiation provided in `TMB`. We first obtain the Jacobian of the inverse of the covariance matrix of coefficients ($\Phi(\theta)^{-1}$), following the [Kenward-Roger calculations](kenward.html#special-considerations-for-mmrm-models). Please note that we only need $P_h$ matrices. Then, to obtain the Jacobian of the covariance matrix of coefficients, following the [algorithm](kenward.html#derivative-of-the-sigma-1), we use $\Phi(\theta)$ estimated in the fit to obtain the Jacobian. The result is a list (of length $k$ where $k$ is the dimension of the variance parameter $\theta$) of matrices of $p \times p$, where $p$ is the dimension of $\beta$. Because only the $P_h$ matrices are needed, the Satterthwaite implementation initializes a first-derivative-only cache: for non-spatial covariance structures, it differentiates the Cholesky factor once and caches the covariance and inverse covariance first derivatives for each observed-visit pattern. It does not run nested automatic differentiation or allocate second-derivative caches. Spatial covariance first derivatives are evaluated analytically on demand. Neither $Q_{hj}$ nor $R_{hj}$ is constructed for the Satterthwaite Jacobian; the gradient and degrees-of-freedom formulas remain unchanged. #### Connection to the scalar Kenward-Roger shortcut For a REML fit, the [Kenward-Roger implementation](kenward.html#one-dimensional-shortcut) can obtain the same scalar degrees of freedom directly from its cached $P_h$ matrices. Write $C = l^\top$, with $l \in \mathbb{R}^p$, and evaluate all quantities at $\hat\theta$, omitting this argument below. The Jacobian identity \[ \frac{\partial\Phi}{\partial\theta_h} = -\Phi P_h\Phi \] implies $f'_h = -l^\top\Phi P_h\Phi l$. Define $a_h = -f'_h/f$. The Satterthwaite formula therefore becomes \[ \hat\nu = \frac{2}{a^\top Wa}. \] For one-dimensional KR, $A_1 = A_2 = a^\top Wa$ and the F scale is exactly $\lambda = 1$. `h_kr_df()` uses this scalar shortcut for both full and linear KR, while Satterthwaite continues to use its cached Jacobian through `h_gradient()`. The equality uses the **unadjusted** covariance $\Phi$ and the same $W$; KR standard errors still use the adjusted covariance $\Phi_A$, so equal scalar degrees of freedom do not imply equal test statistics, p-values, or confidence intervals. For multiple contrasts, KR uses a normalized contrast-space contraction of its moment quantities; the Satterthwaite eigen-decomposition described below remains unchanged. ### Multi-dimensional contrast When $c > 1$ we are testing multiple contrasts at once. Here an F-statistic \[ F = \frac{1}{c} (C\hat\beta)^\top (C \Phi(\hat\theta) C^\top)^{-1} (C\hat\beta) \] is calculated, and we are interested in estimating an appropriate denominator degrees of freedom for $F$, while assuming $c$ are the numerator degrees of freedom. Note that only in special cases, such as orthogonal or balanced designs, the F distribution will be exact under the null hypothesis. In general, it is an approximation. The calculations are described in detail in @Christensen2018, and we don't repeat them here in detail. The implementation is in `h_df_md_sat()` and starts with an eigen-decomposition of the asymptotic variance-covariance matrix of the contrast estimate, i.e. $C \Phi(\hat\theta) C^\top$. This rewrites $cF$ as a sum of squared standardized contrasts. Under the null hypothesis, each standardized contrast is approximated by a $t_{\nu_a}$ distribution, and its square by an $F_{1,\nu_a}$ distribution, where $\nu_a$ is calculated using the one-dimensional Satterthwaite formula, for $a = 1, \dotsc, c$. The eigen-decomposition diagonalizes the estimated contrast covariance; it does not guarantee independence of the studentized statistics, whose denominators are estimated from the data. When the component degrees of freedom exceed two, matching the approximate expectation of the sum to that of $cF_{c,\nu}$ gives the overall denominator degrees of freedom. This expectation calculation uses linearity of expectation and does not require independence. Numerically equal component degrees of freedom are returned directly. For unequal components, the implementation returns two denominator degrees of freedom if any component has at most two. ## Satterthwaite degrees of freedom for empirical covariance In @bell2002bias the Satterthwaite degrees of freedom in combination with a sandwich covariance matrix estimator are described. ### One-dimensional contrast For one-dimensional contrast, following the same notation in [Details of the model fitting in `mmrm`](algorithm.html) and [Details of the Kenward-Roger calculations](kenward.html), we have the following derivation. Let $C = l^\top$ be the contrast, with a column vector $l \in \mathbb{R}^p$. First consider ordinary least squares. Distinguish the model errors $\epsilon = Y - X\beta$ from the fitted residuals $e = Y - X\hat\beta = (I-H)\epsilon$, where \[ H = X(X^\top X)^{-1}X^\top. \] Write $e_i$ for the residuals of subject $i$. The sandwich estimator of the variance of $l^\top\hat\beta$ is \[ v = s l^\top(X^\top X)^{-1}\sum_{i}{X_i^\top A_i e_i e_i^\top A_i X_i} (X^\top X)^{-1} l \] where $s$ takes the value of $\frac{n}{n-1}$, $1$ or $\frac{n-1}{n}$, and $A_i$ takes $I_i$, $(I_i - H_{ii})^{-\frac{1}{2}}$, or $(I_i - H_{ii})^{-1}$ respectively (as in the [empirical covariance with weighted least squares](empirical_wls.html)). Here $I_i$ is the $m_i\times m_i$ identity and $H_{ii}$ is the subject's diagonal block of $H$. Under a working normal model for $\epsilon$, with the design and adjustment matrices treated as fixed, $v$ has the distribution of a weighted sum of independent $\chi_1^2$ variables. The weights are the eigenvalues of the $n\times n$ matrix $\Gamma$ with elements \[ \Gamma_{ij} = \gamma_i^\top V \gamma_j \] where \[ \gamma_i = s^{\frac{1}{2}} (I - H)_i^\top A_i X_i (X^\top X)^{-1} l \] $(I - H)_i$ corresponds to the rows of subject $i$, so that $\gamma_i \in \mathbb{R}^N$ with $N = \sum_i m_i$ observations in total. $V = \operatorname{Var}(\epsilon) \in \mathbb{R}^{N \times N}$ is the working covariance matrix of the model errors. The fitted residuals have covariance $(I-H)V(I-H)^\top$; their projection is already included in $\gamma_i$. In particular, $v = \sum_i(\gamma_i^\top\epsilon)^2$. So the degrees of freedom can be represented as \[ \nu = \frac{(\sum_{i}\omega_i)^2}{\sum_{i}{\omega_i^2}} \] where $\omega_i, i = 1, \dotsc, n$ are the eigenvalues of $\Gamma$. @bell2002bias also suggests that $V$ can be chosen as identity matrix, so $\Gamma_{ij} = \gamma_i^\top \gamma_j$. For generalized least squares, apply these equations to the transformed response and design from the [weighted least squares estimator](algorithm.html#weighted-least-squares-estimator). Specifically, factor $\hat\Omega = L_\Omega L_\Omega^\top$ and use $Y^\dagger = L_\Omega^{-1}Y$, $X^\dagger = L_\Omega^{-1}X$, and $\epsilon^\dagger = L_\Omega^{-1}\epsilon$. The hat matrix and fitted residuals are then computed in these transformed coordinates. With the fitted covariance as the working model, the transformed errors have working covariance $V = I$; the transformation is treated as fixed for this approximation. Below, $X$ and $H$ refer to these transformed coordinates when applying the formulas to `mmrm`. To avoid repeated computation of matrix $A_i$, $H$ etc for different contrasts, we calculate and cache the following \[ \Gamma^\ast_i = (I - H)_i^\top A_i X_i (X^\top X)^{-1} \] which is an $N \times p$ matrix. With different contrasts, we need only calculate the following \[ \gamma_i = s^{\frac{1}{2}} \Gamma^\ast_i l \] to obtain an $N \times 1$ matrix, and $\Gamma$ can be computed with the $\gamma_i$. To obtain the degrees of freedom, and to avoid eigen computation on a large matrix, we can use the following equation \[ \nu = \frac{(\sum_{i}\omega_i)^2}{\sum_{i}{\omega_i^2}} = \frac{\operatorname{tr}(\Gamma)^2}{\sum_{i}{\sum_{j}{\Gamma_{ij}^2}}} \] The common factor $s$ cancels from this degrees-of-freedom ratio, so it need not be included in the calculation. The trace identities used here are proved in the [appendix](#appendix-trace-identities). ### Multi-dimensional contrast For multiple contrasts, we apply the same eigen-decomposition and expectation-matching technique as for asymptotic covariance, using the empirical covariance matrix and the empirical scalar degrees of freedom for each component. ## Appendix: Trace identities ### Cyclic invariance of trace We first show \[ \operatorname{tr}(AB) = \operatorname{tr}(BA) \] Let $A$ have dimension $r\times q$, $B$ have dimension $q\times r$ \[ \operatorname{tr}(AB) = \sum_{i=1}^{r}{(AB)_{ii}} = \sum_{i=1}^{r}{\sum_{j=1}^{q}{A_{ij}B_{ji}}} \] \[ \operatorname{tr}(BA) = \sum_{i=1}^{q}{(BA)_{ii}} = \sum_{i=1}^{q}{\sum_{j=1}^{r}{B_{ij}A_{ji}}} \] so $\operatorname{tr}(AB) = \operatorname{tr}(BA)$ ### Trace and squared eigenvalues We next show \[ \operatorname{tr}(\Gamma) = \sum_{i}(\omega_i) \] and \[ \sum_{i}(\omega_i^2) = \sum_{i}{\sum_{j}{\Gamma_{ij}^2}} \] if $\Gamma = \Gamma^\top$ Following eigen decomposition, we have \[ \Gamma = U \operatorname{diag}(\omega) U^\top \] where $\operatorname{diag}(\omega)$ is the diagonal matrix of the eigenvalues, and $U$ is an orthogonal matrix. Using the previous formula that $\operatorname{tr}(AB) = \operatorname{tr}(BA)$, we have \[ \operatorname{tr}(\Gamma) = \operatorname{tr}(U \operatorname{diag}(\omega) U^\top) = \operatorname{tr}(\operatorname{diag}(\omega) U^\top U) = \operatorname{tr}(\operatorname{diag}(\omega)) = \sum_{i}(\omega_i) \] \[ \operatorname{tr}(\Gamma^\top \Gamma) = \operatorname{tr}(U \operatorname{diag}(\omega) U^\top U \operatorname{diag}(\omega) U^\top) = \operatorname{tr}(\operatorname{diag}(\omega)^2 U^\top U) = \operatorname{tr}(\operatorname{diag}(\omega)^2) = \sum_{i}(\omega_i^2) \] and $\operatorname{tr}(\Gamma^\top \Gamma)$ can be further expressed as \[ \operatorname{tr}(\Gamma^\top \Gamma) = \sum_{i}{(\Gamma^\top \Gamma)_{ii}} = \sum_{i}{\sum_{j}{\Gamma^\top_{ij}\Gamma_{ji}}} = \sum_{i}{\sum_{j}{\Gamma_{ij}^2}} \] # References