# R/CACEnocov.R In experiment: R Package for Designing and Analyzing Randomized Experiments

#### Defines functions covClusterCACEnocov

```###
### Calculate the CACE for individual and group randomized trials
###

CACEnocov <- function(Y, D, Z, data = parent.frame(), grp = NULL,
match = NULL, grp.size = NULL,
method = "standard"){

## get the data
call <- match.call()
Y <- eval(call\$Y, envir = data)
D <- eval(call\$D, envir = data)
Z <- eval(call\$Z, envir = data)
grp <- eval(call\$grp, envir = data)
grp.size <- eval(call\$grp.size, envir = data)
match <- eval(call\$match, envir = data)

## point estimates and variances
ITTY <- ATEcluster(Y = Y, Z = Z, grp = grp, match = match)
#                   grp.size = grp.size)
ITTD <- ATEcluster(Y = D, Z = Z, grp = grp, match = match)
#                   grp.size = grp.size)
CACEest <- ITTY\$est/ITTD\$est

## covariance calculation
if (is.null(grp) && is.null(grp.size)) { # unit randomization
if (is.null(match)) { # without matching
Cov <- cov(Y[Z==0], D[Z==0])/sum(Z==0)+
cov(Y[Z==1], D[Z==1])/sum(Z==1)
CACEvar <- (ITTY\$var*(ITTD\$est^2) + ITTD\$var*(ITTY\$est^2) -
2*Cov*ITTY\$est*ITTD\$est)/(ITTD\$est^4)
} else { # with matching
K <- length(ITTY\$diff)
Cov <- cov(ITTY\$diff, ITTD\$diff)/K
CACEvar <- (var(ITTY\$diff)*(ITTD\$est^2)/K +
var(ITTD\$diff)*(ITTY\$est^2)/K -
2*Cov**ITTY\$est*ITTD\$est)/(ITTD\$est^4)
}
} else if (method %in% c("weighted", "unpooled")) {
if (is.null(match)) { # without matching
if (method == "unpooled") {
Cov <- covCluster(Y[Z==0], D[Z==0], grp[Z==0])\$cov +
covCluster(Y[Z==1], D[Z==1], grp[Z==1])\$cov
} else {
Cov <- covCluster(Y, D, grp)\$cov
}
CACEvar <- (ITTY\$var*(ITTD\$est^2) + ITTD\$var*(ITTY\$est^2) -
2*Cov*ITTY\$est*ITTD\$est)/(ITTD\$est^4)
} else { # with matching
M <- length(ITTY\$diff)
N <- ITTY\$N
Cov <- M*sum((ITTY\$diff*ITTY\$weights - N*ITTY\$est/M) *
(ITTD\$diff*ITTD\$weights - N*ITTD\$est/M))/((M-1)*(N^2))
CACEvar <- (ITTY\$var*(ITTD\$est^2) + ITTD\$var*(ITTY\$est^2) -
2*Cov*ITTY\$est*ITTD\$est)/(ITTD\$est^4)
CovT <- M*sum((ITTY\$diff*ITTY\$weightsT - N*ITTY\$estT/M) *
(ITTD\$diff*ITTD\$weightsT - N*ITTD\$estT/M))/((M-1)*(N^2))
CACEvarT <- (ITTY\$varTw*(ITTD\$estT^2) + ITTD\$varTw*(ITTY\$estT^2) -
2*CovT*ITTY\$estT*ITTD\$estT)/(ITTD\$estT^4)
return(list(est = CACEest, estT = ITTY\$estT/ITTD\$estT,
var = CACEvar, varT = CACEvarT,
ITTd = ITTD, ITTy = ITTY,
cov = Cov, covT = CovT))
}
} else if (method == "standard") {
if (is.null(match)) { # without matching
Cov <- cov(ITTY\$Ysum[ITTY\$Z==0],
ITTD\$Ysum[ITTD\$Z==0])/sum(ITTY\$Z==0) +
cov(ITTY\$Ysum[ITTY\$Z==1],
ITTD\$Ysum[ITTD\$Z==1])/sum(ITTY\$Z==1)
CACEvar <- (ITTY\$var*(ITTD\$est^2) + ITTD\$var*(ITTY\$est^2) -
2*Cov*ITTY\$est*ITTD\$est)/(ITTD\$est^4)
} else { # with matching
stop("this estimator is not available for matched-pair designs.")
}
} else { # unbiased
if (is.null(match)) { # without matching
if (ITTD\$M == (ITTD\$m1*2)) {
Cov <- cov(ITTY\$Ysum[ITTY\$Z==0],
ITTD\$Ysum[ITTD\$Z==0])/sum(ITTY\$Z==0)
+ cov(ITTY\$Ysum[ITTY\$Z==1],
ITTD\$Ysum[ITTD\$Z==1])/sum(ITTY\$Z==1)
CACEvar <- 2*ITTY\$M*((var(ITTY\$Ysum[ITTY\$Z==1]) +
var(ITTY\$Ysum[ITTY\$Z==1]))*ITTD\$est^2 +
(var(ITTD\$Ysum[ITTD\$Z==1]) +
var(ITTD\$Ysum[ITTD\$Z==1]))*ITTY\$est^2 -
2*cov*ITTY\$est*ITTD\$est)/((ITTY\$N^2)*(ITTD\$est^4))
} else {
stop("the treatment and control groups must have the same number of clusters for this estimator")
}
} else { # with matching
Cov <- cov(ITTY\$diff, ITTD\$diff)
CACEvar <- (ITTY\$var*(ITTD\$est^2) +
ITTD\$var*(ITTY\$est^2) -
4*ITTY\$M*Cov*ITTY\$est*ITTD\$est/(ITTY\$N^2))/(ITTD\$est^4)
}
}
return(list(est = CACEest, var = CACEvar, ITTd = ITTD, ITTy = ITTY,
cov = Cov))
}

###
### Calculate the cluster-adjusted covariance for two sample means
### the calculation is based on the formula analogous to the one
### above.
###

covCluster <- function(Y1, Y2, grp) {

n <- c(table(grp))
ugrp <- unique(grp)
k <- length(ugrp)
N <- length(Y1)
Y1bar <- mean(Y1)
Y2bar <- mean(Y2)
Y1gbar <- Y2gbar <- rep(NA, k)
MSw <- 0
## within-group mean squares
for (i in 1:k) {
Y1gbar[i] <- mean(Y1[grp == ugrp[i]])
Y2gbar[i] <- mean(Y2[grp == ugrp[i]])
MSw <- MSw + sum((Y1[grp == ugrp[i]]-Y1gbar[i])*(Y2[grp == ugrp[i]]-Y2gbar[i]))
}
MSw <- MSw / (N-k)
## between-group mean squares
MSb <- sum(n*(Y1gbar-Y1bar)*(Y2gbar-Y2bar))/(k-1)
## intraclass correlation coefficient
n0 <- mean(n)-sum((n-mean(n))^2)/((k-1)*N)
rho <- (MSb - MSw)/(MSb + (n0-1)*MSw)