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

Metin Bulus & Nianbo Dong

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–