Nothing
Getting Started with `SteppedPower`"
In SteppedPower: Power Calculation for Stepped Wedge Designs
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 7,
warning = FALSE
)
library(knitr)
library(SteppedPower)
library(Matrix)
library(plotly)
About SteppedPower
SteppedPower 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
identity 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 package 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.
Methods and Notation
A common approach to model the correlation in
longitudinal studies are random effects [@hussey2007design; @li2020mixed].
Such a model has the form
$$y_{ijk}= T_{ij} \theta + c_i + \mu_j + e_{ijk}$$
with
- $y_{ijk}$ the response in cluster $i$ at time $j$ for individual $k$
- $\mu_j$ the time trend at time $j$
- $T_{ij}$ indicates the treatment status (0 = control, 1 = interventional treatment)
- $\theta$ the treatment effect
- $c_i$ a random cluster effect for cluster $i$ with $c_i \sim N(0,\tau^2)$
- $e_{ijk}$ a normal random error term with $e_{ijk}\sim N(0,\sigma^2)$
For power calculation, the standard deviation of random effects is assumed to be known.
Let's define $\beta:=(\theta,\mu')'$ and $\omega_{ijk}:=c_i+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}= T (\theta_{ij} + d_i) + c_i + \mu_j + t_{ij} + s_{ijk} + e_{ijk}$$
with
- $d_i$ a random treatment effect
- $c_i$ and $d_i$ jointly distributed with $\left(\begin{smallmatrix} \tau^2 & \rho\tau\eta \ \rho\tau\eta & \eta^2\end{smallmatrix}\right)$
- $s_{ijk}$ a random subject specific effect
- $t_{ij}$ a random time effect
Specification of the correlation structure
A different approach to establish a covariance structure is to directly assume cluster cell correlation as proposed by \cite{Preisser2003integrated} and introduced to the context of stepped wedge designs in [@li2018sample]. They use $\boldsymbol{\alpha}:=(\alpha_0, \alpha_1, \alpha_2)$ where
$\alpha_0:=\operatorname{corr}({y_{ijk},y_{ijk'}})$ for $ k\neq k'$ defines the within-period correlation; $\alpha_1:=\operatorname{corr}({y_{ijk},y_{ij'k'}})$ for $j\neq j', \, k\neq k'$ defines the inter-period correlation and $\alpha_0:=\operatorname{corr}({y_{ijk},y_{ij'k}})$ for $j\neq j'$ defines the within-individual correlation.
One receives an easily interpretable correlation matrix, albeit at the cost of some flexibility. We can translate $\boldsymbol{\alpha}$ into a representation with random effects, given a marginal variance $\sigma^2_{\text{marg}}$.
$$
\begin{align}
\tau^2 &:= {\sigma^2_{\text{marg}} \alpha_1} \
\gamma^2 &:= {\sigma^2_{\text{marg}} (\alpha_0-\alpha_1)} \
\psi^2 &:= {\sigma^2_{\text{marg}} (\alpha_2-\alpha_2)} \
\sigma^2_{\text{res}} &:= {\sigma^2_{\text{marg}} (1-\alpha_0 - \alpha_2 + \alpha_1)}
\end{align}
$$
Note that $\boldsymbol{\alpha}$ as defined above cannot specify a random treatment effect.
A quick tour
For most users, the probably most important function is glsPower.
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.
glsPower is essentially just a flexible wrapper for the function
compute_glsPower, which does the actual computation.
compute_glsPower then calls construct_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.
Features
Plot Method
The plot.glsPower method makes up to four plots, visualising the influence
of cluster-period cells, the information content of cluster-period cells,
the treatment matrix $T$ and the covariance matrix
(selectable by which=).
By default, only the first plot is returned. The information content is also
plotted if it the object to be plotted contains the corresponding matrix
(which is usually the case). It is constructed by exploiting the
approximation to a weighted least squares setting. Therefore,
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 $H$ gives an impression of the importance of clusters and
time periods with regard to the estimated coefficients $\hat\beta$.
The first row of $H$ corresponds to the coefficient of the treatment status,
i.e. the treatment effect.
The plot.glsPower 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.
glsPwr <- glsPower(Cl=c(3,2,3), mu0=0, mu1=1, sigma=1, tau=.5, verbose=2)
plot(glsPwr,which=1, show_colorbar=FALSE)$WgtPlot
The information content, as described in [@kasza2019information], looks like this:
plot(glsPwr,which=2, show_colorbar=FALSE)$ICplot
CAVE: The influence plot (above) yields out-of-bag estimates for the influence of
observations on $\hat\theta$, but not for $\text{Var}(\hat\theta)$.
The information content (below) visualises the relative change in
$\text{Var}(\hat\theta)$ when one particular cluster period is omitted.
Find Sample Size for given Power
When the argument power is passed to glsPower,
the sample size needed is calculated, under the assumption of equally sized
clusters and periods.
glsPower(Cl=c(3,3,3), mu0=0, mu1=.2, sigma=1, tau=0, power=.8)
So in this setting, you would need $50$ individuals per period in each cluster to
achieve a target power of $80\%$
Use cases and examples
Comparison of two groups -- Z-Test
This might be a proof of concept rather than an example with
practical relevance, but let's try 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:
glsPower(Cl=c(10,10), mu0=0,mu1=1.2,sigma=1, tau=0, N=1,
dsntype="parallel", timepoints=1)$power
## the same:
glsPower(Cl=c(1,1), mu0=0,mu1=1.2, sigma=1, tau=0, N=10,
dsntype="parallel", timepoints=1)$power
pwr::pwr.norm.test(.6,n=20)$power
A quick note on t-tests: It is much more challenging to use SteppedPower
to 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.
Longitudinal study -- parallel groups
glsPower(Cl=c(10,10),timepoints=5,mu0=0,mu1=.25,
sigma=.5,dsntype="parallel")
glsPower(Cl=c(10,10),timepoints=5,mu0=0,mu1=.25,
sigma=.5,tau=.2,dsntype="parallel")
Stepped Wedge designs with empty sequences (i.e. waves)
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 <- glsPower(Cl=c(1,1,1,0), mu0=0, mu1=1,
sigma=0.4, tau=0, verbose=2)
The treatment matrix is then stored under mod1$DesignMatrix$trtMat:
knitr::kable(mod1$DesignMatrix$trtMat)
Autocorrelated random effects
In longitudinal studies, it can be sensible to assume that correlation within
clusters decreases with increasing time lag.
The argument AR enables the user to specify a AR-1 correlation. AR
must be any value between 0 and 1.
The former corresponds to i.i.d. observations,
the latter to the usual compound symmetry covariance type.
Optionally, you can pass a vector of up to length three to AR. Then, the first element
applies to the cluster intercept, the second to the treatment effect (i.e. slope)
and the third to the subject specific intercept. This is especially useful for
open cohort designs (see below).
An example of a stepped wedge design with 8 clusters in 4 waves,
once with medium autocorrelation (AR=0.6) and
once with high autocorrelation (AR=0.95):
mod2 <- glsPower(Cl=c(2,2,2,2), mu0=0, mu1=1,
sigma=1, N=100, tau=1, AR=.6, verbose=2)
mod3 <- glsPower(Cl=c(2,2,2,2), mu0=0, mu1=1,
sigma=1, N=100, tau=1, AR=.95, verbose=2)
For AR=0.6, the covariance matrix within one cluster then looks like this:
suppressWarnings(knitr::kable(as.matrix(mod2$CovarianceMatrix[1:5,1:5])))
For AR=0.95 it takes the following shape
suppressWarnings(knitr::kable(as.matrix(mod3$CovarianceMatrix[1:5,1:5])))
rm(mod1,mod2,mod3)
Unequal cluster sizes
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
mod4 <- glsPower(Cl=c(1,1,1), mu0=0, mu1=1, N=c(1,3,10),
sigma=1, tau=.5, verbose=2)
plot(mod4, which=2, show_colorbar=FALSE)$ICplot
rm(mod4)
Incomplete Stepped Wedge Designs
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 <- glsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0,mu1=.5, N=80,
incomplete=2, verbose=2)
incompletePwr
A 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
suppressWarnings(knitr::kable(incompleteMat1))
or each row represents a cluster
suppressWarnings(knitr::kable(incompleteMat2))
Now all that's left to do is to plug that into the main function:
incompletePwr1 <- glsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0, mu1=.5, N=80,
incomplete=incompleteMat1, verbose=2)
incompletePwr2 <- glsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0, mu1=.5, N=80,
incomplete=incompleteMat2, verbose=2)
all.equal(incompletePwr,incompletePwr1)
all.equal(incompletePwr,incompletePwr2)
We 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
plot(incompletePwr, show_colorbar=FALSE)$WgtPlot
rm(incompletePwr,incompletePwr1,incompletePwr2,incompleteMat1,incompleteMat2)
The argument incomplete with matrix input works also for other design types,
but makes (supposedly) most sense in the context of stepped wedge designs
Adjustment for Secular Trends
The most usual method for the modelling of potential secular trends is
to take time period as a factor into the analysis model
[@hussey2007design; @hemming2015stepped; @hemming2020reflection].
For diagnostic purposes (or for bold users), some other adjustment options are
implemented.
TimeAdj1 <- glsPower(Cl=rep(2,4), mu0=0, mu1=1, sigma=1, tau=0,
timeAdjust="linear", verbose=2)
TimeAdj2 <- glsPower(Cl=rep(2,4), mu0=0, mu1=1, sigma=1, tau=0,
timeAdjust="factor", verbose=2)
Design matrix of the first cluster with linear adjustment for secular trend:
knitr::kable(head(TimeAdj1$DesignMatrix$dsnmatrix, 5))
Design matrix of the first cluster with categorical adjustment for secular trend:
knitr::kable(head(TimeAdj2$DesignMatrix$dsnmatrix, 5))
Closed cohort SWD
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 <- glsPower(mu0=0, mu1=5, Cl=rep(3,3),
sigma=5, tau=1, psi=3,
N=3, verbose=2)
Closed1
With INDIV_LVL=TRUE, the calculation is done on the individual level.
This yields the same results but is far more computationally expensive and is
mainly intended for diagnostic purposes.
Closed2 <- glsPower(mu0=0, mu1=5, Cl=rep(3,3),
sigma=5, tau=1, psi=3,
N=3, verbose=2, INDIV_LVL = TRUE)
Closed2
plot(Closed2, annotations=FALSE, show_colorbar=FALSE)$WgtPlot
Closed1$power - Closed2$power
Open cohort designs
The AR argument mentioned above specifies the autocorrelation of random effects.
It offers a convenient way to model open cohort designs. The third element of AR
can be seen as the estimated probability that a subject included in period $j$
is also included in period $j+1$, instead of being replaced by a new subject.
Expanding on the closed cohort example above, one could assume that subjects
have a 75% chance to reappear in the next study period.
Open1 <- glsPower(mu0=0, mu1=5, Cl=rep(3,3),
sigma=5, tau=1, psi=3, AR=c(1,1,.75), N=3)
Closed1$power
Open1$power
So this would result in a power loss of roughly
r round( Closed1$power - Open1$power,3)*100%
compared to a closed cohort design in this case.
Covariance structure for open cohort designs
This proposed covariance structure is then a mixture of exchangeable (on cluster level)
and autoregressive (on subject level). Let's have a look at a toy example with
three clusters, each consisting of three subjects. The probability that an
individual reappears in two consecutive periods is estimated to be 60%.
The covariance matrix (on subject level) could be visualised as follows:
Open2Indiv <- glsPower(mu0=0, mu1=10, Cl=c(1,1,1,0),
sigma=1, tau=5, psi=10, AR=c(1,1,.60),
N=3, verbose=2, INDIV_LVL=TRUE)
plot(Open2Indiv, which=4, show_colorbar=FALSE)$CMplot
Note that you can zoom to get a better look at the structure. The autoregressive
structure on the subject level in the small boxes along the main diagonal and the
block exchangeable structure on cluster level becomes evident.
Reproducing common closed formulae
There exist some formulae that do not explicitly construct design and
covariance matrices, most notably [@kasza2020sample] and [@li2020design].
Both include the one presented in [@hussey2007design] as a special case.
These formulae are also implemented in SteppedPower.
As a showcase, here is an example that is loosely based
on the example in [@hussey2007design]:
trtMat <- construct_DesMat(c(6,6,6,6))$trtMat
mu0 <- 0.05 ; mu1 <- 0.032 ; N <- 100
tau <- .025 ; sigma <- sqrt(.041*.959)
gamma <- 0.01 ; psi <- .1 ; chi <- .5 ; AR <- .5
Closed formula for proportional decay structure
VarClosed_Li computes the treatment variance for stepped wedge designs with a
autoregressive (syn. proportional decay) structure. It can be used to compute
a variance for the treatment effect and thus the power for the setting in question:
tmp <- VarClosed_Li(trtMat, tau=tau, psi=psi, N=N, AR=AR)
tTestPwr(mu0-mu1, se=sqrt(tmp), df=Inf)
Explicit calculation yields the same:
a <- SteppedPower::glsPower(Cl=rep(6,4), mu0=mu0, mu1=mu1, AR=AR,
sigma=0, tau=tau, N=N, psi=psi, verbose=1, INDIV_LVL = TRUE)
a
Closed formula for open cohort designs
Kasza et al. offer an extension to open cohort designs by the use of an attrition
or churn rate chi [@kasza2020sample]. It can be understood as the estimated probability
that a subject included in period $j$
is also included in period $j'$ for $j\ne j'$. Note that this differs from the
approach mentioned above, as it results in a mixture of exchangeable covariance matrices.
Since it differs from the way SteppedPower intents to model open cohort designs,
it requires a bit of tinkering to reproduce the results of the closed formula with
this package. It can be done with the help of the random time effect (gamma).
tmp <- VarClosed_Kasza(trtMat, sigma=sigma, tau=tau, gamma=gamma, psi=psi, N=N, chi=0)
tTestPwr(mu0-mu1, se=sqrt(tmp), df = Inf)
glsPower(Cl = rep(6,4), N=N, mu0=mu0, mu1=mu1, verbose=0,
sigma=sigma, tau=tau, gamma=gamma, psi=psi)
tmp <- VarClosed_Kasza(trtMat, sigma=sigma, tau=tau, gamma=gamma, psi=psi, N=N, chi=1)
tTestPwr(mu0-mu1, sqrt(tmp), df = Inf)
glsPower(Cl = rep(6,4), N=N, mu0=mu0, mu1=mu1, verbose=0,
sigma=sigma, tau=tau,
gamma=sqrt(gamma^2+psi^2/N), psi=0)
tmp <- VarClosed_Kasza(trtMat, sigma=sigma, tau=tau, gamma=gamma, psi=psi, N=N, chi=chi)
tTestPwr(mu0-mu1, sqrt(tmp), df = Inf)
glsPower(Cl = rep(6,4), N=N, mu0=mu0, mu1=mu1, verbose=0,
sigma=sigma, tau=tau,
gamma=sqrt(gamma^2+chi*psi^2/N), psi=sqrt(1-chi)*psi)
Session Info {-}
print(sessionInfo(),locale=FALSE)
References {-}
Try the SteppedPower package in your browser
Any scripts or data that you put into this service are public.
SteppedPower documentation built on Sept. 13, 2023, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, warning = FALSE ) library(knitr) library(SteppedPower) library(Matrix) library(plotly)
About SteppedPower
SteppedPower 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 identity 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 package 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.
Methods and Notation
A common approach to model the correlation in longitudinal studies are random effects [@hussey2007design; @li2020mixed]. Such a model has the form
$$y_{ijk}= T_{ij} \theta + c_i + \mu_j + e_{ijk}$$
with
- $y_{ijk}$ the response in cluster $i$ at time $j$ for individual $k$
- $\mu_j$ the time trend at time $j$
- $T_{ij}$ indicates the treatment status (0 = control, 1 = interventional treatment)
- $\theta$ the treatment effect
- $c_i$ a random cluster effect for cluster $i$ with $c_i \sim N(0,\tau^2)$
- $e_{ijk}$ a normal random error term with $e_{ijk}\sim N(0,\sigma^2)$
For power calculation, the standard deviation of random effects is assumed to be known. Let's define $\beta:=(\theta,\mu')'$ and $\omega_{ijk}:=c_i+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}= T (\theta_{ij} + d_i) + c_i + \mu_j + t_{ij} + s_{ijk} + e_{ijk}$$ with
- $d_i$ a random treatment effect
- $c_i$ and $d_i$ jointly distributed with $\left(\begin{smallmatrix} \tau^2 & \rho\tau\eta \ \rho\tau\eta & \eta^2\end{smallmatrix}\right)$
- $s_{ijk}$ a random subject specific effect
- $t_{ij}$ a random time effect
Specification of the correlation structure
A different approach to establish a covariance structure is to directly assume cluster cell correlation as proposed by \cite{Preisser2003integrated} and introduced to the context of stepped wedge designs in [@li2018sample]. They use $\boldsymbol{\alpha}:=(\alpha_0, \alpha_1, \alpha_2)$ where $\alpha_0:=\operatorname{corr}({y_{ijk},y_{ijk'}})$ for $ k\neq k'$ defines the within-period correlation; $\alpha_1:=\operatorname{corr}({y_{ijk},y_{ij'k'}})$ for $j\neq j', \, k\neq k'$ defines the inter-period correlation and $\alpha_0:=\operatorname{corr}({y_{ijk},y_{ij'k}})$ for $j\neq j'$ defines the within-individual correlation.
One receives an easily interpretable correlation matrix, albeit at the cost of some flexibility. We can translate $\boldsymbol{\alpha}$ into a representation with random effects, given a marginal variance $\sigma^2_{\text{marg}}$.
$$ \begin{align} \tau^2 &:= {\sigma^2_{\text{marg}} \alpha_1} \ \gamma^2 &:= {\sigma^2_{\text{marg}} (\alpha_0-\alpha_1)} \ \psi^2 &:= {\sigma^2_{\text{marg}} (\alpha_2-\alpha_2)} \ \sigma^2_{\text{res}} &:= {\sigma^2_{\text{marg}} (1-\alpha_0 - \alpha_2 + \alpha_1)} \end{align} $$
Note that $\boldsymbol{\alpha}$ as defined above cannot specify a random treatment effect.
A quick tour
For most users, the probably most important function is glsPower.
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.
glsPower is essentially just a flexible wrapper for the function
compute_glsPower, which does the actual computation.
compute_glsPower then calls construct_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.
Features
Plot Method
The plot.glsPower method makes up to four plots, visualising the influence
of cluster-period cells, the information content of cluster-period cells,
the treatment matrix $T$ and the covariance matrix
(selectable by which=).
By default, only the first plot is returned. The information content is also plotted if it the object to be plotted contains the corresponding matrix (which is usually the case). It is constructed by exploiting the approximation to a weighted least squares setting. Therefore, 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 $H$ gives an impression of the importance of clusters and
time periods with regard to the estimated coefficients $\hat\beta$.
The first row of $H$ corresponds to the coefficient of the treatment status,
i.e. the treatment effect.
The plot.glsPower 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.
glsPwr <- glsPower(Cl=c(3,2,3), mu0=0, mu1=1, sigma=1, tau=.5, verbose=2) plot(glsPwr,which=1, show_colorbar=FALSE)$WgtPlot
The information content, as described in [@kasza2019information], looks like this:
plot(glsPwr,which=2, show_colorbar=FALSE)$ICplot
CAVE: The influence plot (above) yields out-of-bag estimates for the influence of observations on $\hat\theta$, but not for $\text{Var}(\hat\theta)$. The information content (below) visualises the relative change in $\text{Var}(\hat\theta)$ when one particular cluster period is omitted.
Find Sample Size for given Power
When the argument power is passed to glsPower,
the sample size needed is calculated, under the assumption of equally sized
clusters and periods.
glsPower(Cl=c(3,3,3), mu0=0, mu1=.2, sigma=1, tau=0, power=.8)
So in this setting, you would need $50$ individuals per period in each cluster to achieve a target power of $80\%$
Use cases and examples
Comparison of two groups -- Z-Test
This might be a proof of concept rather than an example with practical relevance, but let's try 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:
glsPower(Cl=c(10,10), mu0=0,mu1=1.2,sigma=1, tau=0, N=1, dsntype="parallel", timepoints=1)$power ## the same: glsPower(Cl=c(1,1), mu0=0,mu1=1.2, sigma=1, tau=0, N=10, dsntype="parallel", timepoints=1)$power
pwr::pwr.norm.test(.6,n=20)$power
A 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.
Longitudinal study -- parallel groups
glsPower(Cl=c(10,10),timepoints=5,mu0=0,mu1=.25, sigma=.5,dsntype="parallel") glsPower(Cl=c(10,10),timepoints=5,mu0=0,mu1=.25, sigma=.5,tau=.2,dsntype="parallel")
Stepped Wedge designs with empty sequences (i.e. waves)
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 <- glsPower(Cl=c(1,1,1,0), mu0=0, mu1=1, sigma=0.4, tau=0, verbose=2)
The treatment matrix is then stored under mod1$DesignMatrix$trtMat:
knitr::kable(mod1$DesignMatrix$trtMat)
Autocorrelated random effects
In longitudinal studies, it can be sensible to assume that correlation within
clusters decreases with increasing time lag.
The argument AR enables the user to specify a AR-1 correlation. AR
must be any value between 0 and 1.
The former corresponds to i.i.d. observations,
the latter to the usual compound symmetry covariance type.
Optionally, you can pass a vector of up to length three to AR. Then, the first element
applies to the cluster intercept, the second to the treatment effect (i.e. slope)
and the third to the subject specific intercept. This is especially useful for
open cohort designs (see below).
An example of a stepped wedge design with 8 clusters in 4 waves,
once with medium autocorrelation (AR=0.6) and
once with high autocorrelation (AR=0.95):
mod2 <- glsPower(Cl=c(2,2,2,2), mu0=0, mu1=1, sigma=1, N=100, tau=1, AR=.6, verbose=2) mod3 <- glsPower(Cl=c(2,2,2,2), mu0=0, mu1=1, sigma=1, N=100, tau=1, AR=.95, verbose=2)
For AR=0.6, the covariance matrix within one cluster then looks like this:
suppressWarnings(knitr::kable(as.matrix(mod2$CovarianceMatrix[1:5,1:5])))
For AR=0.95 it takes the following shape
suppressWarnings(knitr::kable(as.matrix(mod3$CovarianceMatrix[1:5,1:5]))) rm(mod1,mod2,mod3)
Unequal cluster sizes
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
mod4 <- glsPower(Cl=c(1,1,1), mu0=0, mu1=1, N=c(1,3,10), sigma=1, tau=.5, verbose=2) plot(mod4, which=2, show_colorbar=FALSE)$ICplot
rm(mod4)
Incomplete Stepped Wedge Designs
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 <- glsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0,mu1=.5, N=80, incomplete=2, verbose=2) incompletePwr
A 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
suppressWarnings(knitr::kable(incompleteMat1))
or each row represents a cluster
suppressWarnings(knitr::kable(incompleteMat2))
Now all that's left to do is to plug that into the main function:
incompletePwr1 <- glsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0, mu1=.5, N=80, incomplete=incompleteMat1, verbose=2) incompletePwr2 <- glsPower(Cl=rep(2,4), sigma=2, tau=.6, mu0=0, mu1=.5, N=80, incomplete=incompleteMat2, verbose=2) all.equal(incompletePwr,incompletePwr1) all.equal(incompletePwr,incompletePwr2)
We 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
plot(incompletePwr, show_colorbar=FALSE)$WgtPlot
rm(incompletePwr,incompletePwr1,incompletePwr2,incompleteMat1,incompleteMat2)
The argument
incompletewith matrix input works also for other design types, but makes (supposedly) most sense in the context of stepped wedge designs
Adjustment for Secular Trends
The most usual method for the modelling of potential secular trends is to take time period as a factor into the analysis model [@hussey2007design; @hemming2015stepped; @hemming2020reflection].
For diagnostic purposes (or for bold users), some other adjustment options are implemented.
TimeAdj1 <- glsPower(Cl=rep(2,4), mu0=0, mu1=1, sigma=1, tau=0, timeAdjust="linear", verbose=2) TimeAdj2 <- glsPower(Cl=rep(2,4), mu0=0, mu1=1, sigma=1, tau=0, timeAdjust="factor", verbose=2)
Design matrix of the first cluster with linear adjustment for secular trend:
knitr::kable(head(TimeAdj1$DesignMatrix$dsnmatrix, 5))
Design matrix of the first cluster with categorical adjustment for secular trend:
knitr::kable(head(TimeAdj2$DesignMatrix$dsnmatrix, 5))
Closed cohort SWD
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 <- glsPower(mu0=0, mu1=5, Cl=rep(3,3), sigma=5, tau=1, psi=3, N=3, verbose=2) Closed1
With INDIV_LVL=TRUE, the calculation is done on the individual level.
This yields the same results but is far more computationally expensive and is
mainly intended for diagnostic purposes.
Closed2 <- glsPower(mu0=0, mu1=5, Cl=rep(3,3), sigma=5, tau=1, psi=3, N=3, verbose=2, INDIV_LVL = TRUE) Closed2 plot(Closed2, annotations=FALSE, show_colorbar=FALSE)$WgtPlot
Closed1$power - Closed2$power
Open cohort designs
The AR argument mentioned above specifies the autocorrelation of random effects.
It offers a convenient way to model open cohort designs. The third element of AR
can be seen as the estimated probability that a subject included in period $j$
is also included in period $j+1$, instead of being replaced by a new subject.
Expanding on the closed cohort example above, one could assume that subjects have a 75% chance to reappear in the next study period.
Open1 <- glsPower(mu0=0, mu1=5, Cl=rep(3,3), sigma=5, tau=1, psi=3, AR=c(1,1,.75), N=3) Closed1$power Open1$power
So this would result in a power loss of roughly
r round( Closed1$power - Open1$power,3)*100%
compared to a closed cohort design in this case.
Covariance structure for open cohort designs
This proposed covariance structure is then a mixture of exchangeable (on cluster level) and autoregressive (on subject level). Let's have a look at a toy example with three clusters, each consisting of three subjects. The probability that an individual reappears in two consecutive periods is estimated to be 60%. The covariance matrix (on subject level) could be visualised as follows:
Open2Indiv <- glsPower(mu0=0, mu1=10, Cl=c(1,1,1,0), sigma=1, tau=5, psi=10, AR=c(1,1,.60), N=3, verbose=2, INDIV_LVL=TRUE) plot(Open2Indiv, which=4, show_colorbar=FALSE)$CMplot
Note that you can zoom to get a better look at the structure. The autoregressive structure on the subject level in the small boxes along the main diagonal and the block exchangeable structure on cluster level becomes evident.
Reproducing common closed formulae
There exist some formulae that do not explicitly construct design and covariance matrices, most notably [@kasza2020sample] and [@li2020design]. Both include the one presented in [@hussey2007design] as a special case.
These formulae are also implemented in SteppedPower.
As a showcase, here is an example that is loosely based
on the example in [@hussey2007design]:
trtMat <- construct_DesMat(c(6,6,6,6))$trtMat mu0 <- 0.05 ; mu1 <- 0.032 ; N <- 100 tau <- .025 ; sigma <- sqrt(.041*.959) gamma <- 0.01 ; psi <- .1 ; chi <- .5 ; AR <- .5
Closed formula for proportional decay structure
VarClosed_Li computes the treatment variance for stepped wedge designs with a
autoregressive (syn. proportional decay) structure. It can be used to compute
a variance for the treatment effect and thus the power for the setting in question:
tmp <- VarClosed_Li(trtMat, tau=tau, psi=psi, N=N, AR=AR) tTestPwr(mu0-mu1, se=sqrt(tmp), df=Inf)
Explicit calculation yields the same:
a <- SteppedPower::glsPower(Cl=rep(6,4), mu0=mu0, mu1=mu1, AR=AR, sigma=0, tau=tau, N=N, psi=psi, verbose=1, INDIV_LVL = TRUE) a
Closed formula for open cohort designs
Kasza et al. offer an extension to open cohort designs by the use of an attrition
or churn rate chi [@kasza2020sample]. It can be understood as the estimated probability
that a subject included in period $j$
is also included in period $j'$ for $j\ne j'$. Note that this differs from the
approach mentioned above, as it results in a mixture of exchangeable covariance matrices.
Since it differs from the way SteppedPower intents to model open cohort designs,
it requires a bit of tinkering to reproduce the results of the closed formula with
this package. It can be done with the help of the random time effect (gamma).
tmp <- VarClosed_Kasza(trtMat, sigma=sigma, tau=tau, gamma=gamma, psi=psi, N=N, chi=0) tTestPwr(mu0-mu1, se=sqrt(tmp), df = Inf) glsPower(Cl = rep(6,4), N=N, mu0=mu0, mu1=mu1, verbose=0, sigma=sigma, tau=tau, gamma=gamma, psi=psi) tmp <- VarClosed_Kasza(trtMat, sigma=sigma, tau=tau, gamma=gamma, psi=psi, N=N, chi=1) tTestPwr(mu0-mu1, sqrt(tmp), df = Inf) glsPower(Cl = rep(6,4), N=N, mu0=mu0, mu1=mu1, verbose=0, sigma=sigma, tau=tau, gamma=sqrt(gamma^2+psi^2/N), psi=0) tmp <- VarClosed_Kasza(trtMat, sigma=sigma, tau=tau, gamma=gamma, psi=psi, N=N, chi=chi) tTestPwr(mu0-mu1, sqrt(tmp), df = Inf) glsPower(Cl = rep(6,4), N=N, mu0=mu0, mu1=mu1, verbose=0, sigma=sigma, tau=tau, gamma=sqrt(gamma^2+chi*psi^2/N), psi=sqrt(1-chi)*psi)
Session Info {-}
print(sessionInfo(),locale=FALSE)
References {-}
Try the SteppedPower package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.