Here we describe the details of the Satterthwaite degrees of freedom calculations.
In Christensen (2018) 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.
We are also using the same notation as in the Details of the Kenward-Roger calculations. In particular, we assume we have a contrast matrix C ∈ ℝc × p with which we want to test the linear hypothesis Cβ = 0. Further, 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 θ̂. Φ(θ) = {X⊤Ω(θ)−1X}−1 is the asymptotic covariance matrix of β̂.
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(θ̂) = CΦ(θ̂)C⊤ is the scalar in the numerator and we can identify it as the variance estimate for the estimated scalar contrast Cβ̂. The computational challenge is essentially to evaluate the denominator in the expression for ν̂(θ̂), which amounts to computing the k-dimensional gradient f′(θ̂) of f(θ) (for the given contrast matrix C) at the estimate θ̂. We already have the variance-covariance matrix W(θ̂) of the variance parameter vector θ from the model fitting.
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 β. It is better to instead evaluate
the gradient of the matrix valued function Φ(θ), which is therefore
the Jacobian, with regards to θ, 𝒥(θ) = ∇θΦ(θ).
Imagine 𝒥(θ) as the the
3-dimensional array with k
faces of size p × p.
Left and right multiplying each face by C and C⊤ respectively leads to
the k-dimensional gradient
f′(θ) = C𝒥(θ)C⊤.
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
𝒥(θ̂), it is only a matter of
putting the different quantities together to compute the estimate of the
denominator degrees of freedom, ν̂(θ̂).
Currently, we evaluate the gradient of Φ(θ) 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 coefficient (Φ(θ)−1), following the Kenward-Roger calculations. Please note that we only need Ph matrices.
Then, to obtain the Jacobian of the covariance matrix of coefficient, following the algorithm, we use Φ(θ) 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 θ) of matrices of p × p, where p is the dimension of β.
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^\top (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 Christensen (2018), 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Φ(θ̂)C⊤.
The F-statistic can be rewritten as a sum of t2 statistics based on
these eigen-values. The corresponding random variables are independent
(by design because they are derived from the orthogonal eigen-vectors)
and essentially have one degree of freedom each. Hence, each of the
t statistics is treated as
above in the one-dimensional contrast case, i.e. the denominator degree
of freedom is calculated for each of them. Finally, using properties of
the F distribution’s expectation, the denominator degree of freedom for
the whole F statistic is derived.
In Bell and McCaffrey (2002) the Satterthwaite degrees of freedom in combination with a sandwich covariance matrix estimator are described.
For one-dimensional contrast, following the same notation in Details of the model fitting in
mmrm and Details of the
Kenward-Roger calculations, we have the following derivation. For an
estimator of variance with the following term
v = sc⊤(X⊤X)−1∑iXi⊤Aiϵiϵi⊤AiXi(X⊤X)−1c
where s takes the value of $\frac{n}{n-1}$, 1 or $\frac{n-1}{n}$, and Ai takes Ii, $(I_i - H_{ii})^{-\frac{1}{2}}$, or (Ii − Hii)−1 respectively, c is a column vector, then v can be decomposed into the a weighted sum of independent χ12 distribution, where the weights are the eigenvalues of the n × n matrix G with elements Gij = gi⊤Vgj
where
$$ g_i = s^{\frac{1}{2}} (I - H)_i^\top A_i X_i (X^\top X)^{-1} c $$ H = X(X⊤X)−1X⊤
(I − H)i corresponds to the rows of subject i.
So the degrees of freedom can be represented as $$ \nu = \frac{(\sum_{i}\lambda_i)^2}{\sum_{i}{\lambda_i^2}} $$
where λi, i = 1, …, n are the eigenvalues of G. Bell and McCaffrey (2002) also suggests that V can be chosen as identify matrix, so Gij = gi⊤gj.
Following Weighted Least Square Estimator, we can transform the original X into x̃ to use the above equations.
To avoid repeated computation of matrix Ai, H etc for different contrasts, we calculate and cache the following
Gi* = (I − H)i⊤AiXi(X⊤X)−1 which is a ∑imi × p matrix. With different contrasts, we need only calculate the following gi = Gi*c to obtain a ∑imi × 1 matrix, G can be computed with gi.
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}\lambda_i)^2}{\sum_{i}{\lambda_i^2}} = \frac{tr(G)^2}{\sum_{i}{\sum_{j}{G_{ij}^2}}} $$
The scale parameter is not used throughout the package.
The proof is as following
Let A has dimension p × q, B has dimension q × p $$ tr(AB) = \sum_{i=1}^{p}{(AB)_{ii}} = \sum_{i=1}^{p}{\sum_{j=1}^{q}{A_{ij}B_{ji}}} $$
$$ tr(BA) = \sum_{i=1}^{q}{(BA)_{ii}} = \sum_{i=1}^{q}{\sum_{j=1}^{p}{B_{ij}A_{ji}}} $$
so tr(AB) = tr(BA)
Following eigen decomposition, we have G = QΛQ⊤ where Λ is diagonal matrix, Q is orthogonal matrix.
Using the previous formula that tr(AB) = tr(BA), we have
tr(G) = tr(QΛQ⊤) = tr(ΛQ⊤Q) = tr(Λ) = ∑i(λi)
tr(G⊤G) = tr(QΛQ⊤QΛQ⊤) = tr(Λ2Q⊤Q) = tr(Λ2) = ∑i(λi2)
and tr(G⊤G) can be further expressed as
tr(G⊤G) = ∑i(G⊤G)ii = ∑i∑jGij⊤Gji = ∑i∑jGij2
For multi-dimensional contrast we use the same technique for multi-dimensional contrast for asymptotic covariance.