# Bound Constrained Optimal Design of Multilevel Regression Discontinuity Designs and Randomized Controlled Trials

#### 2019-08-24

To install and load the package:

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

The following examples demonstrate how to perform bound constrained optimal sample allocation (BCOSA) 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.crd2r2(order = 0,
constrain = "cost", cost = 12500,
cn1 = c(5, 2), cn2 = c(50, 20),
es = .20,  power = .80, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = NULL, n1 = 20,  n2 = NULL)
## Solution converged with LBFGS
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2     p [cost]  mdes 95%lcl 95%ucl power
##    20 132 0.386  12510 0.209  0.063  0.356 0.763
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.209 (with power = 80)
## power = 0.763 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## Random assignment design w/o score
# comparisons to CRDs
# CRD w/ linear score variable
crd1 <- cosa.crd2r2(order = 1,
constrain = "cost", cost = 12500,
cn1 = c(5, 2), cn2 = c(50, 20),
es = .20,  power = .80, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = .386, n1 = 24,  n2 = NULL)
## Solution converged with LBFGS
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2   [p] [cost]  mdes 95%lcl 95%ucl power
##    24 116 0.388  12478 0.356  0.106  0.605 0.351
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.356 (with power = 80)
## power = 0.351 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## RD design w/ linear functional form
# CRD w/ linear + quadratic score variable
crd2 <- cosa.crd2r2(order = 2,
constrain = "cost", cost = 12500,
cn1 = c(5, 2), cn2 = c(50, 20),
es = .20,  power = .80, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = .386, n1 = 24,  n2 = NULL)
## Solution converged with LBFGS
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2   [p] [cost]  mdes 95%lcl 95%ucl power
##    24 116 0.388  12478 0.368   0.11  0.627 0.331
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.368 (with power = 80)
## power = 0.331 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## RD design w/ linear + quadratic functional form
# 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, .50), xlim = c(10, 200),
ylab = "MDES (with Power = .80)", xlab = "Number of Clusters",
main = expression(CRT), locate = TRUE)
plot(crd1, ypar = "mdes", xpar = "n2",
ylim = c(.10, .50), xlim = c(10, 200),
ylab = "MDES (with Power = .80)", xlab = "Number of Clusters",
main = expression(CRD(S)), locate = TRUE)
plot(crd2, ypar = "mdes", xpar = "n2",
ylim = c(.10, .50), xlim = c(10, 200),
ylab = "MDES (with Power = .80)", xlab = "Number of Clusters",
main = expression(CRD (S + S^2)), locate = TRUE)
# compare statistical power
plot(crt, ypar = "power", xpar = "n2",
ylim = c(.10, .85), xlim = c(10, 200),
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, 200),
ylab = "Power (for ES = .20)", xlab = "Number of Clusters",
main = expression(CRD(S)), locate = TRUE)
plot(crd2, ypar = "power", xpar = "n2",
ylim = c(.10, .85), xlim = c(10, 200),
ylab = "Power (for ES = .20)", xlab = "Number of Clusters",
main = expression(CRD (S + S^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 200 clusters for the CRD). Precise number of clusters can be found via placing the primary constraint on either the effect size or power rate.

## Primary Constraint on the Effect Size

# cost constrained - optimize p and n2
# CRT('order = 0' or 'rhots = 0')
cosa.crd2r2(order = 0,
constrain = "es", es = .20,
cn1 = c(5, 2), cn2 = c(50, 20),
power = .80, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = NULL, n1 = 20,  n2 = NULL)
## Solution converged with SLSQP
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2     p  cost [mdes] 95%lcl 95%ucl power
##    20 144 0.389 13680    0.2   0.06   0.34   0.8
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.2 (with power = 80)
## power = 0.8 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## Random assignment design w/o score
# comparisons to CRDs
# CRD w/ linear score variable
cosa.crd2r2(order = 1,
constrain = "es", es = .20,
cn1 = c(5, 2), cn2 = c(50, 20),
power = .80, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = .386, n1 = 24,  n2 = NULL)
## Solution converged with SLSQP
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2   [p]  cost [mdes] 95%lcl 95%ucl power
##    24 363 0.386 38964    0.2   0.06   0.34   0.8
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.2 (with power = 80)
## power = 0.8 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## RD design w/ linear functional form
# CRD w/ linear + quadratic score variable
cosa.crd2r2(order = 2,
constrain = "es", es = .20,
cn1 = c(5, 2), cn2 = c(50, 20),
power = .80, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = .386, n1 = 24,  n2 = NULL)
## Solution converged with SLSQP
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2   [p]  cost [mdes] 95%lcl 95%ucl power
##    24 389 0.386 41752    0.2   0.06   0.34   0.8
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.2 (with power = 80)
## power = 0.8 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## RD design w/ linear + quadratic functional form

## Primary Constraint on the Power Rate

# cost constrained - optimize p and n2
# CRT('order = 0' or 'rhots = 0')
cosa.crd2r2(order = 0,
constrain = "power", power = .80,
cn1 = c(5, 2), cn2 = c(50, 20),
es = .20, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = NULL, n1 = 20,  n2 = NULL)
## Solution converged with SLSQP
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2     p  cost mdes 95%lcl 95%ucl [power]
##    20 144 0.389 13680  0.2   0.06   0.34     0.8
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.2 (with power = 80)
## power = 0.8 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## Random assignment design w/o score
# comparisons to CRDs
# CRD w/ linear score variable
cosa.crd2r2(order = 1,
constrain = "power", power = .80,
cn1 = c(5, 2), cn2 = c(50, 20),
es = .20, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = .386, n1 = 24,  n2 = NULL)
## Solution converged with SLSQP
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2   [p]  cost mdes 95%lcl 95%ucl [power]
##    24 363 0.386 38964  0.2   0.06   0.34     0.8
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.2 (with power = 80)
## power = 0.8 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## RD design w/ linear functional form
# CRD w/ linear + quadratic score variable
cosa.crd2r2(order = 2,
constrain = "power", power = .80,
cn1 = c(5, 2), cn2 = c(50, 20),
es = .20, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = .386, n1 = 24,  n2 = NULL)
## Solution converged with SLSQP
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2   [p]  cost mdes 95%lcl 95%ucl [power]
##    24 389 0.386 41752  0.2   0.06   0.34     0.8
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.2 (with power = 80)
## power = 0.8 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## RD design w/ linear + quadratic functional form

## Balanced Case in CRTs

# fix 'p = .50' to inpsect cost for the balanced case
cosa.crd2r2(order = 0,
constrain = "power", power = .80,
cn1 = c(5, 2), cn2 = c(50, 20),
es = .20, rho2 = .20,
g2 = 5, r21 = .20, r22 = .30,
p = .50, n1 = 20,  n2 = NULL)
## Solution converged with SLSQP
##
## Rounded solution:
## ---------------------------------------------------
##  [n1]  n2   [p]  cost mdes 95%lcl 95%ucl [power]
##    20 137 0.496 14340  0.2   0.06   0.34     0.8
## ---------------------------------------------------
## Per unit marginal costs:
##  Level 1 treatment: 5
##  Level 1 control: 2
##  Level 2 treatment: 50
##  Level 2 control: 20
## ---------------------------------------------------
## MDES = 0.2 (with power = 80)
## power = 0.8 (for ES = 0.2)
## ---------------------------------------------------
## []: point constrained (fixed)
## <<: bound constrained
## Random assignment design w/o score

–o–