SteppedPowerSteppedPowerSteppedPower offers tools for power and sample size calculation as well as design diagnostics for longitudinal mixed model settings, with a focus on stepped wedge designs. Other implemented study design types are parallel, parallel with baseline period(s) and crossover designs. Further design types can be easily defined by the user.
Currently, normal outcomes and binomial outcomes with logit link are implemented. The following random effects can be specified: random cluster intercept, random treatment effect, random subject specific intercept and random time effect. The covariance structure can be compound symmetry or autoregressive.
This code is modularised in order to be flexible and easy to use (and hopefully to maintain as well). At the same time, the use of sparse matrix classes of the matrix package makes computation of large designs feasible.
A common approach to model the correlation in longitudinal studies are random effects.(Hussey and Hughes 2007; Li et al. 2020) Such a model has the form
\[y_{ijk}= \mu_0 + \alpha_i + T_{ij} \theta + b_j + e_{ijk}\]
with
For power calculation, the standard deviation of random effects is assumed to be known. Lets define \(\beta:=(\mu,\alpha',\theta)'\) and \(\omega_{ijk}:=b_j+e_{ijk}\). This leads to a compact and more general notation of the above equation:
\[\begin{align} y_{ijk}&= X_{ij}\beta + \omega_{ijk}\\ \text{or, in matrix notation:} \qquad \\ y&=X\beta + \omega \end{align}\]
Where \(X\) is the corresponding design matrix and \(\omega\sim N(0,\Omega)\), where \(\Omega\) is a compound-symmetry (syn. exchangeable) variance matrix defined by \(\tau\) and \(\sigma\). We are thus in a weighted least squares setting, so the variance of \(\beta\) is
\[ \text{Var}(\hat\beta) = {(X'\Omega^{-1}X)^{-1}}\]
We can then calculate the power of a z-test
\[ \text{power} = \Phi\left(\frac{\theta_A-\theta_0}{\sqrt{\text{Var}(\hat \theta)}}- Z_{1-\frac{\alpha}{2}}\right) \]
where \(\text{Var}(\hat \theta)\) is the diagonal element of \(\Omega\) that corresponds to \(\hat\theta\).
Extensions to the above formula implemented in this package are
with leads to the following extended model formula:
\[y_{ijk}= g\big( \mu + \alpha_i + X (\theta_{ij} + c_j) + b_j + e_{ijk}\big)\] with
For most users, the probably most important function is wlsPower. It calls several auxiliary functions which will be shortly discussed here. This section is not essential for the usage of SteppedPower, it might be helpful to design non-standard user defined settings.
wlsPower is essentially just a flexible wrapper for the function compute_wlsPower, which does the actual computation.
compute_wlsPower then calls constuct_DesMat and construct_CovMat.
construct_DesMat builds the design matrix which consists of the treatment status, usually built byconstruct_trtMat and the time adjustment, usually built by construct_timeadjust. There is also the option to pass a user defined definition of the treatment status to construct_DesMat. If not specified, the number of timepoints is guessed as the fewest number of periods (timepoints) possible with the given design, i.e. two for cross-over designs or the number of waves plus one for stepped wedge designs.
construct_CovMat builds the covariance matrix (explicitly). It uses construct_CovBlk to construct the blocks for each cluster which are then combined to a block diagonal matrix.
In the weighted least squares setting, the estimator \(\hat \beta\) is a linear function of the data \(y\)
\[\hat\beta = \underbrace{(X'\Omega^{-1}X)^{-1}(X'\Omega^{-1})}_{=:\text{M}}\cdot y\]
with \(X\) the design matrix and \(\Omega\) the covariance matrix as above. The matrix \(M\) gives an impression of the importance of clusters and time periods with regard to the estimated coefficients \(\hat\beta\). The first row of \(M\) corresponds to the coefficient of the treatment status, i.e. the treatment effect.
The plot.wlsPower method visualises this first row of \(M\) as a matrix where rows and columns correspond to clusters and time periods, respectively.
Furthermore, to give a rough comparison of importance between clusters (or between time periods), the sum of absolute weights per row (or per column) is also shown.
CAVE: These are out-of-bag estimates for the influence of observations on \(\hat\theta\), but not for \(\text{Var}(\hat\theta)\) !
When the argument Power is passed to wlsPower, the sample size needed is calculated, under the assumption of equally sized clusters and periods.
This might be a proof of concept rather than an example with practical relevance, but let’s to compare the mean in two groups. For two groups of 10 observations each, the power of a Z-test can be calculated as follows:
wlsPower(Cl=c(10,10), mu0=0,mu1=.6,sigma=1, tau=0, N=1,
dsntype="parallel", timepoints=1)
#> Power = 0.2687
#> Significance level (two sided) = 0.05
## the same:
wlsPower(Cl=c(1,1), mu0=0,mu1=.6, sigma=1, tau=0, N=10,
dsntype="parallel", timepoints=1)
#> Power = 0.2687
#> Significance level (two sided) = 0.05A quick Note on t-tests: It is much more challenging to use
SteppedPowerto reproduce settings in which the variance is assumed to be unknown, most prominently the well known t-test. In this package, you find implemented some (experimental) heuristics for guessing the denominator degrees of freedom, but they yield rather scaled Wald tests than t tests. The main difference is that the distribution under the alternative is assumed to be symmetric, whereas the t-test assumes a non-central (hence skewed) t-distribution.
Periods in which no cluster switches to the intervention are specified by inserting zeros into the Cl argument, i.e. Cl=c(4,4,4,0).
mod1 <- wlsPower(Cl=c(1,1,1,0), mu0=0, mu1=1,
sigma=0.4, tau=0, verbose=2)
knitr::kable(mod1$DesignMatrix$trtMat)| 0 | 1 | 1 | 1 | 1 |
| 0 | 0 | 1 | 1 | 1 |
| 0 | 0 | 0 | 1 | 1 |
The argument N defines the cluster size. N can be * a scalar, if all clusters have the same assumed size, which is also constant over time * a vector, if the size differs between clusters but is assumed to be constant over time * a matrix where each row corresponds to either a cluster or a wave of clusters and each column corresponds to a timepoint
Suppose you do not plan to observe all clusters over the whole study period. Rather, clusters that switch early to the intervention are not observed until the end. Analogous, observation starts later in clusters that switch towards the end of the study. This is sometimes called ‘incomplete SWD’ [hemming2015stepped].
There are two ways to achieve this in SteppedPower, both by using the incomplete argument. One can either scalar, which then defines the number of observed periods before and after the switch from control to intervention in each cluster.
If for example the study consists of eight clusters in four sequences (i.e. five timepoints), and we observe two timepoints before and after the switch, then we receive
incompletePwr <- wlsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0,mu1=.5, N=80,
incomplete=2, verbose=2)
incompletePwr
#> Power = 0.8221
#> Significance level (two sided) = 0.05A slightly more tedious, but more flexible way is to define a matrix where each row corresponds to either a cluster or a wave of clusters and each column corresponds to a timepoint. If a cluster is not observed at a specific timepoint, set the value in the corresponding cell to 0. For the example above, such a matrix would look like this:
TM <- toeplitz(c(1,1,0,0))
incompleteMat1 <- cbind(TM[,1:2],rep(1,4),TM[,3:4])
incompleteMat2 <- incompleteMat1[rep(1:4,each=2),]A matrix where each row represents a wave of clusters
| 1 | 1 | 1 | 0 | 0 |
| 1 | 1 | 1 | 1 | 0 |
| 0 | 1 | 1 | 1 | 1 |
| 0 | 0 | 1 | 1 | 1 |
or each row represents a cluster
| 1 | 1 | 1 | 0 | 0 |
| 1 | 1 | 1 | 0 | 0 |
| 1 | 1 | 1 | 1 | 0 |
| 1 | 1 | 1 | 1 | 0 |
| 0 | 1 | 1 | 1 | 1 |
| 0 | 1 | 1 | 1 | 1 |
| 0 | 0 | 1 | 1 | 1 |
| 0 | 0 | 1 | 1 | 1 |
Now all that’s left to do is to plug that into the main function:
incompletePwr1 <- wlsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0, mu1=.5, N=80,
incomplete=incompleteMat1, verbose=2)
incompletePwr2 <- wlsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0, mu1=.5, N=80,
incomplete=incompleteMat2, verbose=2)
all.equal(incompletePwr,incompletePwr1)
#> [1] TRUE
all.equal(incompletePwr,incompletePwr2)
#> [1] TRUEWe can also have a quick look at the projection matrix where we see that the clusters have a weight of exactly zero at the timepoints where they are not observed
The argument
incompletewith matrix input works also for other design types, but makes (supposedly) most sense in the context of stepped wedge designs
The most usual method for the modelling of potential secular trends is to take time period as a factor into the analysis model (Hussey and Hughes 2007; Hemming et al. 2015; Hemming and Taljaard 2020).
For diagnostic purposes (or for bold users), some other adjustment options are implemented.
TimeAdj1 <- wlsPower(Cl=rep(2,4), mu0=0, mu1=1, tau=0,
timeAdjust="linear", verbose=2)
TimeAdj2 <- wlsPower(Cl=rep(2,4), mu0=0, mu1=1, tau=0,
timeAdjust="factor", verbose=2)Design matrix of the first cluster with linear adjustment for secular trend:
| trt | ||
|---|---|---|
| 0 | 1 | 0.2 |
| 1 | 1 | 0.4 |
| 1 | 1 | 0.6 |
| 1 | 1 | 0.8 |
| 1 | 1 | 1.0 |
Design matrix of the first cluster with categorical adjustment for secular trend:
| trt | |||||
|---|---|---|---|---|---|
| 0 | 1 | 0 | 0 | 0 | 0 |
| 1 | 1 | 1 | 0 | 0 | 0 |
| 1 | 1 | 0 | 1 | 0 | 0 |
| 1 | 1 | 0 | 0 | 1 | 0 |
| 1 | 1 | 0 | 0 | 0 | 1 |
In a closed cohort the patients are observed over the whole study period. The same correlation structure arises in cross sectional stepped wedge designs if subclusters exist (such as wards within clinics). The argument psi denotes the standard deviation of a random subject (or subcluster) specific intercept.
The power is calculated on aggregated cluster means:
Closed1 <- wlsPower(mu0=0, mu1=5, Cl=rep(3,3), sigma=5, tau=1, psi=2, gamma=1,
N=3, verbose=2)
a <- plot(Closed1$DesignMatrix)WithINIDV_LVL=TRUE, the calculation is done on the individual level. This yields the same results but is far more comutationally expensive and is mainly intended for diagnostic purposes.
sessionInfo()
#> R version 4.0.3 (2020-10-10)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19041)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=C LC_CTYPE=German_Germany.1252
#> [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
#> [5] LC_TIME=German_Germany.1252
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] Matrix_1.2-18 SteppedPower_0.1.0 knitr_1.30
#>
#> loaded via a namespace (and not attached):
#> [1] highr_0.8 pillar_1.4.7 compiler_4.0.3 tools_4.0.3
#> [5] digest_0.6.27 viridisLite_0.3.0 jsonlite_1.7.2 evaluate_0.14
#> [9] lifecycle_0.2.0 tibble_3.0.4 gtable_0.3.0 lattice_0.20-41
#> [13] pkgconfig_2.0.3 rlang_0.4.10 crosstalk_1.1.1 yaml_2.2.1
#> [17] xfun_0.20 pwr_1.3-0 stringr_1.4.0 dplyr_1.0.2
#> [21] httr_1.4.2 generics_0.1.0 vctrs_0.3.6 htmlwidgets_1.5.3
#> [25] grid_4.0.3 tidyselect_1.1.0 glue_1.4.2 data.table_1.13.6
#> [29] R6_2.5.0 plotly_4.9.3 rmarkdown_2.6 farver_2.0.3
#> [33] tidyr_1.1.2 ggplot2_3.3.3 purrr_0.3.4 magrittr_2.0.1
#> [37] scales_1.1.1 ellipsis_0.3.1 htmltools_0.5.1.1 colorspace_2.0-0
#> [41] stringi_1.5.3 lazyeval_0.2.2 munsell_0.5.0 crayon_1.3.4Hemming, Karla, Terry P Haines, Peter J Chilton, Alan J Girling, and Richard J Lilford. 2015. “The Stepped Wedge Cluster Randomised Trial: Rationale, Design, Analysis, and Reporting.” Bmj 350: h391.
Hemming, Karla, and Monica Taljaard. 2020. “Reflection on Modern Methods: When Is a Stepped-Wedge Cluster Randomized Trial a Good Study Design Choice?” International Journal of Epidemiology.
Hussey, Michael A, and James P Hughes. 2007. “Design and Analysis of Stepped Wedge Cluster Randomized Trials.” Contemporary Clinical Trials 28 (2): 182–91.
Li, Fan, James P Hughes, Karla Hemming, Monica Taljaard, Edward R Melnick, and Patrick J Heagerty. 2020. “Mixed-Effects Models for the Design and Analysis of Stepped Wedge Cluster Randomized Trials: An Overview.” Statistical Methods in Medical Research, 0962280220932962.