SteppedPowerFor a binomial outcome we keep the model from the Getting Started vignette (Hussey and Hughes 2007; Li et al. 2020), but with an identity link:
\[y_{ijk}= T_{ij} \theta + c_i + \mu_j + e_{ijk}\]
with \(T_{ij}\) the treatment
status, \(\theta\) the treatment
effect, and \(c_i \sim N(0,\tau^2)\) a
random cluster effect. Because the link is the identity, the GLS
machinery of glsPower() applies directly. The key
difference to the Gaussian case is that the residual variance is tied to
the mean:
\[\sigma^2_{ij} = p_{ij}(1-p_{ij}), \qquad p_{ij} = \mu_0 + T_{ij}\,\theta .\]
\(\sigma^2\) is not a free
parameter, but is rather specified by the cell prevalences \(p_{ij}\). Passing sigma for a
binomial family triggers a warning and is ignored.
By default, mu0 and mu1 are the prevalence
conditional on the random effects being zero. With
marginal_mu = TRUE they are read as
marginal prevalences and back-transformed via numerical
integration. This only matters when \(\tau\) is non-negligible.
Note that
glsPower()derives \(\sigma^2\) cell-by-cell and not just once using the mean prevalence \(p_{\text{mid}}=(\mu_0+\mu_1)/2\). This may cause differences to other implementations.
SteppedPower offers two ways to specify the
within-cluster correlation, as described in the
Getting Started vignette. For linear outcomes, the two
routes are equivalent. For binomial outcomes, they differ because the
residual variance changes accross cells.
In this approach, the standard deviations of the random terms are
specified directly: tau (cluster intercept),
gamma (time effect), psi (subject intercept),
eta (random treatment effect), and AR (AR(1)
decay). This route is required when you need eta,
rho (intercept-slope correlation), or AR.
The second implemented approach follows notation from (Li et al. 2018). Here the correlation structure is described by three scalars, \(\alpha_0\), \(\alpha_1\), and \(\alpha_2\), which provide an alternative parameterization of the ICC, CAC, and IAC correlations (see the Getting Started vignette).
Given a marginal variance \(\sigma^2_{\text{marg}}\), these translate into random effects via
\[ \begin{align*} \tau^2 &:= \sigma^2_{\text{marg}} \cdot \alpha_1 \\ \gamma^2 &:= \sigma^2_{\text{marg}} \cdot (\alpha_0 - \alpha_1) \\ \psi^2 &:= \sigma^2_{\text{marg}} \cdot (\alpha_2 - \alpha_1) \\ \sigma^2_{\text{res}} &:= \sigma^2_{\text{marg}} \cdot (1 - \alpha_0 - \alpha_2 + \alpha_1) \end{align*} \]
If alpha_0_1_2 of length 2 is passed, a cross-sectional
design is assumed and \(\alpha_2\) is
set equal to \(\alpha_1\). As with the
ICC/CAC/IAC parameterization, \(\boldsymbol{\alpha}\) cannot specify a
random treatment effect.
Suppose we are planning a stepped wedge trial with four sequences with one cluster each, and we expect a baseline prevalence of \(p=0.2\) and a treatment effect of \(\theta=0.15\). We want to target an ICC of \(0.2\).
Figure 1: Design matrix for 4 sequences with 1 cluster each. 5 periods.
The two calls below target the same ICC of \(0.2\) but yield different power. With \(p=0.2\) and \(\theta=0.15\), the binomial variance ranges
from \(0.16\) (control) to \(0.2275\) (treatment). Under random effects,
tau is a fixed scalar, usually relative to the midpoint
prevalence \(p_{\text{mid}}=0.275\).
Hence, the implied ICC varies from \(0.31\) to \(0.22\) across cells. Under
alpha_0_1_2, tau is recalculated cell by cell
and the ICC stays at \(0.2\)
throughout.
p <- 0.2; effect <- 0.15; icc <- 0.2
p_mid <- p + effect * 0.5
DM <- construct_DesMat(Cl = rep(1,4), N = 50)
pow_re <- glsPower(DesMat = DM, mu0 = p, mu1 = p + effect,
tau = sqrt(p_mid * (1 - p_mid) * icc / (1 - icc)),
family = "binomial", verbose = 2)
#> The assumed odds ratio is 2.1538
pow_al <- glsPower(DesMat = DM, mu0 = p, mu1 = p + effect,
alpha_0_1_2 = c(icc, icc),
family = "binomial", verbose = 2)
#> Since length of alpha_0_1_2 is 2, a cross-sectional design isassumed. Hence, alpha2 is set to alpha1.
#> The assumed odds ratio is 2.1538
c(power_randomEffects = pow_re$power,
power_alpha012 = pow_al$power)
#> power_randomEffects power_alpha012
#> 0.8478889 0.7944555plotly::subplot(
plot(pow_re, which = 4, show_colorbar = FALSE)$CMplot,
plot(pow_al, which = 4, show_colorbar = FALSE)$CMplot,
nrows = 1, margin = 0.05)Figure 2: Covariance matrix under random effects (left) and alpha_0_1_2 (right), both targeting ICC = 0.2.
CAVE: \(\sigma^2_{ij}=p_{ij}(1-p_{ij})\) is cell-specific – \(\sigma^2\) moves with \(p\), so the same \(\tau\) maps to a different ICC in control vs. intervention periods. The alpha parameterization makes this explicit (alphas are ratios to \(\sigma^2_{\text{marg}}\)); random effects keeps \(\tau\) fixed and lets the implied ICC drift.
The direction reverses when \(p\) and \(p_1\) are symmetric around \(0.5\) (here \(p=0.25\), \(p_1=0.75\)), so both arms share the same variance \(\sigma^2=0.1875\) and there is no ICC drift. The difference now comes from the anchoring: random effects derives \(\tau\) from \(p_{\text{mid}}=0.5\), whose variance \(0.25\) exceeds the actual cell variance \(0.1875\), inflating \(\tau\) and lowering power.
The contour plots below shows the difference in treatment effect variance estimators across the \((p_0, p_1)\) plane for the example setting above. The diagonal (\(p_0=p_1\), no treatment effect) is zero by construction.
Figure 3: Relative difference of variance estimators. ICC=0.2, N=50, 4 sequences. blue: midpoint method has higher variance than cell-by-cell method.
Figure 4: Relative difference of variance estimators. ICC=0.2, N=50, 4 sequences. blue: alpha method has higher variance than random effects method.
Count outcomes are handled very similary in
SteppedPower. It uses a Poisson distribution with an
identity link:
\[y_{ijk}= T_{ij} \theta + c_i + \mu_j + e_{ijk}\]
As with the binomial case, the residual variance is tied to the mean via \(\sigma^2 = \mu\):
\[\sigma^2_{ij} = \mu_{ij} = \mu_0 + T_{ij}\,\theta .\] Again, \(\sigma^2\) is derived cell-by-cell. The same two parametrisations of the covariance structure are available with the same implications.