Design and analysis considerations for cohort stepped wedge cluster randomized trials with a decay correlation structure

3.1 |. Population-averaged models

The population-averaged model relates the marginal mean, μ_ijt, to the time trend and the intervention effect by

where g is the link function and β_t is the tth period effect.⁵ Furthermore, X_it is the intervention status, which is equal to 1 or 0 depending on whether cluster i receives intervention during period t, and δ describes the intervention effect on the link function scale. Model (1) can be regarded as the marginal counterpart of a number of existing random-effects models, such as those proposed in the works of Hemming et al,⁷ Hooper et al⁸ and Kasza et al.^10,11 Like these random-effects models, our marginal model (1) does not specify treatment by period interaction and so δ should be interpreted as the average intervention effect across periods. We write the collection of model parameters as θ = (β₁, …, β, δ)′, the collection of intervention status for cluster i (a sequence of ones preceded by zeros) as X_i = (X_i1, …, X_iT)′, and define v(μ_ijt) as a known variance function. To allow for potential correlation decay over time, we define the proportional decay correlation structure similar to the work of Lefkopoulou et al.¹⁶ Specifically, we define the within-period correlation as the correlation between outcomes for two distinct individuals from the same cluster during the same period, ie, corr(y_ijt, y_ij′t) = τ for j ≠ j′. The same definition is prevalently used in the posttest-only parallel designs.¹ We then assume a first-order autoregressive structure for the set of within-individual repeated measurements. Therefore, the within-individual correlation, which describes the association between outcomes measured at time t and t′ of the same individual, is corr(y_ijt, y_ijt′) = ρ^{|t − t′|}, t ≠ t′. Finally, we define the between-period correlation as the correlation between outcome measured at time t for individual j and outcome measured at time t′ for individual j′, and assumes a decay structure as corr(y_ijt, y_ij′t′) = τρ^|t−t′| for j ≠ j′, t ≠ t′.

In a closed-cohort design, typically, the between-period correlation is much smaller than the within-individual correlation for any two fixed periods, since the former is defined for measurements from two distinct individuals, while the latter is defined for those from the same individual. This ordering of magnitude is reflected in the proportional decay structure because τ < 1 and corr(y_ijt, y_ij′t′) = τρ^|t−t′| < corr(y_ijt, y_ijt′) = ρ^|t−t′|. By comparison, the exponential decay structure by Kasza et al¹⁰ is based on cross-sectional designs and therefore obviates the need for modeling the within-individual correlation from repeated measurements. In fact, the exponential decay structure assumes corr(y_ijt, y_ij′t′) = τρ^|t−t′| irrespective of whether j = j′ since two different sets of participants are included in two different periods for a cross-sectional design. In this respect, parameters τ and ρ have similar interpretations in both structures. We provide a visual comparison of these two structures in Table 2. In summary, the proportional decay correlation structure is defined through two parameters, τ and ρ, with the former resembling the traditional ICC definition in a parallel design and the latter controlling for the degree of correlation decay. Of note, the works Shults and Morrow¹⁷ and Liu et al¹⁸ also adopted the same proportional decay structure in longitudinal CRTs with a parallel assignment, and we extend the discussion of this decay structure to CRTs with a staggered assignment.

TABLE 2.

Examples of the proportional decay structure for cohort designs and the exponential decay structure for cross-sectional designs. The illustration is based on a stepped wedge trial with T = 3 periods and N_i = 2 measurements per cluster period. Define y_i = (y_i11, y_i12, y_i13, y_i21, y_i22, y_i23)′

	Proportional decay structure	Exponential decay structure
corr(yi)

Open in a new tab

3.2 |. Quasi-least squares analysis

We use the QLS approach introduced by Shults and Morrow¹⁷ to simultaneously estimate the intervention effect in model (1) and the correlation parameters in a cohort stepped wedge design. In particular, the QLS approach and the traditional GEE approach share the same estimating equations for the marginal mean parameters, whereas the former provides a flexible and convenient way to estimate nonstandard correlation structures. Furthermore, both the QLS and GEE approaches are robust to correlation misspecification, namely, estimators for the marginal intervention effect remain consistent even if the working correlation model is misspecified. In sufficiently large samples (usually I ≥ 30), the robust sandwich variance could also be used to adequately quantify the uncertainty of the intervention effect estimate even under correlation misspecification. We refer the readers to the textbook of Shults and Hilbe²³ for a full exposition on the advantages of QLS over the traditional GEE.

Write y_ij = (y_ij1, …, y_ijT)′, μ_ij = (μ_ij1, …, μ_ijT)′, $y_i={\left(y_{i1}^{\prime },\dots ,y_{i{N}_i}^{\prime }\right)}^{\prime }$, and ${\mu }_i={\left({\mu }_{i1}^{\prime },\dots ,{\mu }_{i{N}_i}^{\prime }\right)}^{\prime }$. Furthermore, define D_i(θ) =∂μ_i ∕∂θ′, and let the working covariance of y_i be $V_i=\varphi A_i^{1/2}(\theta )R_i\left({\alpha }_0,{\alpha }_1\right)A_i^{1/2}(\theta )$, where ϕ > 0 is the dispersion parameter, $A_i(\theta )=\text{diag}\left\{A_{i1}(\theta ),\dots ,A_{i{N}_i}(\theta )\right\}$, A_ij(θ) = diag{v(μ_ij1), …, v(μ_ijT)}, and R_i(α₀, α₁) is a positive definite working correlation parameterized by α₀ and α₁. We assume the true correlation structure among elements of y_i is the proportional decay structure, denoted by R_i(τ, ρ). In matrix notation, we can verify that the proportional decay structure induces separability between τ and ρ in that R_i(τ, ρ) = G_i(τ) ⊗ F(ρ), where G_i(τ) is a N_i × N_i exchangeable correlation matrix (with respect to τ) and F(ρ) is a T × T first-order autoregressive correlation matrix as follows:

We could verify that the determinant

Therefore, valid correlation values that ensure positive definite R_i(τ, ρ) should be contained in the triangular region

Finally, the inverse of the R_i also exists in closed form and is given by $R_i^{-1}(\tau ,\rho )=G_i^{-1}(\tau )\otimes F_i^{-1}(\rho )$, where

I_i is a N_i × N_i identity matrix, J_i is a N_i × N_i matrix of ones, C₂ = diag(0, 1, …, 1, 0), and C₁ is a T × T tridiagonal matrix with zeros on the main diagonal and ones on the two subdiagonals.

To introduce the QLS estimating equations, we further define $r_{ij}(\theta )=A_{ij}^{-1/2}(\theta )\left(y_{ij}-{\mu }_{ij}\right)$, and write $r_i(\theta )={\left(r_{i1}^{\prime }(\theta ),\dots ,r_{i{N}_i}^{\prime }(\theta )\right)}^{\prime }$. The first-stage QLS estimates for θ, α₀, and α₁ are obtained by alternating between the following estimating equations until convergence:

In particular, (3) is the usual GEE coupled with the proportional decay structure, and (4) and (5) are scalar equations for the first-stage correlation estimates. Further, closed-form solutions exist for ${\widehat{\alpha }}_0$ and ${\widehat{\alpha }}_1$ (within an iterative step) and are provided in Web Appendix A. Chaganty and Shults²¹ showed that ${\widehat{\alpha }}_0$ and ${\widehat{\alpha }}_1$ are asymptotically biased for τ and ρ. To eliminate the large-sample bias in the first-stage correlation estimates, Chaganty and Shults²¹ provided the following second-stage estimating equations to obtain: $\widehat{\tau }$, $\widehat{\rho }$

The closed-form solution for (6) and (7) are provided by the work of Shults and Morrow¹⁷ as

and $\widehat{\rho }=2{\widehat{\alpha }}_1/\left(1+{\widehat{\alpha }}_1^2\right)$.

3.3|. Bias-corrected correlation estimation

Although the correlation estimates obtained from the second-stage QLS estimating equations are unbiased in large samples, they could be subject to finite-sample bias. The finite-sample bias is a typical consideration in stepped wedge CRTs,⁵ as they usually include a small number of clusters (I ≤ 30). We refine the QLS approach with finite-sample bias corrections to the correlation estimating equations by utilizing the cluster leverage,²⁴ defined as $H_i=D_i(\theta ){\left({\sum }_{i=1}^I{D}_i^{\prime }(\theta )V_i{D}_i(\theta )\right)}^{-1}{D}_i^{\prime }(\theta )V_i$. Specifically, when I is small, the estimated residual $y_i-{\widehat{\mu }}_i$ tends to be biased toward zero, and following the work of Preisser et al,²⁵ we have $E\left[\left(y_i-{\widehat{\mu }}_i\right){\left(y_i-{\widehat{\mu }}_i\right)}^{\prime }\right]\approx \left(I-H_i\right)\text{cov}\left(y_i\right)=\varphi \left(I-H_i\right)A_i^{1/2}\text{corr}\left(y_i\right)A_i^{1/2}$ and, therefore, $E\left[r_i(\widehat{\theta })r_i^{\prime }(\widehat{\theta })\right]\approx \varphi A_i^{-1/2}\left(I-H_i\right)A_i^{1/2}\text{corr}\left(y_i\right)$. This last equation suggests that ${\varphi }^{-1}{A}_i^{-1/2}{\left(I-H_i\right)}^{-1}{A}_i^{1/2}{r}_i(\widehat{\theta })r_i^{\prime }(\widehat{\theta })$ is a better estimator for corr(y_i) compared to the simple cross-product ${\varphi }^{-1}{r}_i(\widehat{\theta })r_i^{\prime }(\widehat{\theta })$, since the former accounts for finite-sample bias in a multiplicative fashion. Furthermore, observe that

for all values of α₀, α₁, and θ; we propose to replace the first-stage estimating equations (4) and (5) by

where ${\tilde{R}}_i(\theta )={\varphi }^{-1}{A}_i^{-1/2}{\left(I-H_i\right)}^{-1}{A}_i^{1/2}{r}_i(\theta )r_i^{\prime }(\theta )$ represents the matrix-adjusted estimator for the correlation structure. The solutions obtained from (8) and (9) could effectively reduce the finite-sample bias in ${\widehat{\alpha }}_0$ and ${\widehat{\alpha }}_1$, which would in turn decrease the finite-sample bias in the QLS estimators for τ and ρ. Of note, similar finite-sample matrix adjustment was developed by Preisser et al²⁵ for the Prentice-type GEE,²⁶ and we extend this finite-sample bias-correction approach to the QLS estimating equations. Accurately estimating the correlation parameters in the analysis stage has practical implications since these estimates could be used to guide the planning of future trials.¹ Additional details of the matrix-adjusted estimating equations, (8) and (9), along with the closed-form updates are provided in Web Appendix B.

3.4 |. Bias-corrected covariance estimation

The availability of a small number of clusters may also have implications for estimating the variance using GEE-based approaches.²⁷ In general, the variance of the marginal mean model parameter $\widehat{\theta }$ can be estimated using the model-based variance ${\mathrm{\Omega }}_1^{-1}={\left({\sum }_{i=1}^I{D}_i^{\prime }{V}_i{D}_i\right)}^{-1}$ or the sandwich variance ${\mathrm{\Omega }}_1^{-1}{\mathrm{\Omega }}_0{\mathrm{\Omega }}_1^{-1}$, where

and both Ω₀ and Ω₁ are evaluated at ($\widehat{\theta }$, $\widehat{\tau }$, $\widehat{\rho }$). When both C_i and B_i are identity matrices, Equation (10) reduces to the uncorrected sandwich estimator of Liang and Zeger,¹⁹ which we denote as BC0. BC0 provides valid inference regardless of the correct specification of the working correlation R_i, as long as the number of clusters is sufficiently large (I ≥ 30), while the consistency of the model-based variance requires the correct specification of the correlation structure. As $r_i(\widehat{\theta })$ is biased toward zero with a limited number of clusters, BC0 is likely to underestimate the variance and alternative choices of matrices C_i and B_i maybe necessary to provide a partial correction to the finite-sample bias.²⁷ We consider three popular approaches for bias corrections summarized as follows and also in Table 3: the finite-sample correction due to the work of Kauermann and Carroll,²⁸ or BC1, is given by C_i = I and B_i = (I − H_i)^−1/2; the finite-sample correction due to the work of Mancl and DeRouen,²⁹ or BC2, is given by C_i = I and B_i = (I − H_i)⁻¹; and the finite-sample correction due to the work of Fay and Graubard,³⁰ or BC3, given by $C_i=\text{diag}\left\{{\left(1-\text{min}\left\{\zeta ,{\left[D_i^{\prime }{V}_i^{-1}{D}_i{\mathrm{\Omega }}_1^{-1}\right]}_{jj}\right\}\right)}^{-1/2}\right\}$ and B_i = I, where the bound parameter ζ is a user-defined constant (< 1) with a default value 0.75. Because the matrix elements of the cluster leverage are between 0 and 1, we generally have BC0 < BC1 < BC2 in terms of the amount of correction.²⁵ Furthermore, Scott et al³¹ have shown that BC3 tends to be close to BC1. Of note, the estimation of dispersion parameter should only affect the model-based variance. Similar to the work of Liang and Zeger,¹⁹ we propose to consistently update the dispersion parameter from iteration s to s + 1 by ${\widehat{\varphi }}^{(s+1)}={\widehat{\varphi }}^{(s)}{\sum }_{i=1}^I\text{tr}\left({\tilde{R}}_i\right)/\left\{\sum _{i=1}^{I}T{N}_i-(T+1)\right\}$.

4.1 |. Closed-form variance of the intervention effect

Under the null hypothesis H₀: δ = δ₀, the large-sample variance of $\sqrt{n}\left(\widehat{\delta }-{\delta }_0\right)$ is provided by the (T + 1, T + 1) element of the large-sample covariance matrix of $\sqrt{n}\left(\widehat{\theta }-{\theta }_0\right)$. Since the QLS estimator $\widehat{\delta }$ is asymptotically normal, we could use the z-test statistic $\widehat{\delta }/\sqrt{\text{var}(\widehat{\delta })}$ to test the null of no intervention effect H₀: δ = 0, and the power to detect an intervention effect of size δ ≠ 0 with a prescribed type I error rate α is approximately

where Φ is the standard normal cumulative distribution function and z_α/2 is the normal quantile such that Φ(z_α/2) = 1 − α∕2. Because there is uncertainty in estimating the asymptotic variance $\text{var}(\widehat{\delta })$_, an alternative two-sided test uses the same statistic but refers to the t-distribution. We consider two choices of degrees of freedom (DoF). The first DoF dates back to the work of Mancl and DeRouen²⁹ and equal to the number of clusters minus the number of regression parameters; this DoF has been previously used in the GEE analyses of parallel CRTs,³² three-level CRTs,³³ crossover CRTs,⁶ stepped wedge CRTs,^5,34 and shown to have test size not exceeding the nominal level. The second DoF was suggested in the PhD dissertation of Ford³⁵ and specifies DoF = I − 2. This choice of DoF was found to provide excellent control of type I error rate for GEE analyses of both parallel CRTs and stepped wedge CRTs in the work of Ford and Westgate.^35,36 With the same effect size δ and prescribed type I error rate α, the power of the t-test is approximately

where Ψ_DoF is the t distribution function and the quantile t_α/2 is chosen such that Ψ_DoF(t_α/2) = 1 − α∕2. We notice that, because the t-distribution has a heavier tail compared with the standard normal distribution, the QLS z-test ismore likely to result in an inflated type I error rate with the use of BC0 than the corresponding QLS t-test. As the bias-corrected sandwich variance estimators (BC1, BC2, and BC3) provide different degrees of inflation relative to the uncorrected variance BC0, an investigation of both the z- and t-tests coupled with the collection of alternative variance estimators could help inform the practical choice among the analytical options for stepped wedge CRTs.

To assist the design of stepped wedge trials allowing for correlation decay, we derive a new closed-form variance expression for $\widehat{\delta }$ assuming the outcome is continuous and g is the identity link function. We will return to categorical outcomes and nonlinear link functions in Section 7. To do so, we follow the work of Shih³⁷ and assume the covariance of Y_i to be known as var(Y_i) = Vi. Therefore, $\text{var}(\widehat{\delta })$ is the (T+1, T+1) element of themodel-based variance ${\mathrm{\Omega }}_1^{-1}$. We further assume a balanced design such that an equal number of participants will be recruited in each cluster prior to the first period, so that N_i = N. Such a simplification assumption is routinely made in designing stepped wedge trials.

Under a balanced design, we could write the design matrix corresponding to cluster i as Zi = 1N ⊗ (I_T, X_i), where 1_N is a N-vector of ones. Then, the variance of the intervention effect $\widehat{\delta }$ is equal to the lower-right element of $\varphi {\left\{{\sum }_{i=1}^I{Z}_i^{\prime }{R}_i^{-1}(\tau ,\rho )Z_i\right\}}^{-1}$, where ϕ is the marginal variance. We show in Web Appendix C that a closed-form variance expression for $\widehat{\delta }$ is

where the design constant $U={\sum }_{i=1}^I{\sum }_{t=1}^T{X}_{it}$ is the total number of cluster periods exposed under the intervention condition, $W={{\sum }_{t=1}^T\left({\sum }_{i=1}^I{X}_{it}\right)}^2$ is the squared number of clusters receiving the intervention summed across periods, and $V={\sum }_{i=1}^I{\sum }_{t=1}^T{X}_{it}{X}_{i,t+1}$ and $Q={\sum }_{t=1}^{T-1}\left({\sum }_{i=1}^I{X}_{it}\right)\left({\sum }_{i=1}^I{X}_{i,t+1}\right)$ are cross-product terms resulting from the decay correlation structure. It is interesting to see that this variance expression does not depend on the magnitude of the period effect β_t as long as they are controlled for in the marginal mean model. Noticeably, the QLS-based variance (13) extends the formula due to the work of Liu et al¹⁸ to longitudinal cluster designs with staggered randomization. Furthermore, as the cohort size N becomes large, the variance expression converges to

which is a finite constant since |ρ| < 1 and τ > 0 for large N, according to (2). Therefore, the limit of the variance is a positive constant determined by available design resources I, T, and two correlation values, τ and ρ, and cannot be made arbitrarily small. In other words, the power of the stepped wedge design may not be increased to one by solely increasing the cohort size, which is consistent with the known results for parallel cluster randomized designs.¹ For this reason, when N is large, variance (14) could be used in the design stage to approximate the variance (13). Finally, given hypothesized values for I, N, T, and correlation parameters, variance expression (13) or (14) can be used in Equations (11) and (12) to obtain the predicted power.

4.2 |. The design effect

For determining the required sample size based on Equations (11) and (12), it is straightforward to solve N by fixing the required number of clusters I but not the other way around. However, with a predetermined cohort size N for each cluster, we could postulate a series of values for I and find the smallest value such that the estimated power is at least equal to the prescribed level. Additionally, in the following case studied by Woertman et al,²² we could derive a simple expression for the design effect (DE) relative to an individually randomized trial to simplify sample size calculation. Specifically, we assume that an equal number of clusters switch to intervention at each step so that m_s = m, and, furthermore, an equal number of measurements are taken between steps such that c_s = s for all s = 1, …, S. We then write the total number of clusters I = Sm and total number of periods T = b + Sc, and the design constants become

Plugging the design constants back into the variance formula (13) and dividing by the variance of the two-sample mean difference 4ϕ∕(NSm), we obtain

The aforementioned DE allows us to easily study how the design resources affect the statistical efficiency relative to individual randomization and how the correlation parameters affect the statistical power. For example, since the DE is free of b, the relative design efficiency does not change according to the number of baseline periods. However, for fixed values of the correlation parameters, increasing the number of steps S and number of measurements between steps c decreases the DE and increases the efficiency. On the other hand, for fixed design resources, larger values of the within-period correlation τ increases the DE, confirming that τ functions as the traditional ICC of a parallel CRT. By contrast, the role of correlation parameter ρ is characterized by f(p) = (1 − ρ²)∕[(S + 1)c(1 − ρ)² + 6ρ], which is monotonically increasing on (−1, r) and decreasing on (r, 1), where

For convenience, we can also define the decay parameter d = 1 − ρ so that d = 0 and d = 1 correspond to no decay and total decay, respectively. Since it is more plausible that ρ ∈ (0, 1), the aforementioned result suggests that, with an increasing level of correlation decay, the DE first increases to its largest value and then decreases, with the maximum DE obtained at ρ = r. A numerical illustration of the DE as a function of the decay parameter is provided in Web Figure 2.

5.1 |. Simulation design

We carry out a simulation study (1) to compare the correlation estimates from the uncorrected QLS and the proposed matrix-adjusted QLS (MAQLS), and (2) to evaluate the utility of the proposed power formula for QLS-based analyses of stepped wedge CRTs. For the second objective, we first determine the empirical type I error rates for the QLS-based tests coupled with alternative variance estimators, and then identify valid tests (those with a close-to-nominal type I error rate) whose empirical power corresponds well with the predicted power from the proposed formula. Findings specific to the second objective are informative for practical data analysis since we prefer tests that maintain a valid size and meanwhile demonstrate empirical power that is at least the magnitude of the analytical prediction.

Within-cluster correlated continuous outcomes were generated from a multivariate normal distribution with mean given by μ_ijt = β_t+X_itδ and covariance ϕR(τ, ρ), where R(τ, ρ) is the proportional decay structure defined in Section 3.1. We set the marginal variance ϕ = 1 and assumed a gently increasing period effect such that β₁ = 0 and β_t+1−β_t = 0.1×(0.5)^t−1 for t ≥ 1. As discussed before, the predicted power should be insensitive to the magnitude of the period effects as long as they are accounted for in the QLS analyses. We fix the effect size δ∕ϕ^1/2 at zero for studying test size and choose δ∕ϕ^1/2 from {0.2, 0.3, 0.4, 0.5} for studying power. We choose the within-period correlation τ ∈ {0.03, 0.1}, which represent typical ICC values reported in the parallel CRTs.¹ We further chose ρ ∈ {0.2, 0.8}, representing large and moderate degree of correlation decay over time. The number of clusters is varied from 9 to 24 as stepped wedge CRTs usually include a limited number of clusters. We specify the number of periods 3 ≤ T ≤ 8 as these values are frequently reported in practice, according to recent reviews by the works of Martin et al³⁸ and Grayling et al.³⁹ Finally, the cohort size is chosen as 5 ≤ N ≤ 24 to ensure that the predicted power is at least 80%. For illustration, we focus on standard stepped wedge designs so that there is only one baseline period and the number of steps S = T − 1. In other words, an equal number of I∕S clusters cross over to intervention during each step, and the outcome is measured only once for each individual between consecutive steps. For each scenario, 10 000 data replications were generated and analyzed using both QLS and MAQLS. For the first objective, we report the percent relative bias in estimating τ and ρ. In general, an unbiased approach for estimating the correlation parameters is preferred since accurate reporting of correlations is critical for planning future trials. Web Tables 1 and 2 provide a summary of the simulation scenarios along with the convergence rates. The convergence rates are similar between QLS and MAQLS, and all exceed 97%. For the second objective, we consider both the z-tests and the t-tests for testing the null hypothesis of no intervention effect, coupled with five different variance estimators for $\widehat{\delta }$, namely, the model-based variance, BC0, BC1, BC2, and BC3. The nominal type I error rate is held fixed at 5%, and we consider an empirical type I error rate between 4.5% and 5.5% to be acceptable based on the margin of error from a binomial model with 10 000 replications. By a similar reasoning, since the predicted power in each scenario is at least 80%, we consider an empirical power that differs by no more than 0.8% from the predicted value to be acceptable.

5.2 |. Results

For the first objective, we summarize in Table 4 and Web Table 3 the percent relative bias in estimating the correlations with QLS and MAQLS. It is evident that the percent bias in estimating the within-period correlation τ is much larger than that in estimating the correlation parameter ρ, without respect to the incorporation of matrix adjustment to the first-stage estimating equations (4) and (5). However, the QLS estimator for τ exhibits noticeable downward bias, especially when the number of clusters is not large. By contrast, MAQLS substantially reduces such finite-sample bias and improves the estimation of τ. On the other hand, the QLS estimator for the parameter ρ seems more accurate in that the absolute percent bias only occasionally exceeds one. Nevertheless, MAQLS still mildly improves the estimation of ρ in that the absolute percent bias is always maintained under one. The comparative findings between QLS and MAQLS are consistent regardless of the magnitude of intervention effect δ. Therefore, MAQLS is the preferred approach because it provides much less biased estimates for the correlation parameters; these more accurate correlation estimates will eventually facilitate accurate estimation of sample size and power for future cohort stepped wedge trials.

For the second objective, we present the empirical type I error rates of the z-tests and t-tests for the QLS and MAQLS analyses in Web Tables 4 to 6 and Web Figure 3. Overall, we observe that the matrix adjustment to the correlation estimation mildly affects the tests with the model-based variance but has little impact on the tests with the sandwich variance. This is in accordance with the work of Lu et al,⁴⁰ who observed the same results for the GEE analyses of pretest-posttest CRTs. Since MAQLS provides more accurate estimation of the correlations, we will focus on this approach. Table 5 summarizes the empirical type I error rates of the MAQLS z-tests and t-tests with DoF = I − 2. We leave the results for t-tests with DoF = I − (T + 1) to Web Table 6 as these tests are conservative in many cases. From Table 5, we observe that MAQLS z-tests are more liberal than the corresponding MAQLS t-tests. The type I error rates of the MAQLS z-tests coupled with the model-based variance or BC2 are close to nominal when I ≥ 20, while the MAQLS z-tests coupled with BC0, BC1, or BC3 are always liberal. By contrast, only the MAQLS t-tests with BC0 remain liberal, while MAQLS t-tests with BC1 or BC3 maintain close-to-nominal size and MAQLS t-tests with model-based variance or BC2 are conservative. Overall, the t-tests with DoF = I − 2 and BC1 demonstrate test sizes that consistently agree with the nominal level.

Web Tables 7 to 9 and Web Figure 4 present the predicted and empirical power for each simulation scenario. Because we are only interested in tests that maintain close-to-nominal sizes, we summarize in Table 6 the differences between the empirical and the predicted power only for the MAQLS z-tests and t-tests with DoF = I − 2. Among the z-tests, only the choice of BC2 provides substantially lower power than predicted. While the choices of model-based variance BC1 and BC3 provide adequate power for the z-tests in a number of scenarios, one should be cautious in adopting these tests with a small number of clusters since they may carry an inflated test size. On the other hand, the empirical power for MAQLS t-tests coupled with the model-based variance, BC1 or BC3 corresponds reasonably well with the analytical prediction from the proposed formula even for as few as nine clusters, while the empirical power for MAQLS t-tests with BC2 still tends to be substantially lower than predicted. Interestingly, the MAQLS t-test with DoF = I − (T + 1) coupled with model-based variance, BC1 or BC3 also demonstrates empirical power fairly close to prediction, even though these tests are more conservative under the null. These tests may not be preferred over the t-tests with DoF = I − 2 since they are less powerful.

6.1 |. The AEP study

We illustrate the proposed sample size procedure to design a cohort stepped wedge CRT that aims to study the effect of an exercise intervention on the physical function of patients with end-stage renal disease.⁴¹ The intervention was an accredited exercise physiologist (AEP) coordinated resistance exercise program, offered at hemodialysis clinics to improve the quality of life for dialysis patients. During the planning phase, it was determined that I = 15 clinics (clusters) were available, and would be randomized over T = 4 periods evenly spaced across 48 weeks. At baseline, no exercise programs were offered to any clinic. At week 12, 36 and 48, a random subset of 5 clinics cross over from control to intervention. A closed cohort of patients was recruited at baseline, and would be followed up during the study period. The primary patient-level outcome was the 30-second sit-to-stand test, recording the number of times a patient could rise from and return to a seated position in a 30-second time frame. The 30-second sit-to-stand test was conducted at the end of each period, resulting in four outcome measurements per patient. Based on a prior study within a similar context, a conservative estimate of the effect size was given by δ∕ϕ = 0.325,⁴² and the within-period correlation was estimated to be τ = 0.03.⁴³ With I = 15 clusters, the simulations suggest the MAQLS t-test with DoF = I − 2 = 13 could maintain nominal size and adequate power; we illustrate the power calculation based on the t-test statistic.

Given this is a standard stepped wedge design where an equal number of clinics switch to intervention at each step, we can show that U = IT∕2, W = I²T(2T − 1)∕{6(T − 1)}, V = I(T − 2)∕2 and Q = I²T(T − 2)∕{3(T − 1)}. The variance expression (13) is then simplified to

If we anticipate large correlation decay so that ρ = 0.2, the power is estimated using Equations (11) and (16) to be 79.4% if N = 21 and 80.5% if N = 22. Therefore, at least N = 22 patients should be recruited in each clinic to achieve 80% power under the proportional decay structure. On the other hand, we could arrive at the same results by using the DE formula (15). For example, in an individual randomized study, 348 patients would be required for the hypothesized effect size. Assuming 21 patients will be included in each clinic, the DE is approximately 0.92, indicating a total of 320 patients in approximately 15.2 clinics would be required. Since the study affords to randomize only 15 clinics, we increase the cohort size to N = 22, resulting in a DE 0.94. Therefore, 326 patients are required for a total of 326∕22 ≈ 14.8 clinics, and we conclude that 22 patients in 15 clinics ensured 80% power.

While the within-period correlation estimate was available from prior studies, published estimates of the correlation parameter ρ (or decay parameter d = 1 − ρ) are currently rare. For this reason, we carry out a sensitivity analysis on the power and present the results in panel (A) of Figure 2, where we fix the design resources but vary τ ∈ (0.03, 0.06) and ρ ∈ (0, 1). Note that the upper bound of the within-period correlation τ was reported by the work of Littenberg and MacLean⁴³ and is used in this assessment. As expected, larger values of the within-period correlation reduce the study power, and furthermore, given a certain value of the within-period correlation, a greater magnitude of decay (smaller ρ or larger d) generally reduces the study power unless ρ ≈ 0 (or d ≈ 1). For the hypothesized τ = 0.03, the study power remains at least close to 80% regardless of the correlation decay. On the other hand, the amount of correlation decay could result in further power loss if the within-period correlation τ increases. Nevertheless, the power loss is at most around 10% even if the within-period correlation τ approximates the upper bound 0.06.

6.2 |. The CORE study

We next illustrate the proposed sample size procedure by designing the CORE stepped wedge trial. The CORE study aims to evaluate the patient-centered service design in health providers to improve the psychosocial recovery outcomes for people with severe mental illness in Australia.⁴⁴ The new service design intervention adopted the experience-based co-design to identify users’ positive and negative experiences of the service, and involved patients’ participation to co-design solutions to the negative experiences. A total of I = 11 teams from four health service providers would be participating the study; each team involved a number of service users who will be affected by the intervention. A stepped wedge design was considered appropriate for the study due to logistical constraint in simultaneously introducing the intervention to more than a few teams. The experience-based co-design intervention will be delivered to the clusters in three waves, each with a duration of 9 months. Four teams will start the intervention in wave 1 and wave 2, respectively, while the remaining three teams receive the intervention in the final wave. In other words, the study includes four periods, with a baseline period lasting about 6 months.

The outcome of interest is the improvement in psychosocial recovery measured by the recovery assessment scale revised,⁴⁵ and was measured for each user at the end of baseline period and each of the three follow-up period. The standardized effect size on the psychosocial recovery outcome was estimated to be 0.35, and the within-period correlation was assumed to be τ = 0.1.⁴⁴ Since the study affords to randomize only 11 clusters, there may be a risk of inflated type I error rate with a z-test. As the t-test with DoF = I − 2 = 9 performs best with respect to empirical size and power in the simulations, we determine the required cohort size based on a test using expressions (12) and (13). Assuming the correlation decay is only moderate so that ρ = 0.8, power is estimated to be 0.79 for N = 8 and 0.81 for N = 9, barring drop out. Therefore, N = 9 is required to ensure 80% power given a 5% test size. We further conducted a sensitivity analysis to see how power changes according to the degree of correlation decay, and presented the power contour in panel (B) of Figure 2. Due to the small sample size and the heavy tail of the t distribution, the study is sensitive to correlation decay when τ = 0.1, and remains so even if τ approaches zero. On the other hand, the actual study planned to recruit N = 30 users in each team. With this larger cohort size, panel (C) of Figure 2 suggests that the power becomes less sensitive to the correlation decay, especially as the within-period correlation approaches zero. For example, if τ ≤ 0.02, the study power remains at least around 80% regardless of the amount of correlation decay.

6.3 |. Consequences of specifying a nondecay correlation structure

For GEE analyses of cohort stepped wedge designs, Li et al⁵ recently developed the block exchangeable correlation structure for estimating the sample size and power. With a continuous outcome and identity link function, the block exchangeable correlation structure is also implied by the linear mixed effects model for cohort studies discussed in the works of Hooper et al,⁸ and Girling and Hemming.⁴⁶ To understand the implications of alternative correlation models, we compare the variance of the intervention effect estimator obtained under the proportional decay correlation structure to that obtained under the block exchangeable correlation structure. Specifically, both the proportional decay structure and the block exchangeable correlation structure assume a constant within-period correlation, ie, corr(y_ijt, y_ij′t) = τ for j ≠ j′. However, the latter correlation model also assumes constant between-period and within-individual correlations such that $\text{corr}\left(y_{ijt},y_{i{j}^{\prime }{t}^{\prime }}\right)={\alpha }_1^{\text{BE}}$ for j ≠ j′, t ≠ t′ and $\text{corr}\left(y_{ijt},y_{ij{t}^{\prime }}\right)={\alpha }_2^{\text{BE}}$, t ≠ t′. To focus ideas, we carry out the comparisons based on the standard stepped wedge designs with a single baseline period and an equal number of clinics switching to intervention at each step. We refer to the variance of $\widehat{\delta }$ obtained under the proportional decay structure as ${\text{var}}^{\text{PD}}(\widehat{\delta })$, whose expression is provided in (16). We refer to the variance of $\widehat{\delta }$ obtained under the block exchangeable structure as ${\text{var}}^{\text{BE}}(\widehat{\delta })$, whose expression is derived in the work of Li et al⁵ as

with ${\lambda }_3=1+(N-1)\left(\tau -{\alpha }_1^{\text{BE}}\right)-{\alpha }_2^{\text{BE}}$ and ${\lambda }_4=1+(N-1)\tau +(T-1)(N-1){\alpha }_1^{\text{BE}}+(T-1){\alpha }_2^{\text{BE}}$ as the two distinct eigenvalues of the block exchangeable matrix. If we define function $h\left({\alpha }_1^{\text{BE}},{\alpha }_2^{\text{BE}}\right)=(N-1){\alpha }_1^{\text{BE}}+{\alpha }_2^{\text{BF}}$, we can write the relative variance as

For each value of the within-period correlation τ and each value of the decay parameter d = 1 − ρ in the proportional decay model, there may exist pairs of values for (${\alpha }_1^{\text{BE}}$, ${\alpha }_2^{\text{BE}}$) that result in the same variance of the intervention effect. Obtaining the equal variance correspondence is equivalent to finding the straight line $h\left({\alpha }_1^{\text{BE}},{\alpha }_2^{\text{BE}}\right)=\eta$ that solves ${\text{var}}^{\text{PD}}(\widehat{\delta })/{\text{var}}^{\text{BE}}(\widehat{\delta })=1$. Because the relative variance is quadratic in $h\left({\alpha }_1^{\text{BE}},{\alpha }_2^{\text{BE}}\right)$, such a line may not always exist within the plausible range of (${\alpha }_1^{\text{BE}}$, ${\alpha }_2^{\text{BE}}$) that ensures a positive definite block exchangeable correlation matrix. We confirm this observation by plotting the relative variance for varying correlation parameters. Fixing τ = 0.03, T = 4 and N = 20 similar to the accredited exercise physiologist study, we present in Figure 3 the contour of ${\text{var}}^{\text{PD}}(\widehat{\delta })/{\text{var}}^{\text{BE}}(\widehat{\delta })$ over the (${\alpha }_1^{\text{BE}}$, ${\alpha }_2^{\text{BE}}$) plane for each level of decay d ∈ {0, 1, …, 0.9}. The dashed thick line indicates the equal variance correspondence. We also present the contour plots by specifying N = 100, T = 8, and τ = 0.1 in Web Figures 5 to 7. It is evident that the existence and location of the equal variance line on the (${\alpha }_1^{\text{BE}}$, ${\alpha }_2^{\text{BE}}$) plane depends on the value of correlation decay, cohort size and number of periods. If the equal variance line exists, the value of τ and number of periods only affects its location while cohort size N further affects its slope, as expected from inspecting expression (18). Apart from the equal variance correspondence, the two variances are generally different. Depending on the gray shades (colored shades in the online version), either the proportional decay or the block exchangeable structure will lead to a larger variance of the intervention effect and require a larger sample size. As a result, using the block exchangeable correlation structure in the presence of correlation decay could lead to an overestimation or underestimation of the true intervention effect variance, and vice versa. Particularly, the variances returned from these two correlation models become close when ${\alpha }_1^{\text{BE}}$, ${\alpha }_2^{\text{BE}}$ approximate zero and when the decay parameter, d = 1 − ρ, approximates one. These two restrictions result in many nearly-zero entries in the block exchangeable and proportional decay correlation matrices, increasing the dependence of sample size estimation on the within-period correlation τ. Therefore, it is anticipated that the variances from the two models become similar in this particular scenario, even though in general there could be large differences between ${\text{var}}^{\text{PD}}(\widehat{\delta })$ and ${\text{var}}^{\text{BE}}(\widehat{\delta })$. To summarize, our key message for cohort stepped wedge designs echo the principal findings in the work of Kasza et al¹⁰ for cross-sectional designs. It is possible to grossly overestimate or underestimate the variance of the intervention effect if the correlation model is misspecified, except in restrictive scenarios. In practice, researchers could investigate the sensitivity of sample size estimates to misspecification of the correlation structure when there is limited preliminary data at the design stage of the trial.