A power approximation for the Kenward and Roger Wald test in the linear mixed model

Abstract

We derive a noncentral power approximation for the Kenward and Roger test. We use a method of moments approach to form an approximate distribution for the Kenward and Roger scaled Wald statistic, under the alternative. The result depends on the approximate moments of the unscaled Wald statistic. Via Monte Carlo simulation, we demonstrate that the new power approximation is accurate for cluster randomized trials and longitudinal study designs. The method retains accuracy for small sample sizes, even in the presence of missing data. We illustrate the method with a power calculation for an unbalanced group-randomized trial in oral cancer prevention.

About the Authors:

Sarah M. Kreidler

Contributed equally to this work with: Sarah M. Kreidler, Brandy M. Ringham, Keith E. Muller, Deborah H. Glueck

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Department of Biostatistics and Informatics, University of Colorado Denver, Aurora, CO, United States of America

Brandy M. Ringham

Contributed equally to this work with: Sarah M. Kreidler, Brandy M. Ringham, Keith E. Muller, Deborah H. Glueck

Roles Writing – review & editing

* E-mail: [email protected]

Affiliation: LEAD Center, University of Colorado Denver, Aurora, CO, United States of America

ORCID logo https://orcid.org/0000-0001-5677-373X

Keith E. Muller

Contributed equally to this work with: Sarah M. Kreidler, Brandy M. Ringham, Keith E. Muller, Deborah H. Glueck

Roles Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – review & editing

Affiliation: Department of Health Outcomes and Biomedical Informatics, University of Florida, Gainesville, FL, United States of America

Deborah H. Glueck

Contributed equally to this work with: Sarah M. Kreidler, Brandy M. Ringham, Keith E. Muller, Deborah H. Glueck

Roles Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing

Affiliation: Department of Pediatrics, University of Colorado Denver, Aurora, CO, United States of America

Abstract

We derive a noncentral power approximation for the Kenward and Roger test. We use a method of moments approach to form an approximate distribution for the Kenward and Roger scaled Wald statistic, under the alternative. The result depends on the approximate moments of the unscaled Wald statistic. Via Monte Carlo simulation, we demonstrate that the new power approximation is accurate for cluster randomized trials and longitudinal study designs. The method retains accuracy for small sample sizes, even in the presence of missing data. We illustrate the method with a power calculation for an unbalanced group-randomized trial in oral cancer prevention.

Figures

Fig 3

Table 2

Fig 1

Fig 2

Table 1

Fig 3

Table 2

Fig 1

Fig 2

Table 1

Citation: Kreidler SM, Ringham BM, Muller KE, Glueck DH (2021) A power approximation for the Kenward and Roger Wald test in the linear mixed model. PLoS ONE 16(7): e0254811. https://doi.org/10.1371/journal.pone.0254811

Editor: Lei Shi, Yunnan University of Finance and Economics, CHINA

Received: January 7, 2021; Accepted: July 2, 2021; Published: July 21, 2021

Copyright: © 2021 Kreidler et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Source code, data and instructions for reproducing the manuscript results are available at http://github.com/SampleSizeShop/mixedPower.

Funding: This study was supported by The National Institute of Dental and Craniofacial Research (www.nih.gov) in the form of a grant awarded to KEM and DHG (NIDCR 1 R01 DE020832-01A1), The National Institute of General Medical Sciences (www.nih.gov) in the form of a grant awarded to KEM and DHG (NIGMS 9R01GM121081-05), and the Office of the Director (www.nih.gov) in the form of a grant awarded to Dana Dabelea, PI (OD 5UG3OD023248-02). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have read the journal’s policy and have the following competing interests: SMK is an employee of Sunrun, Inc., but was not affiliated with the company at the time the study was conducted. This does not alter our adherence to PLOS ONE policies on sharing data and materials. There are no patents, products in development or marketed products associated with this research to declare.

1 Introduction

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, FR, under the alternative. The reference distribution of FR 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 Notation, models, and hypothesis testing

2.1 Notation

For i ∈ {1, …, n}, let a = {ai} denote an n × 1 column vector. Furthermore, for i ∈ {1, …, n} and j ∈ {1, …, m}, let A = {aij} indicate an n × m matrix with transpose A′ = {aji}. Let Id be a (d × d) identity matrix. For a matrix A = [a1 a2 … an], let . Define the Kronecker product of two matrices A and B as A ⊗ B = {aij B} [9, Section 1.3].

Extend the direct sum operator [9, Section 1.3] to sets of arbitrarily sized matrices as follows. Let {A1, …, AJ} be a set of matrices such that Aj has dimension (rj × cj). Let be an (ri × cj) matrix of zeros. Define the direct sum of {A1, …, AJ} as(1)

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

Let E0(u) and EA(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 Ψ = Σ−1, write to indicate that W−1 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 pi be the number of observations for the ith independent sampling unit, with p = maxi(pi). For the ith independent sampling unit, let yi be the (pi × 1) vector of observed outcomes, Xi be the (pi × r) fixed effects design matrix of rank r, and ei be the (pi × 1) vector of random errors. Assume that for i ≠ j, ei ⊥ ej and yi ⊥ yj. Let Σi be a (pi × pi) 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 ys 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 θ0 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 E0(wR) and from observed data. Kenward and Roger [1] substituted E0(wR) and into method of moments approximations for λ and the reference distribution of FR under the null. With ,(15)(16)and(17)

3 Power approximation for the Kenward-Roger test in the linear mixed model

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, FR, under the alternative. The reference distribution of FR 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.1.1 The conditional distribution of .

The weighted least squares estimate [12] of β is(18)With ,(19)

3.1.2 The approximate distribution of .

We approximate the distribution of(20)with a single central Wishart. The result follows from Theorems 1, 2 and 3 in A. The theorems provide an approximate distribution for a positive definite sum of potentially singular quadratic forms in independent inverse central Wishart matrices.

The accuracy of the approximation depends on the degrees of freedom of the component quadratic forms. To ensure sufficient degrees of freedom, we make the following homoscedasticity assumptions. Recall p = maxi(pi). With Σmax a symmetric, positive definite matrix, assume Σi ≡ Σmax for all i ∈ {1, …, N} such that pi = p. Let Nd indicate the number of independent sampling units with observation pattern Rd. Note . For independent sampling units with observation pattern Rd, assume(21)Without loss of generality, permute the independent sampling units in Eq 8 so that(22)Estimate Σs with(23)

The following thought experiment gives reasonable approximations for the distribution of each . All independent sampling units with observed data pattern Rd have pd observations. For each Rd, suppose we form a complete, balanced mixed model containing only the independent sampling units with observed data pattern Rd. For each balanced mixed model, assume that Xs includes the full time by treatment interaction. This permits recasting each balanced mixed model as an equivalent general linear multivariate model [9, Chapter 14]. For cluster randomized designs, we assume that the mixed model is recast as a two-stage model of cluster means [13, Chapter 4], a special case of the multivariate model.

For the dth multivariate model, let q be the rank of the multivariate design matrix and be the (Nd × pd) matrix of residuals. Assume Nd > (q + pd + 1). Then an unbiased, consistent estimate of Σd, , can be formed using known results for the multivariate model. Thus,(24)with distribution(25)

Recall that in the Wald statistic (Eq 12),(26)Using Eq 25 and Theorem 3 in Appendix, approximate the distribution of with a single inverse central Wishart,(27)From the linear properties of Wishart matrices [11, p. 111, Theorem 3.4.1],(28)

3.1.3 Combining and to form an approximate .

We now combine and as described in Sections 3.1.1 and 3.1.2 to form a Wald statistic,(29)We assume that w ≈ wR. From Eq 19, is approximately Gaussian. From Eq 28, is approximately Wishart.

For conciseness of notation, write μ = (θ − θ0), with estimate , and . Define and . Assume that . The assumption rests on the following logic. If we had estimated both Σs and β using multivariate techniques, independence would follow [14, p. 291, Theorem 8.2.2]. Applying Theorem 4 in Appendix,(30)where(31)and(32)

From Eq 30, we calculate E0(w), EA(w), and , using standard results for central and noncentral distributions [10].

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, FR, 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, FR = λwR, 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 wR.

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 θ0. For i ∈ {1, …, N}, specify Xi and Rd. Calculate a, ν, and γ as described in Section 3.2. Form the reference distribution of . Using the approximate reference distribution of FR under the null, , find the critical value(42)Finally, using the approximate reference distribution of FR under the alternative, , calculate power as(43)

4 Simulation study

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 θ0. For i ∈ {1, …, N}, we specified Xi and Rd. We generated 10, 000 replicates of es and computed ys 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.1.1 Cluster randomized designs.

We compared approximate and empirical power for 36 cluster randomized designs. We assumed that each design had a single Gaussian outcome. Half of the clusters were assumed to have complete data, with the remaining clusters assumed to have some amount of missing data. We varied the number of treatment groups, t ∈ {2, 4}, the number of clusters randomized to each treatment, Ntreatment ∈ {10, 40}, the total number of participants in a complete cluster, p ∈ {5, 50} and the ratio of the incomplete cluster size to the complete cluster size s ∈ {0.6, 0.8, 1}. We only included designs which met the assumption that Nd > (q + pd + 1) for all Rd.

For each design, we repeated the simulations for several intraclass correlation values ρ ∈ {0.04, 0.1, 0.2, 0.5}, with(44)The β matrix had the form(45)for designs with 2 treatments and(46)for designs with 4 treatments. The scale factor b was selected so that the approximate power was roughly 0.2, 0.5 or 0.8. In each scenario, we calculated power for the null hypothesis of no difference among treatment groups at α = 0.05. We used the Wald test with denominator degrees of freedom as described by Kenward and Roger [1].

4.1.2 Longitudinal designs.

We calculated approximate and empirical power for 36 longitudinal study designs. Each design had 5 repeated measures and 50 participants per treatment group. We varied the number of treatment groups, t ∈ {2, 4}, the pattern of missing data, either monotone (missing the 4th and 5th observations), or non-monotone (missing the 2nd and 4th observations), and the number of participants in each treatment group with some amount of missing data, Nincomplete ∈ {0, 10, 20}. For observations within a given participant, we assumed a first-order auto-regressive correlation structure [12, p. 99], with ρ = 0.4 and σ2 = 1. The β matrix had the form(47)for designs with 2 treatments and(48)for designs with 4 treatments. The scale factor and hypothesis testing were as described for the cluster randomized designs with one exception: we calculated power for the null hypothesis of no time by treatment interaction.

4.1.3 Performance criteria.

For each design, we computed the deviation as approximate power minus empirical power. We produced box plots summarizing the deviations overall, within all cluster randomized trials, and within all longitudinal designs. For the cluster randomized trials, we produced box plots stratified by the number of treatment groups, the cluster size, and the ratio of the incomplete cluster size to the complete cluster size. For the longitudinal designs, we produced box plots summarizing the deviations stratified by the number of treatment groups, the pattern of missing observations, and the number of incomplete independent sampling units per treatment.

Positive deviations indicated that the approximate power values were larger than the empirical power values. Negative deviations indicated that the approximate power values were smaller than the empirical power values.

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).

[Figure omitted. See PDF.]

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.

[Figure omitted. See PDF.]

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

[Figure omitted. See PDF.]

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.

[Figure omitted. See PDF.]

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

5 Applied example

We demonstrate a power calculation for an unbalanced cluster-randomized trial of an intervention to reduce oral cancer risk behaviors. The example is based on a hypothetical study examining the impact of workplace smoking cessation programs on tobacco use. We used a synthetic, rather than a real example, so that the power calculation is easy to follow. In a real power calculation, values of differences in means, standard deviations and intra-class correlation coefficients could be drawn from the literature, as described in Guo et al. [17].

For our demonstration, we assume that 80 worksites will be randomized to 2 smoking cessation programs, with 40 sites per treatment condition. Of the 40 sites randomized to each smoking cessation program, 25 worksites will have 30 participants, and the remaining 15 will have 20 participants. The outcome for the analysis will be urinary cotinine. We wish to detect a difference of 25 ng/ml. We assume a standard deviation of 125 ng/ml, and an intraclass correlation of 0.04. We will calculate power for the Kenward and Roger [1] test of the smoking cessation program effect. We set α = 0.05.

To begin the calculation, we first identify the patterns of observations in the study, including complete clusters with 30 participants, and incomplete clusters with 20 participants. Table 2 summarizes the design matrices and patterns of observations by cluster size and treatment assignment.

[Figure omitted. See PDF.]

Table 2. Design matrices and patterns of observations for proposed study of smoking cessation programs.

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

In addition, we define(49)(50)(51)and(52)At an α level of 0.05, the approximate power to detect a treatment difference of 25 ng/ml was 0.87 for the Wald test with Kenward and Roger [1] denominator degrees of freedom.

6 Discussion

We describe a power approximation for the Kenward and Roger (1997) test of fixed effects in the linear mixed model. The method was accurate to within about ±0.06 for all designs, with the best accuracy observed for longitudinal designs. We note that Kenward and Roger (2009) have since described a refinement which improves estimation of the non-linear covariance structures in small samples. We have restricted our discussion to the Kenward and Roger (1997) approach, since it is most commonly used in statistical practice.

The method has several limitations. The assumption of Nd > (q + pd + 1) may be too restrictive for multilevel designs with large cluster sizes. In addition, we assume that the pattern of missing data is known. The method does not apply to repeated covariates, which often appear in biomedical studies. However, the method does apply to baseline covariates, a common study design. We make a strong homoscedasticity assumption of equal variance for each independent sampling unit. This assumption means that the power computations are not appropriate for random regression, for models with group differences in variance, or for certain spatial-temporal applications. Nevertheless, the assumption of homoscedasticity is widely made for randomized controlled clinical trials, laboratory studies, and observational studies, which makes the method useful for a variety of cases. Lastly, the method has not been evaluated for binary or Poisson data.

The analytic results from this manuscript suggest several future extensions. We may be able to calculate power for linear mixed models with random missing data patterns by invoking conditional distribution theory and calculating expected power across patterns of missingness. In addition, the approach used to form the distribution of provides the first step towards a non-iterative alternative to restricted maximum likelihood estimation for some mixed models. For big data applications, such a non-iterative approach may facilitate highly parallel computation of parameter estimates in mixed models.

Our power approximation provides a general, flexible, accurate and rapid method to calculate power for the Kenward and Roger (1997) test. For studies in which the Kenward and Roger (1997) test is the planned method of data analysis, our power approximation should be used. By aligning power analysis with the planned data analysis, researchers can more accurately assess power for biomedical studies. Accurate power analysis is an ethical imperative for research with human participants.

7 Appendix

A Appendix: Theorems and proofs

Theorem 1. For m ∈ {1, …, k}, let pm ∈ {1, 2, …, p}, Nm > (pm + 3) and define Ψm = {ψmij} to be a (pm × pm) 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 Rm ⊂ {1, 2, …, p} of cardinality pm, define Xm to be a (pm × qp) matrix of rank pm < qp with the form (53) If for each i ∈ {1, …, q}, there exists at least one m such that Xm = Iq({i}) ⊗ Ip, then (54) is positive definite.

Proof. Let Qi = {Xm: Xm = Iq({i}) ⊗ Ip(Rm)}. Then(55)Note that for i ∈ {1, 2, ‥, q}, Iq({i})′ Iq({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 Qi, there exists a ci 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}, Rm ⊂ {1, 2, …, p} of cardinality pm, Xm = Iq({i}) ⊗ Ip(Rm) a (pm × qp) matrix of rank pm < qp, Nm > (pm + 3), Ψm a (pm × pm) 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 (pm × 1) vector of eigenvalues and V is the (pm × pm) orthogonal matrix of eigenvectors of . Then(61)

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

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

Since λ0 contains only zero eigenvalues and using the definition of the trace,(63)Theorem 3. For m ∈ {1, …, k}, let pm ∈ {1, …, p}, Nm > (pm+ 3) and let Ψm = {ψmij} be a (pm × pm) 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 Rm ⊂ {1, …, p} of cardinality pm, define Xm to be a (pm × qp) matrix of rank pm < qp with the form (64) Under the assumption that for each i ∈ {1, …, q}, there exists at least one m such that Xm = Iq({i}) ⊗ Ip, it can be shown that (65) is approximately distributed as .

Proof. Theorem 1 in Appendix demonstrates that Q−1 is positive definite under the restriction that for each i ∈ {1, …, q}, there exists at least one m such that Xm = Iq({i}) ⊗ Ip.

To derive an approximate distribution for Q−1, 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 χ2 random variables [9, Theorem 9.5, p. 176]. Approximate the distribution of U with a single noncentral χ2, so that . Using the approach described by Kim et al. [8], obtain values for λu, nu 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−1 x,and the result follows.

Acknowledgments

A portion of this paper was submitted to the University of Colorado Denver in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Biostatistics for Dr. Sarah M. Kreidler.

Citation: Kreidler SM, Ringham BM, Muller KE, Glueck DH (2021) A power approximation for the Kenward and Roger Wald test in the linear mixed model. PLoS ONE 16(7): e0254811. https://doi.org/10.1371/journal.pone.0254811

References

1. Kenward MG, Roger JH. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics. 1997;53(3):983–997. pmid:9333350

2. Helms RW. Intentionally incomplete longitudinal designs: I. Methodology and comparison of some full span designs. Statistics in medicine. 1992;11(14-15):1889–1913. pmid:1480880

3. Stroup WW. Mixed Model Procedures to Assess Power, Precision, and Sample Size in the Design of Experiments. 1999 Proceedings of the Biopharmaceutical Section, Alexandria, VA: American Statistical Association. 1999; p. 15–24.

4. Tu XM, Kowalski J, Zhang J, Lynch KG, Crits-Christoph P. Power analyses for longitudinal trials and other clustered designs. Statistics in medicine. 2004;23(18):2799–2815. pmid:15344187

5. Tu XM, Zhang J, Kowalski J, Shults J, Feng C, Sun W, et al. Power analyses for longitudinal study designs with missing data. Statistics in medicine. 2007;26(15):2958–2981. pmid:17154250

6. Shieh G. A unified approach to power calculation and sample size determination for random regression models. Psychometrika. 2007;72(3):347–360.

7. Chi YY, Glueck DH, Muller KE. Power and Sample Size for Fixed-Effects Inference in Reversible Linear Mixed Models. The American Statistician. in press;. pmid:32042203

8. Kim HY, Gribbin MJ, Muller KE, Taylor DJ. Analytic, Computational, and Approximate Forms for Ratios of Noncentral and Central Gaussian Quadratic Forms. Journal of computational and graphical statistics: a joint publication of American Statistical Association, Institute of Mathematical Statistics, Interface Foundation of North America. 2006;15(2):443–459. pmid:23843686

9. Muller KE, Stewart PW. Linear model theory: univariate, multivariate, and mixed models. Hoboken, New Jersey: John Wiley and Sons; 2006.

10. Johnson NL, Kotz S, Balakrishnan N. Continuous univariate distributions. Wiley & Sons; 1995.

11. Gupta AK, Nagar DK. Matrix variate distributions. Boca Raton, FL: Chapman & Hall; 2000.

12. Verbeke G, Molenberghs G. Linear mixed models for longitudinal data. New York: Springer; 2009.

13. Murray DM. Design and Analysis of Group- Randomized Trials. 1st ed. Oxford University Press, USA; 1998.

14. Anderson TW. An Introduction to Multivariate Statistical Analysis. 2nd ed. Wiley Series in Probability and Statistics. Wiley; 1984.

15. R Development Core. R: A Language and Environment for Statistical Computing. Vienna, Austria; 2010. Available from: http://www.R-project.org/.

16. SAS Institute Inc. SAS 9.3 Software, Version 9.3. Cary, NC; 2013. Available from: http://www.sas.com/software/sas9/.

17. Guo Y, Logan HL, Glueck DH, Muller KE. Selecting a sample size for studies with repeated measures. BMC Medical Research Methodology. 2013;13(1). pmid:23902644

18. Mathai AM, Provost SB. Quadratic Forms in Random Variables: Theory and Applications. Marcel Dekker Incorporated; 1992.

19. Arnold SF. The theory of linear models and multivariate analysis. New York: Wiley; 1981.