1.1 Motivation

Linear mixed models are widely used in biomedical research for inference in analyses with missing data. Kenward and Roger [1] described a scaled Wald statistic and null case reference distribution for tests of fixed effects in the linear mixed model. Despite the widespread use of the Kenward and Roger [1] method for data analysis, no general methods are available to calculate power for the Kenward and Roger [1] test.

Several authors have described power approximations for related tests and models. Helms [2] described a noncentral power approximation for a Wald test. Helms used a different null case reference distribution than the one derived by Kenward and Roger. Stroup [3] suggested an “exemplary data” approach for calculating power for mixed models with missing data. Tu et al. [4, 5] developed an asymptotic power approximation based on generalized estimating equations. Shieh [6] provided non-central power approximations for multivariate models with random covariates and no missing data. Chi, Glueck, and Muller [7] demonstrated that power methods for the general linear multivariate model may be used in complete, balanced, homoscedastic mixed models.

We derive a noncentral power approximation for the Kenward and Roger [1] test for a broad range of models. We use a method of moments approach [8] to form an approximate distribution of the Kenward and Roger [1] scaled Wald statistic, F_R, under the alternative. The reference distribution of F_R under the alternative depends on the approximate moments of the unscaled Wald statistic.

The remainder of the manuscript is organized as follows. In Section 2, we introduce notation for the general linear mixed model and briefly review the methods of Kenward and Roger [1]. In Section 3, we describe a noncentral power approximation for the Kenward and Roger [1] test. In Section 4, we summarize the Monte Carlo simulation study used to evaluate the power approximation. In Section 5, we demonstrate a power calculation for a longitudinal trial in oral cancer prevention. In Section 6, we provide concluding remarks.

2.1 Notation

For i ∈ {1, …, n}, let a = {a_i} denote an n × 1 column vector. Furthermore, for i ∈ {1, …, n} and j ∈ {1, …, m}, let A = {a_ij} indicate an n × m matrix with transpose A′ = {a_ji}. Let I_d be a (d × d) identity matrix. For a matrix A = [a₁ a₂ … a_n], let . Define the Kronecker product of two matrices A and B as A ⊗ B = {a_ij B} [9, Section 1.3].

Extend the direct sum operator [9, Section 1.3] to sets of arbitrarily sized matrices as follows. Let {A₁, …, A_J} be a set of matrices such that A_j has dimension (r_j × c_j). Let be an (r_i × c_j) matrix of zeros. Define the direct sum of {A₁, …, A_J} as (1)

For δ ∈ {1, …, (2^p − 1)} and d ∈ {1, …, δ}, define the set R_d where R_d ⊆ {1, …, p} of cardinality 1 ≤ p_d ≤ p. For every R_d, let D_p,d, a deletion matrix, be the (p_d × p) submatrix of I_p formed by keeping each row i of I_p such that i ∈ R_d. For example, given a (p × p) matrix A and R_d = {1, 3}, (2) and (3)

Let E₀(u) and E_A(u) be the expectations of the random variable u under the null and alternative hypotheses, respectively. Similarly, let and indicate the variance under the null and alternative hypotheses. For random matrix variates, denote the covariance under the null and alternative hypotheses as and , respectively.

Let indicate that random variable X follows a distribution D exactly, while indicates that distribution is followed approximately. Let indicate that the random variable F follows a noncentral distribution [10] with numerator degrees of freedom ν_n, denominator degrees of freedom ν_d, and noncentrality parameter γ. For γ = 0, F is said to follow a central distribution, written . Define such that for 0 ≤ b ≤ 1 (4)

Use to indicate that the (N × p) matrix Y follows a matrix Gaussian distribution, with M an (N × p) matrix of means, Ξ an (N × N) symmetric, positive definite column covariance matrix, and Σ a (p × p) symmetric, positive definite row covariance matrix [9, Chapter 8]. Write to indicate that the (p × p) matrix W follows a central Wishart distribution of dimension p, degrees of freedom N, on covariance Σ. For Ψ = Σ⁻¹, write to indicate that W⁻¹ follows a central inverse Wishart distribution of dimension p, degrees of freedom N+ p+ 1, and precision matrix Ψ [11, p. 111, Theorem 3.4.1].

2.2 The general linear mixed model

We describe the general linear mixed model for Gaussian outcomes using the notation of Muller and Stewart [9, Chapter 5]. Let i ∈ {1, …, N} indicate the ith independent sampling unit [9, Chapter 5]. An independent sampling unit may be a single participant, as in a clinical trial, or a group of participants, as in a cluster-randomized study. Observations from two different independent sampling units are statistically independent. Observations within an independent sampling unit may be correlated. For example, for a particpant in a longitudinal trial, repeated measurements over time will be correlated.

Let p_i be the number of observations for the ith independent sampling unit, with p = max_i(p_i). For the ith independent sampling unit, let y_i be the (p_i × 1) vector of observed outcomes, X_i be the (p_i × r) fixed effects design matrix of rank r, and e_i be the (p_i × 1) vector of random errors. Assume that for i ≠ j, e_i ⊥ e_j and y_i ⊥ y_j. Let Σ_i be a (p_i × p_i) symmetric, positive definite matrix, with (5) Let β be the (r × 1) vector of regression parameters. The linear mixed model for the ith independent sampling unit is (6)

Let . Define the (n × 1) vectors and . Stack the fixed effect design matrices into the (n × r) matrix (7) Throughout, we assume that predictor values are not allowed to change within an independent sampling unit, i.e., that there are no repeated covariates. In addition, we assume that all predictor values are fixed as part of the study design. The population-averaged form of the linear mixed model is (8) Define (9) The distribution of y_s is (10)

2.3 Tests for fixed effects in mixed models

Let α be the Type I error rate. Let C be the (a × r) matrix of fixed effects contrasts. Define the (a × 1) matrix θ = Cβ, and let θ₀ be the (a × 1) matrix of null values. The general linear hypothesis may be stated as (11)

In order to conduct power analysis for the general linear hypothesis in the mixed model, we must consider the target estimation method. Several estimation methods have been described for mixed models [12, Chapter 5]. Common estimation methods include restricted maximum likelihood and maximum likelihood.

Let m indicate the estimation method. Let and be the estimates of Σ_s and β obtained from method m. Define . The Wald statistic for the linear mixed model is (12)

The distribution of the Wald statistic is not known exactly for any m. Various reference distributions have been suggested for each estimation method m. In general, the distributions share a common form, with (13) Under the null hypothesis, γ_m = 0 and .

2.4 The Kenward-Roger test for fixed effects

Kenward and Roger [1] suggested using restricted maximum likelihood estimation (m = R) and a scaled Wald statistic. (14) Kenward and Roger [1] used Taylor expansion to estimate E₀(w_R) and from observed data. Kenward and Roger [1] substituted E₀(w_R) and into method of moments approximations for λ and the reference distribution of F_R under the null. With , (15) (16) and (17)

3.1 The approximate moments of the Wald statistic

We derive a noncentral power approximation for the Kenward and Roger [1] test. The method of moments approach [8] is used to form an approximate distribution of the Kenward and Roger [1] scaled Wald statistic, F_R, under the alternative. The reference distribution of F_R under the alternative depends on the approximate moments of the unscaled Wald statistic.

We demonstrate that the Wald statistic has an approximately noncentral reference distribution under the alternative and a central reference distribution under the null. The result depends on approximate distributional results for both and . Because distributional results are, in general, not available for restricted maximum likelihood estimation, we instead use distributional results based on other techniques.

Let m = W indicate weighted least squares, and m = M denote multivariate methods. Approximate by , which is Gaussian, conditional on Σ_s. The term can be approximated by . We show that is approximately Wishart. Finally, under the assumption of independence, we combine the terms to obtain an approximate distribution.

3.2 A three-moment approximation for the distribution of the Kenward and Roger scaled Wald statistic under the alternative hypothesis

We use a method of moments approach [8] to form the approximate distribution of Kenward and Roger [1] scaled Wald statistic, F_R, under the alternative. The parameters of the distribution depend on the approximate Wald moments derived in Section 3.1. We approximate the distribution of the Kenward and Roger [1] statistic, F_R = λw_R, by the distribution of F = λw, where . Thus (33) To obtain values for λ, ν, and γ under the alternative, we match three moments, setting (34) (35) and (36) With (37) we obtain (38) (39) and (40) When γ = 0, Eq 39 reduces to (41) which shares the same form as the result obtained by Kenward and Roger (Eq 16). The exact values of ρ, and hence ν, will differ due to the disparate techniques used to obtain moments for the Wald statistics, w and w_R.

3.3 Power calculation for the Kenward and Roger test

We calculate power for the Kenward and Roger test as follows. Define α, Σ_max, β, C and θ₀. For i ∈ {1, …, N}, specify X_i and R_d. Calculate a, ν, and γ as described in Section 3.2. Form the reference distribution of . Using the approximate reference distribution of F_R under the null, , find the critical value (42) Finally, using the approximate reference distribution of F_R under the alternative, , calculate power as (43)

4.1 Methods

We compared approximate power values, calculated as in Section 3.3, with empirical power for two types of study designs: unbalanced, cluster randomized trials and longitudinal studies with known dropout patterns. Approximate power was calculated using our mixedPower package for R version 4.0.2 [15].

Empirical power was calculated by Monte Carlo simulation in SAS [16, version 9.4]. We defined α, Σ_max, β, C and θ₀. For i ∈ {1, …, N}, we specified X_i and R_d. We generated 10, 000 replicates of e_s and computed y_s as in Eq 8. For each replicate, we tested the linear contrast C using SAS PROC MIXED with the DDFM = KenwardRoger flag to request Kenward and Roger [1] denominator degrees of freedom. Empirical power was estimated as the proportion of replicates for which the null hypothesis was rejected. Source code is available at http://github.com/SampleSizeShop/mixedPower.

4.2 Results

Fig 1 summarizes the deviations between the approximate and the empirical power values. The three box plots show results for all designs, for cluster randomized trials, and for longitudinal studies. Overall, the median deviation between the approximate and the empirical power values was 0.010 (min: −0.010, 1st quartile: 0.005, 3rd quartile: 0.015, max: 0.064). For cluster randomized trials, the median deviation was 0.011, (min: −0.001, 1st quartile: 0.006, 3rd quartile: 0.017, max: 0.064). For longitudinal designs, the median deviation was 0.003, (min: −0.010, 1st quartile: 0.000, 3rd quartile: 0.009, max: 0.016).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Power deviations for all designs, cluster randomized designs only, and longitudinal designs only.

(center line, median; box limits, 1st and 3rd quartiles; whiskers, minimum and maximum).

https://doi.org/10.1371/journal.pone.0254811.g001

Further details for cluster-randomized designs are shown in Fig 2. The accuracy of the power approximation improved with larger cluster sizes. The approximation retained accuracy regardless of the ratio of incomplete to complete cluster sizes. As shown in Table 1, accuracy was similar across ICC values, with slight improvements with increasing correlation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Power deviations for cluster randomized designs.

(center line, median; box limits, 1st and 3rd quartiles; whiskers, minimum and maximum).

https://doi.org/10.1371/journal.pone.0254811.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Deviations between approximate and empirical power in cluster randomized designs by ICC.

https://doi.org/10.1371/journal.pone.0254811.t001

Results for longitudinal designs are shown in Fig 3. The power approximation was highly accurate for all longitudinal designs tested.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Power deviations in longitudinal designs.

(center line, median; box limits, 1st and 3rd quartiles; whiskers, minimum and maximum).

https://doi.org/10.1371/journal.pone.0254811.g003

A Appendix: Theorems and proofs

Theorem 1. For m ∈ {1, …, k}, let p_m ∈ {1, 2, …, p}, N_m > (p_m + 3) and define Ψ_m = {ψ_mij} to be a (p_m × p_m) symmetric, positive definite matrix. Define a set of k ≥ 2, independent, non-identically distributed inverse central Wishart random matrices, such that for m ∈ {1, …, k}, . For i ∈ {1, …, q} and R_m ⊂ {1, 2, …, p} of cardinality p_m, define X_m to be a (p_m × qp) matrix of rank p_m < qp with the form (53) If for each i ∈ {1, …, q}, there exists at least one m such that X_m = I_q({i}) ⊗ I_p, then (54) is positive definite.

Proof. Let Q_i = {X_m: X_m = I_q({i}) ⊗ I_p(R_m)}. Then (55) Note that for i ∈ {1, 2, ‥, q}, I_q({i})′ I_q({i}) is a (q × q) matrix for which the ith diagonal element is 1 and all remaining elements are 0. Therefore, Eq 55 can be equivalently expressed as a direct sum. (56)

From Mathai and Provost [18, p.18, Theorem 2.2b.1], it follows that each is positive semi-definite. By assumption, for each Q_i, there exists a c_i such that such that . Then (57) Because is positive definite and the remaining are positive semi-definite for i ∈ {1, …, q}, then (58) is positive definite.

Since is a block matrix, the eigenvalues of are the eigenvalues of all of the blocks. Since each block (Eq 58) is positive definite and hence has positive eigenvalues, it follows that must also be positive definite.

Theorem 2. For m ∈ {1, …, k}, i ∈ {1, …, q}, R_m ⊂ {1, 2, …, p} of cardinality p_m, X_m = I_q({i}) ⊗ I_p(R_m) a (p_m × qp) matrix of rank p_m < qp, N_m > (p_m + 3), Ψ_m a (p_m × p_m) symmetric, positive definite matrix, and , (59) Proof. Let Dg(x) indicate a square matrix with the elements of the vector x on the diagonal.

Since is positive definite and has full rank, then by Lemma 1.24 (a) of Muller and Stewart [9], it has the spectral decomposition (60) where λ is the (p_m × 1) vector of eigenvalues and V is the (p_m × p_m) orthogonal matrix of eigenvectors of . Then (61)

Since X_m has deficient rank p_m < qp, then by Lemma 1.25 of of Muller and Stewart [9] it must have qp − p_m zero eigenvalues. Let λ₀ be the (qp − p_m × 1) vector of zero eigenvalues and V₀ the [qp × (qp − p_m)] matrix of corresponding eigenvectors. Then (62)

Selecting V₀ such that , and X_m V₀ = 0, ensures that is orthogonal. Then Eq 62 is the spectral decomposition of , with eigenvalues .

Since λ₀ contains only zero eigenvalues and using the definition of the trace, (63) Theorem 3. For m ∈ {1, …, k}, let p_m ∈ {1, …, p}, N_m > (p_m+ 3) and let Ψ_m = {ψ_mij} be a (p_m × p_m) symmetric, positive definite matrix. Define a set of k ≥ 2, independent, non-identically distributed inverse central Wishart random matrices, such that for m ∈ {1, …, k}, . For i ∈ {1, …, q} and R_m ⊂ {1, …, p} of cardinality p_m, define X_m to be a (p_m × qp) matrix of rank p_m < qp with the form (64) Under the assumption that for each i ∈ {1, …, q}, there exists at least one m such that X_m = I_q({i}) ⊗ I_p, it can be shown that (65) is approximately distributed as .

Proof. Theorem 1 in Appendix demonstrates that Q⁻¹ is positive definite under the restriction that for each i ∈ {1, …, q}, there exists at least one m such that X_m = I_q({i}) ⊗ I_p.

To derive an approximate distribution for Q⁻¹, we match the expectation of the sum of the Wishart matrices and the variance of the trace of the sum of the Wishart matrices. Set (66) and (67)

From Theorem 2 in Appendix and the independence of the , (68)

Then the approximate parameters for are (69) and (70) where (71) (72) (73) (74) (75) and (76) The method of moments approximation yields an asymptotic approximation for the sum, as desired.

Theorem 4. Let n and p be positive integers, μ be a (p × 1) vector of means, and Σ_x ≠ Σ_W be symmetric and positive definite (p × p) matrices. Suppose independently of . Then (77) with (78) (79) and (80) Proof. Define . Define . Using Lemma 17.10 in Arnold [19, p. 319], it follows that . Hence, V ⊥ x, which implies V ⊥ U.

The expression U is a weighted sum of noncentral χ² random variables [9, Theorem 9.5, p. 176]. Approximate the distribution of U with a single noncentral χ², so that . Using the approach described by Kim et al. [8], obtain values for λ_u, n_u and δ_u by matching the following three moments: (81) (82) and (83) The moments of U are [9, Corollary 9.6.3, p. 179], (84) (85) and (86) Then the approximate parameters of U are (87) (88) and (89) Since , , and V ⊥ U, Because U/V = x′W⁻¹ x, and the result follows.