"
In cosa: Bound Constrained Optimal Sample Size Allocation

knitr::opts_chunk$set(echo = TRUE)
library(cosa)

To install and load the package:

install.packages("cosa")
library(cosa)

The following examples demonstrate how to perform bound constrained optimal sample size allocation (BCOSSA) for a two-level cluster randomized trial (CRT) and for the corresponding cluster-level regression discontinuity (CRD) design under three primary constraints. Note: NULL arguments are provided for clarity, otherwise they do not have to be explicit.

Primary Constraint on the Total Cost

# cost constrained - optimize p and n2
# CRT('order = 0' or 'rhots = 0')
crt <- cosa.crd2(order = 0,
                 constrain = "cost", cost = 12500, 
                 cn1 = c(5, 2), cn2 = c(50, 20),
                 es = .20,  power = .80, rho2 = .20,
                 g2 = 1, r21 = .20, r22 = .30, 
                 p = NULL, n1 = 20,  n2 = NULL)
# comparisons to CRDs
# CRD w/ linear score variable interacting with the treatment
crd1 <- cosa.crd2(order = 1, interaction = TRUE,
                  constrain = "cost", cost = 12500,
                  cn1 = c(5, 2), cn2 = c(50, 20),
                  es = .20,  power = .80, rho2 = .20,
                  g2 = 1, r21 = .20, r22 = .30,
                  p = .386, n1 = 24,  n2 = NULL)
# CRD w/ quadratic score variable interacting with the treatment
crd2 <- cosa.crd2(order = 2, interaction = TRUE,
                  constrain = "cost", cost = 12500,
                  cn1 = c(5, 2), cn2 = c(50, 20),
                  es = .20,  power = .80, rho2 = .20,
                  g2 = 1, r21 = .20, r22 = .30,
                  p = .386, n1 = 24,  n2 = NULL)

# example plots
par(mfrow = c(2, 3), mai = c(.6, .6, .6, .2))
# compare minimum detectable effect size and 95% CI
plot(crt, ypar = "mdes", xpar = "n2",
     ylim = c(.10, .90), xlim = c(10, 800), 
     ylab = "MDES (with Power = .80)", xlab = "Number of Clusters",
     main = expression(CRT), locate = TRUE)
plot(crd1, ypar = "mdes", xpar = "n2",
     ylim = c(.10, .90), xlim = c(10, 800),
     ylab = "MDES (with Power = .80)", xlab = "Number of Clusters",
     main = expression(CRD (S + TS)), locate = TRUE)
plot(crd2, ypar = "mdes", xpar = "n2",
     ylim = c(.10, .90), xlim = c(10, 800),
     ylab = "MDES (with Power = .80)", xlab = "Number of Clusters",
     main = expression(CRD (S + S^2 + TS + TS^2)), locate = TRUE)
# compare statistical power 
plot(crt, ypar = "power", xpar = "n2",
     ylim = c(.10, .85), xlim = c(10, 800), 
     ylab = "Power (for ES = .20)", xlab = "Number of Clusters",
     main = expression(CRT), locate = TRUE)
plot(crd1, ypar = "power", xpar = "n2",
     ylim = c(.10, .85), xlim = c(10, 800), 
     ylab = "Power (for ES = .20)", xlab = "Number of Clusters",
     main = expression(CRD(S + TS)), locate = TRUE)
plot(crd2, ypar = "power", xpar = "n2",
     ylim = c(.10, .85), xlim = c(10, 800), 
     ylab = "Power (for ES = .20)", xlab = "Number of Clusters",
     main = expression(CRD (S + S^2 + TS + TS^2)), locate = TRUE)

As seen from the MDES and power plots, a little less than 150 clusters are needed to obtain the benchmark power rate of 80% for the CRT, more than 350 cluster are needed for CRD with linear score variable interacting with the treatment, and more than 650 clusters are needed for CRD with quadratic score variable interacting with the treatment. Precise number of clusters can be found via placing the primary constraint on either the effect size or the power rate.

Primary Constraint on the Effect Size

# cost constrained - optimize p and n2
# CRT('order = 0' or 'rhots = 0')
cosa.crd2(order = 0,
          constrain = "es", es = .20,  
          cn1 = c(5, 2), cn2 = c(50, 20),
          power = .80, rho2 = .20,
          g2 = 1, r21 = .20, r22 = .30, 
          p = NULL, n1 = 20,  n2 = NULL)
# comparisons to CRDs
# CRD w/ linear score variable interacting with the treatment
cosa.crd2(order = 1, interaction = TRUE,
          constrain = "es", es = .20,  
          cn1 = c(5, 2), cn2 = c(50, 20),
          power = .80, rho2 = .20,
          g2 = 1, r21 = .20, r22 = .30,
          p = .389, n1 = 24,  n2 = NULL)
# CRD w/ quadratic score variable interacting with the treatment
cosa.crd2(order = 2, interaction = TRUE,
          constrain = "es", es = .20,  
          cn1 = c(5, 2), cn2 = c(50, 20),
          power = .80, rho2 = .20,
          g2 = 1, r21 = .20, r22 = .30,
          p = .389, n1 = 24,  n2 = NULL)

Primary Constraint on the Power Rate

# cost constrained - optimize p and n2
# CRT('order = 0' or 'rhots = 0')
cosa.crd2(order = 0,
          constrain = "power", power = .80, 
          cn1 = c(5, 2), cn2 = c(50, 20),
          es = .20, rho2 = .20,
          g2 = 1, r21 = .20, r22 = .30, 
          p = NULL, n1 = 20,  n2 = NULL)
# comparisons to CRDs
# CRD w/ linear score variable interacting with the treatment
cosa.crd2(order = 1, interaction = TRUE,
          constrain = "power", power = .80, 
          cn1 = c(5, 2), cn2 = c(50, 20),
          es = .20, rho2 = .20,
          g2 = 1, r21 = .20, r22 = .30,
          p = .389, n1 = 24,  n2 = NULL)
# CRD w/ quadratic score variable interacting with the treatment
cosa.crd2(order = 2, interaction = TRUE,
          constrain = "power", power = .80,  
          cn1 = c(5, 2), cn2 = c(50, 20),
          es = .20, rho2 = .20,
          g2 = 1, r21 = .20, r22 = .30,
          p = .389, n1 = 24,  n2 = NULL)