--- title: "Binomial and Count Outcomes in `SteppedPower`" author: - name: Philipp Mildenberger affiliation: Institute of Medical Biostatistics, Epidemiology and Informatics ([IMBEI, Mainz](https://www.unimedizin-mainz.de/imbei/imbei/welcome-page/)) date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true number_sections: true bibliography: references.bib vignette: > %\VignetteIndexEntry{Binomial and Count Outcomes in `SteppedPower`} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, warning = FALSE ) library(knitr) library(SteppedPower) library(Matrix) library(plotly) assign("knit_print.plotly", function(x, ...) { x <- plotly::partial_bundle(x) NextMethod() }, envir = .GlobalEnv) # A function for captioning and referencing images fig <- local({ i <- 0 ref <- list() list( cap=function(refName, text) { i <<- i + 1 ref[[refName]] <<- i paste("Figure ", i, ": ", text, sep="") }, ref=function(refName) { paste("(Fig.", ref[[refName]],")") }) }) ``` # Methods and Notation For a binomial outcome we keep the model from the [Getting Started](Getting_Started.html) vignette [@hussey2007design; @li2020mixed], 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. # Covariance Structure Specifications `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. ## Random effects 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`. ## $\alpha_0$, $\alpha_1$, $\alpha_2$ The second implemented approach follows notation from [@li2018sample]. 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](Getting_Started.html) 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. # Example Usage :::::: {.columns} ::: {.column width="40%"} 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$. ::: ::: {.column width="60%"} ```{r design-plot, fig.cap=fig$cap("design","Design matrix for 4 sequences with 1 cluster each. 5 periods."), fig.width=2, fig.height=1, echo=FALSE} dsn_bin <- construct_DesMat(Cl = rep(1, 4), N = 50) plot(dsn_bin) ``` ::: :::::: 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. ```{r} 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) pow_al <- glsPower(DesMat = DM, mu0 = p, mu1 = p + effect, alpha_0_1_2 = c(icc, icc), family = "binomial", verbose = 2) c(power_randomEffects = pow_re$power, power_alpha012 = pow_al$power) ``` ```{r, fig.cap=fig$cap("covmat","Covariance matrix under random effects (left) and alpha_0_1_2 (right), both targeting ICC = 0.2.")} plotly::subplot( plot(pow_re, which = 4, show_colorbar = FALSE)$CMplot, plot(pow_al, which = 4, show_colorbar = FALSE)$CMplot, nrows = 1, margin = 0.05) ``` >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. ## Differences in the Variance estimator $\text{Var}(\hat\theta)$ the $(p_0, p_1)$ Plane 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. ```{r variance computation, echo=FALSE} # data-raw/binom_outcomes_grid.R for the computation. # vals <- readRDS("inst/vignettes/BinomOutcomes_vals.rds") vals <- readRDS(system.file("vignettes", "BinomOutcomes_vals.rds", package = "SteppedPower")) p_vals <- unique(vals$p0) reMat <- matrix(vals[, "re"], length(p_vals), length(p_vals)) alMat <- matrix(vals[, "al"], length(p_vals), length(p_vals)) linMat <- matrix(vals[, "lin"], length(p_vals), length(p_vals)) ``` ### Cell-by-cell calculation of $\sigma^2$ versus midpoint method ```{r contour_diff, fig.width=6.5, fig.height=5.5, echo=FALSE, fig.cap=fig$cap("contour_midpoint","Relative difference of variance estimators. ICC=0.2, N=50, 4 sequences. blue: midpoint method has higher variance than cell-by-cell method.")} diff_mat <- (reMat - linMat)/linMat * 100 maxabs <- max(abs(diff_mat), na.rm = TRUE) plot_ly( x = p_vals, y = p_vals, z = diff_mat, type = "contour", colorscale = list(list(0, "steelblue"), list(0.5, "white"), list(1, "firebrick") ), zmin = -maxabs, zmax = maxabs, contours = list(coloring = "fill", showlabels = TRUE), line = list(smoothing = 1.3) ) |> add_trace(x = range(p_vals), y = range(p_vals), type = "scatter", mode = "lines", line = list(color = "gray", dash = "dash"), showlegend = FALSE, hoverinfo = "none") |> layout( xaxis = list(title = "p0"), yaxis = list(title = "p1") ) ``` ### Random effects versus alpha notation ```{r contour_diff2, fig.width=6.5, fig.height=5.5, echo=FALSE, fig.cap=fig$cap("contour_alpha","Relative difference of variance estimators. ICC=0.2, N=50, 4 sequences. blue: alpha method has higher variance than random effects method.")} diff_mat <- (reMat - alMat)/alMat * 100 maxabs <- max(abs(diff_mat), na.rm = TRUE) plot_ly( x = p_vals, y = p_vals, z = diff_mat, type = "contour", colorscale = list(list(0, "steelblue"), list(0.5, "white"), list(1, "firebrick") ), zmin = -maxabs, zmax = maxabs, contours = list(coloring = "fill", showlabels = TRUE), line = list(smoothing = 1.3)) |> add_trace(x = range(p_vals), y = range(p_vals), type = "scatter", mode = "lines", line = list(color = "gray", dash = "dash"), showlegend = FALSE, hoverinfo = "none") |> layout( xaxis = list(title = "p0"), yaxis = list(title = "p1") ) ``` # Count (Poisson) Outcomes 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. # References {-}