--- title: "Stratified Designs" author: "Gosuke Homma" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{Stratified Designs} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4, fig.align = "center" ) ``` ```{r load} library(bbssr) ``` ## Design and notation A trial is stratified into $J \ge 2$ strata, for example by centre or by a baseline factor, and stratum $j$ has $N_{1j}$ patients in group 1 and $N_{2j}$ in group 2, both fixed by the design (random sizes are treated in [Random stratum and group sizes](#random-stratum-and-group-sizes)). Let $x_{1j}$ and $x_{2j}$ be the numbers of responders of the two groups and $z_j = x_{1j} + x_{2j}$ the number of responders of the stratum. All tests of the `Strat*` functions use $S = \sum_j x_{1j}$, the number of responders of group 1 summed over the strata. Given $z = (z_1, \ldots, z_J)$, the null expectation of $S$ is $E = \sum_j N_{1j} z_j / N_j$ with $N_j = N_{1j} + N_{2j}$, and each test rejects when $S$ is at most a lower or at least an upper critical value. The functions follow those for two independent groups: `StratRR()`, `StratPower()` and `StratSampleSize()` for the fixed design, and `StratPowerBSSR()`, `StratGridBSSR()`, `StratBSSR()`, `StratTypeIErrorBSSR()` and `StratAlphaAdjBSSR()` for blinded sample size re-estimation. ## Tests Five stratified tests are available. `'CMH'` is the Cochran-Mantel-Haenszel statistic $(S - E) / \sqrt{V}$ with the hypergeometric variance $V = \sum_j N_{1j} N_{2j} z_j (N_j - z_j) / \{N_j^2 (N_j - 1)\}$, and `'Cochran'` is Cochran's statistic with the binomial variance $\sum_j N_{1j} N_{2j} z_j (N_j - z_j) / N_j^3$, the statistic of Woolson et al. (1986). `'CMH-cc'` and `'Cochran-cc'` reduce $|S - E|$ by 0.5 when it is at least 0.5, the continuity correction of `stats::mantelhaen.test()`. Nam (1992) reduces $S - E$ by 0.5 in every outcome; the two corrections reject the same outcomes at one-sided levels up to 0.078, since an outcome with $|S - E| < 0.5$ has a larger one-sided p-value. `'Strat-Fisher'` is the stratified exact conditional test of Jung (2014), which refers $S$ to its distribution given $z$, the convolution of the hypergeometric distributions of the strata. Three unstratified tests, which analyse the table of totals and ignore the strata, are available for comparison (see [Analyses that ignore the strata](#analyses-that-ignore-the-strata)). Jung (2014, Table 2) analyses the response to thymosin in three strata of patients with bronchogenic carcinoma. ```{r jung-example} thymosin <- array(c(10, 12, 1, 1, 9, 11, 0, 1, 8, 7, 0, 3), c(2, 2, 3)) p.exact <- stats::mantelhaen.test(thymosin, alternative = 'greater', exact = TRUE)$p.value p.mh <- stats::mantelhaen.test(thymosin, alternative = 'greater', correct = FALSE)$p.value round(c(exact = p.exact, mh = p.mh), 4) ``` The one-sided p-values `r sprintf('%.4f', p.exact)` and `r sprintf('%.4f', p.mh)` are those of Jung (2014) for the stratified exact test and the Mantel-Haenszel test, which `'Strat-Fisher'` and `'CMH'` compute in the same way. `StratRR()` lists the critical values of a design for every vector of stratum totals. ```{r rr} rr <- StratRR(N1 = c(10, 12), N2 = c(10, 12), alpha = 0.025, Test = 'Strat-Fisher') print(rr, n = 6) ``` ## Power and sample size The power is the probability of the rejection region when the numbers of responders are binomial in every group and stratum. It is computed exactly for any number of strata and any response probabilities, so the odds ratio need not be common to the strata. The normal approximation to the total sample size is formula (2.8) of Woolson et al. (1986), and with a continuity correction formula (2.3) of Nam (1992). A design is described by the fractions of the strata and the allocation ratio, and the exact search runs over the total sample size within this family. ```{r sample-size} p2 <- c(0.3, 0.5) p1 <- p2 * 2.5 / (1 - p2 + p2 * 2.5) ss.std <- StratSampleSize(p1, p2, frac = c(0.4, 0.6), alpha = 0.025, tar.power = 0.8, Test = 'CMH', method = 'standard') ss.cmh <- StratSampleSize(p1, p2, frac = c(0.4, 0.6), alpha = 0.025, tar.power = 0.8, Test = 'CMH') ss.exact <- StratSampleSize(p1, p2, frac = c(0.4, 0.6), alpha = 0.025, tar.power = 0.8, Test = 'Strat-Fisher') c(standard = ss.std$N, CMH = ss.cmh$N, exact = ss.exact$N) ``` With a common odds ratio of 2.5, the normal approximation gives a total of `r ss.std$N` patients, at which the exact power of the CMH test is `r sprintf('%.3f', ss.std$Power)`. The exact search gives `r ss.cmh$N` patients for the CMH test and `r ss.exact$N` for the stratified exact test, whose conditional control of the type I error rate makes it conservative. ## Blinded sample size re-estimation At an interim analysis the numbers of responders of every stratum can be counted without unblinding. The blinded estimate of the pooled response probability of stratum $j$ is $s_j$ divided by the interim size of the stratum. It is divided between the groups with the assumed effect `Delta.A`, common to the strata on the scale `effect`, and the total sample size is re-estimated for these probabilities and divided between the strata as in the planned design. With two strata the final sizes depend on the interim data only through $(s_1, s_2)$, and the power is computed exactly by enumeration. With a common odds ratio the planning probabilities of a stratum whose pooled estimate is close to 0 or 1 are close to each other, so the examples bound the final total by twice the planned total with `N.max`. ```{r bssr} plan <- StratSampleSize(p1, p2, frac = c(0.5, 0.5), alpha = 0.025, tar.power = 0.8, Test = 'CMH', method = 'standard') pooled <- (p1 + p2) / 2 res <- StratPowerBSSR(p = rbind(pooled, pooled + 0.1), Delta.A = 2.5, Delta.T = 2.5, N1 = plan$N1, N2 = plan$N2, omega = 0.5, r = 1, alpha = 0.025, tar.power = 0.8, Test = 'CMH', effect = 'OR', ss.method = 'standard', N.max = 2 * plan$N) res ``` When the pooled response probabilities are those of the planning (the first row), the power of the re-estimation design is `r sprintf('%.3f', res$power.BSSR[1])` with an expected total of `r sprintf('%.1f', res$E.N[1])`. When they are 0.1 higher in both strata, the re-estimation adjusts the sample size and the power is `r sprintf('%.3f', res$power.BSSR[2])`, while that of the fixed design is `r sprintf('%.3f', res$power.TRAD[2])`. ## Type I error rate Under the null hypothesis both groups of stratum $j$ respond with the probability $\theta_j$. Given the interim and the second-stage numbers of responders of a stratum, the numbers of responders of group 1 are hypergeometric whatever $\theta_j$, so with two strata the type I error rate is a polynomial in $(\theta_1, \theta_2)$. `StratTypeIErrorBSSR()` computes its coefficients in the products of the Bernstein bases of the two strata, evaluates it on a grid, and certifies its largest value over the rectangle spanned by the grid: the rectangle is subdivided until the largest coefficient of every piece exceeds the largest value found by at most $10^{-12}$. ```{r tie} tie <- StratTypeIErrorBSSR(Delta.A = 2.5, N1 = c(20, 20), N2 = c(20, 20), omega = 0.5, r = 1, alpha = 0.025, tar.power = 0.8, Test = 'CMH', effect = 'OR', ss.method = 'standard', N.max = 160, theta = seq(0.05, 0.95, by = 0.05)) tie mx <- attr(tie, 'max') ``` The largest type I error rate of the re-estimation design with 20 patients per group and stratum is `r sprintf('%.5f', mx$TIE[1])`, at $(\theta_1, \theta_2) =$ (`r sprintf('%.3f', mx$theta.S1[1])`, `r sprintf('%.3f', mx$theta.S2[1])`), and no rate in the rectangle exceeds `r sprintf('%.5f', mx$bound[1])`. That of the fixed design is `r sprintf('%.5f', mx$TIE[2])`. ```{r tie-plot} plot(tie) ``` `StratAlphaAdjBSSR()` lowers the nominal level until the certified largest rate does not exceed `alpha`. ```{r alpha-adj} adj <- StratAlphaAdjBSSR(Delta.A = 2.5, N1 = c(20, 20), N2 = c(20, 20), omega = 0.5, r = 1, alpha = 0.025, tar.power = 0.8, Test = 'CMH', effect = 'OR', ss.method = 'standard', N.max = 160, theta = seq(0.05, 0.95, by = 0.05)) adj ``` `StratCondRejectBSSR()` gives, under the null hypothesis, the probability of rejection given the pooled numbers of responders of the strata at the interim analysis and in the second stage. These probabilities do not depend on $(\theta_1, \theta_2)$, and averaging them over the second stage gives the type I error rate given the interim outcome. Below they are evaluated at the point of the largest rate found above. ```{r crp} crp <- StratCondRejectBSSR(Delta.A = 2.5, N1 = c(20, 20), N2 = c(20, 20), omega = 0.5, r = 1, alpha = 0.025, tar.power = 0.8, Test = 'CMH', effect = 'OR', ss.method = 'standard', N.max = 160, theta = cbind(mx$theta.S1[1], mx$theta.S2[1])) crp by.s <- attr(crp, 'by.s') head(by.s[order(-by.s$contribution), ], 5) ``` The contributions of the interim outcomes sum to `r sprintf('%.5f', attr(crp, 'TIE')$TIE)`, the largest type I error rate found above. ## Analyses that ignore the strata `'Unstrat-Chisq'` is Pearson's chi-squared test of the table of totals, with $\sum_j N_{1j}$ and $\sum_j N_{2j}$ patients and $\sum_j z_j$ responders, which is the comparison of two proportions with the pooled variance; `'Unstrat-Chisq-cc'` adds the continuity correction, and `'Unstrat-Fisher'` is Fisher's exact test of that table. Their critical values depend on $z$ only through $\sum_j z_j$, but their power and type I error rate are computed with the response probabilities of each stratum. An unstratified test controls the level when the strata are distributed alike in the two groups, and its level can be distorted otherwise (Nam, 1992; Jung, 2014). In the next design group 1 holds a third of the first stratum, whose patients respond rarely, and three quarters of the second. ```{r unstrat} th <- rbind(c(0.2, 0.6), c(0.4, 0.4)) lev <- sapply(c('CMH', 'Strat-Fisher', 'Unstrat-Chisq', 'Unstrat-Fisher'), function(Test) { StratPower(th, th, N1 = c(20, 60), N2 = c(40, 20), alpha = 0.025, Test = Test)$Power }) rownames(lev) <- c('theta = (0.2, 0.6)', 'theta = (0.4, 0.4)') round(lev, 4) ``` With the response probabilities 0.2 and 0.6 in the two strata, the one-sided type I error rate of the unstratified chi-squared test is `r sprintf('%.3f', lev[1, 'Unstrat-Chisq'])` instead of 0.025, because the table of totals compares groups made of different mixtures of the strata. When the strata respond alike, the rate is `r sprintf('%.4f', lev[2, 'Unstrat-Chisq'])`. The re-estimation can ignore the strata as well: with `ss.strata = 'pooled'` the sample size is re-estimated from the proportion of responders among all the interim patients, used for every stratum, so the final sizes depend only on the total number of interim responders. The final analysis is still the test given by `Test`. ```{r bssr-pooled} res.pooled <- StratPowerBSSR(p = rbind(pooled, pooled + 0.1), Delta.A = 2.5, Delta.T = 2.5, N1 = plan$N1, N2 = plan$N2, omega = 0.5, r = 1, alpha = 0.025, tar.power = 0.8, Test = 'CMH', effect = 'OR', ss.method = 'standard', N.max = 2 * plan$N, ss.strata = 'pooled') rbind(stratified = res$E.N, pooled = res.pooled$E.N) ``` The rows give the expected final totals of the two estimates in the two scenarios of the previous section. ## Random stratum and group sizes When the patients are enrolled one at a time, the size of each stratum is known only at the end of the enrolment, and the sizes of the groups within a stratum depend on the randomization. With two strata, `StratPower()` describes such a design by its total `N`, the fractions of the strata `frac` and the allocation ratio `r`, `StratSampleSize()` searches the total for the same description, and `alloc` chooses how the sizes are drawn. With `alloc = 'stratum'` the size of stratum 1 is binomial with the probability `frac[1]`, and within each stratum the groups are allocated by permuted blocks of size `block`, as in stratified block randomization. With `'group'` the size of stratum 1 is fixed at `floor(N * frac[1])`, and each patient is in group 1 with the probability $r / (1 + r)$ (Jung, 2014, Section 3.2, where the size of stratum 1 is the round-off of `N * frac[1]`; the integer part used here reproduces all but one of the sample sizes of the column 'n fixed' of Table 4 of the article). With `'both'` the stratum sizes are binomial as well (Jung, 2014, Section 3.1). The power is the sum over the stratum and group sizes of their probability times the power of the fixed design they give. Sizes whose probability is below $10^{-14}$ are left out, and the attribute `skipped` gives their total probability, a bound on the error of the power. ```{r alloc} alloc.power <- function(alloc) { StratPower(p1, p2, alpha = 0.025, Test = 'CMH', alloc = alloc, N = 80, frac = c(0.5, 0.5))$Power } pw <- c(fixed = StratPower(p1, p2, N1 = c(20, 20), N2 = c(20, 20), alpha = 0.025, Test = 'CMH')$Power, stratum = alloc.power('stratum'), group = alloc.power('group')) round(pw, 4) ``` For the design of the section on the type I error rate, with 20 patients in each group and stratum and 80 in all, the power of the CMH test with the planning probabilities is `r sprintf('%.4f', pw[['fixed']])` with fixed group sizes, `r sprintf('%.4f', pw[['stratum']])` with permuted blocks in random strata and `r sprintf('%.4f', pw[['group']])` with random groups in fixed strata. Permuted blocks keep the power within 0.002 of that with fixed sizes, while binomial group sizes, which vary around the planned sizes, lower it by about 0.02. The number of configurations of the sizes grows most quickly with `'both'`, whose power takes far longer to compute than with the other models. `StratPowerBSSR()`, `StratTypeIErrorBSSR()`, `StratAlphaAdjBSSR()` and `StratBSSR()` take `alloc` and `block` as well. With two strata the interim analysis has the first patients enrolled, as many as the interim total of the design, and the final analysis the first `N` patients, where `N` is the total re-estimated from the interim data with the fractions of the strata among the interim patients (`ss.frac = 'observed'`, the default) or with those of the plan (`'planned'`). Given an allocation sequence, that is the strata and the groups of the patients in the order of enrolment, the design has fixed group sizes and is computed exactly. The power and the type I error rate are averaged over `nseq` sequences drawn with the random numbers of R, and `se.BSSR` is the Monte Carlo standard error of the average. The power of the fixed design with the same allocation is exact. ```{r bssr-alloc} bssr.small <- function(...) { StratPowerBSSR(p = rbind(pooled, pooled + 0.1), Delta.A = 2.5, Delta.T = 2.5, N1 = c(20, 20), N2 = c(20, 20), omega = 0.5, r = 1, alpha = 0.025, tar.power = 0.8, Test = 'CMH', effect = 'OR', ss.method = 'standard', N.max = 160, ...) } res.small <- bssr.small() set.seed(2) res.alloc <- bssr.small(alloc = 'stratum', nseq = 20) res.small[, c('power.BSSR', 'power.TRAD', 'E.N')] res.alloc[, c('power.BSSR', 'se.BSSR', 'power.TRAD', 'E.N')] ``` With fixed group sizes the re-estimation raises the power from `r sprintf('%.4f', res.small$power.TRAD[1])`, that of the fixed design with 80 patients, to `r sprintf('%.4f', res.small$power.BSSR[1])` in the first scenario, with an expected total of `r sprintf('%.1f', res.small$E.N[1])` close to the bound of 160. With permuted blocks in random strata and 20 sequences, the power of the re-estimation design is `r sprintf('%.4f', res.alloc$power.BSSR[1])` and `r sprintf('%.4f', res.alloc$power.BSSR[2])` in the two scenarios, with the standard errors `r sprintf('%.5f', res.alloc$se.BSSR[1])` and `r sprintf('%.5f', res.alloc$se.BSSR[2])`, against `r sprintf('%.4f', res.small$power.BSSR[1])` and `r sprintf('%.4f', res.small$power.BSSR[2])` with fixed group sizes, and the expected totals are `r sprintf('%.1f', res.alloc$E.N[1])` and `r sprintf('%.1f', res.alloc$E.N[2])`, against `r sprintf('%.1f', res.small$E.N[1])` and `r sprintf('%.1f', res.small$E.N[2])`. For this design the random sizes change the power of the re-estimation design by less than 0.01 and its expected total by less than one patient. The default `nseq = 100` divides the standard errors by about 2.2. ## Several designs `StratGridBSSR()` evaluates several re-estimation designs over common scenarios, as `BinaryGridBSSR()` does for two groups (see `vignette('bbssr-design-grid')`). Each row of a data frame is a design, whose columns are arguments of `StratPowerBSSR()`; an argument with one value per stratum is given by the columns with the suffixes `.S1`, `.S2`, and so on, or as a vector shared by all designs. Below, two tests and two timings of the interim analysis are compared for the design of the section on the type I error rate, with the pooled response probability of stratum 2 varying between the scenarios. ```{r grid} design <- expand.grid(Test = c('CMH', 'Cochran'), omega = c(0.3, 0.5), stringsAsFactors = FALSE) q <- cbind(0.4, c(0.5, 0.6, 0.7)) grid <- StratGridBSSR(design, p = q, Delta.A = 2.5, effect = 'OR', N1 = c(20, 20), N2 = c(20, 20), r = 1, alpha = 0.025, tar.power = 0.8, ss.method = 'standard', N.max = 160) summary(grid)[, c('Test', 'omega', 'power.BSSR.min', 'power.BSSR.max', 'E.N.max')] ``` The summary gives the range of the power of each design over the scenarios and its largest expected total sample size. The plot draws the power against the probability that varies between the scenarios, here that of stratum 2, with the fixed design as dashed lines. ```{r grid-plot} plot(grid, colour.by = 'Test', facet.by = 'omega') ``` ## Three or more strata With three or more strata the number of interim outcomes grows as a product over the strata, and `StratPowerBSSR()` and `StratTypeIErrorBSSR()` estimate the rates of the re-estimation design by simulation, with the random numbers of R and a Monte Carlo standard error. With fixed group sizes the rates of the fixed design remain exact. With `alloc = 'stratum'` or `'both'` the strata and the groups of the patients are simulated with the responses, in the fixed design as well, whose standard error is `se.TRAD`. The model `'group'` needs two strata. ```{r three-strata} set.seed(1) q2 <- c(0.2, 0.35, 0.5) q1 <- q2 * 2.5 / (1 - q2 + q2 * 2.5) plan3 <- StratSampleSize(q1, q2, frac = c(0.3, 0.3, 0.4), alpha = 0.025, tar.power = 0.8, Test = 'CMH', method = 'standard') res3 <- StratPowerBSSR(p = (q1 + q2) / 2, Delta.A = 2.5, Delta.T = 2.5, N1 = plan3$N1, N2 = plan3$N2, omega = 0.5, r = 1, alpha = 0.025, tar.power = 0.8, Test = 'CMH', effect = 'OR', ss.method = 'standard', N.max = 2 * plan3$N, nsim = 2e4) res3[, c('power.BSSR', 'se.BSSR', 'power.TRAD', 'E.N')] ``` ## References Jung, S. H. (2014). Stratified Fisher's exact test and its sample size calculation. *Biometrical Journal*, 56, 129-140. Nam, J. (1992). Sample size determination for case-control studies and the comparison of stratified and unstratified analyses. *Biometrics*, 48, 389-395. Woolson, R. F., Bean, J. A. and Rojas, P. B. (1986). Sample size for case-control studies using Cochran's statistic. *Biometrics*, 42, 927-932.