R example 4: the Children’s television workshop experiment.

The Children’s television workshop experiment.

## Data input

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls

W <- rep(c(0, 1), 8)
G <- as.vector(rbind(1:8, 1:8))
X <- c(12.9, 12, 15.1, 12.3, 16.8, 17.2, 15.8, 18.9, 13.9, 15.3, 14.5, 16.6, 17, 16, 15.8, 20.1)
Y <- c(54.6, 60.6, 56.5, 55.5, 75.2, 84.8, 75.6, 101.9, 55.3, 70.6, 59.3, 78.4, 87, 84.2, 73.7, 108.6)
data <- data.frame(G, W, X, Y)
data

library(ggpubr)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

temp <- data.frame("Control" = data$Y[data$W == 0], "Treatment" = data$Y[data$W == 1])
temp

ggpaired(data = temp, cond1 = "Control", cond2 = "Treatment", fill = "condition", palette = "jco", line.color = "gray", 
         ylab = "Post-test Score")

`gather_()` was deprecated in tidyr 1.2.0.

temp <- data.frame("Control" = data$X[data$W == 0], "Treatment" = data$X[data$W == 1])
ggpaired(data = temp, cond1 = "Control", cond2 = "Treatment", fill = "condition", palette = "jco", line.color = "gray", 
         ylab = "Pre-test Score")

• Compute Fisher’s exact p-value

testStat <- function(Y0, Y1, w, pair_idx, method = "mean diff") {
 n <- length(Y0)
 if (method == "mean diff")
   val <- abs(mean(Y0[w == 1]) - mean(Y0[w == 0])) ## replace abs(mean(Y1[w == 1]) - mean(Y0[w == 0])-tau)
 else if (method == "rank diff"){
  r <- rank(Y0)
  val <- abs(mean(r[w == 1]) - mean(r[w == 0]))
 } else if (method == "within-pair rank diff") {
   pairs <- unique(pair_idx)
   v <- sapply(pairs, function(idx) {
     Yt <- Y0[pair_idx == idx & w == 1]
     Yc <- Y0[pair_idx == idx & w == 0]
     return((Yt > Yc) - (Yt < Yc))
   })
   val <- abs(2 * sum(v)/length(Y0))
 }
  return(val)
}

find.p.value <- function(tau, yobs, w, pair_idx, seed = 0, method = "mean diff", 
                         B = 10000, error = 1e-10){
  ## error is to fix our problems with floating point errors
  n <- length(yobs)
  y0 <- rep(NA, n)
  y1 <- rep(NA, n)
  
  y1[w == 1] <- yobs[w == 1]
  y0[w == 0] <- yobs[w == 0]
  
  y1[w == 0] <- yobs[w == 0] + tau
  y0[w == 1] <- yobs[w == 1] - tau
  
  t.obs <-  testStat(y0, y1, w, pair_idx, method = method)
  
  
  set.seed(seed)
  pairs <- unique(pair_idx)
  t.stat <- sapply(1:B, function(i){
    w.iter <- rep(NA, n)
    for (i in 1:length(pairs)) {
      temp.idx <- which(pair_idx == pairs[i])
      w.iter[sample(temp.idx)] <- c(0, 1)
    }
    return(testStat(y0, y1, w.iter, pair_idx, method = method))
  })

  # subtracting error fixes floating point errors
  p.val <- sum(t.stat >= (t.obs - error))/B
  res <- c(t.obs, p.val)
  names(res) <- paste(method, c("test stats", "p-value"))
  return(res)
}


find.p.value(0, data$Y, data$W, data$G)

mean diff test stats    mean diff p-value 
             13.4250               0.0332 

find.p.value(0, data$Y, data$W, data$G, method = "rank diff")

rank diff test stats    rank diff p-value 
              3.7500               0.0332 

find.p.value(0, data$Y, data$W, data$G, method = "within-pair rank diff")

within-pair rank diff test stats    within-pair rank diff p-value 
                             0.5                              0.3 

• What if we use ignore the pairwise randomization mechanism and treat it as complete randomization?

find.p.value.wrong <- function(tau, yobs, w, pair_idx, seed = 0, method = "mean diff", 
                         B = 10000, error = 1e-10){
  ## error is to fix our problems with floating point errors
  n <- length(yobs)
  y0 <- rep(NA, n)
  y1 <- rep(NA, n)
  
  y1[w == 1] <- yobs[w == 1]
  y0[w == 0] <- yobs[w == 0]
  
  y1[w == 0] <- yobs[w == 0] + tau
  y0[w == 1] <- yobs[w == 1] - tau
  
  t.obs <-  testStat(y0, y1, w, pair_idx, method)
  
  
  set.seed(seed)
  pairs <- unique(pair_idx)
  t.stat <- sapply(1:B, function(i){
    w.iter <- sample(w)
    return(testStat(y0, y1, w.iter, pair_idx, method = method))
  })

  # subtracting error fixes floating point errors
  p.val <- sum(t.stat >= (t.obs - error))/B
  res <- c(t.obs, p.val)
  names(res) <- paste(method, c("test stats", "p-value"))
  return(res)
}


find.p.value.wrong(0, data$Y, data$W, data$G)

mean diff test stats    mean diff p-value 
             13.4250               0.1058 

find.p.value.wrong(0, data$Y, data$W, data$G, method = "rank diff")

rank diff test stats    rank diff p-value 
              3.7500               0.1282 

We get larger p-values from complete randomization mechanism.

• Compare 95% CI of Fisher V.S. Neyman

## CI from Fisher
pvals <- sapply(0:30, function(tau) find.p.value(tau, data$Y, data$W, data$G)[2])
names(pvals) <- 0:30
pvals

     0      1      2      3      4      5      6      7      8      9     10     11     12 
0.0332 0.0423 0.0568 0.0740 0.0921 0.1251 0.1618 0.1993 0.3063 0.3722 0.4869 0.6235 0.7679 
    13     14     15     16     17     18     19     20     21     22     23     24     25 
0.9322 0.8945 0.7365 0.5933 0.4632 0.3318 0.2616 0.2060 0.1564 0.1232 0.0909 0.0649 0.0649 
    26     27     28     29     30 
0.0469 0.0329 0.0245 0.0245 0.0245 

LB <- sapply(seq(1, 2, 0.1), function(tau) find.p.value(tau, data$Y, data$W, data$G)[2])
names(LB) <- seq(1, 2, 0.1)
LB

     1    1.1    1.2    1.3    1.4    1.5    1.6    1.7    1.8    1.9      2 
0.0423 0.0423 0.0423 0.0423 0.0423 0.0423 0.0487 0.0487 0.0487 0.0487 0.0568 

UB <- sapply(seq(25, 26, 0.1), function(tau) find.p.value(tau, data$Y, data$W, data$G)[2])
names(UB) <- seq(25, 26, 0.1)
UB

    25   25.1   25.2   25.3   25.4   25.5   25.6   25.7   25.8   25.9     26 
0.0649 0.0649 0.0547 0.0547 0.0547 0.0547 0.0469 0.0469 0.0469 0.0469 0.0469 

The 95% CI from Fisher’s exact p-values is [1.9, 25.6]

## CI from Neyman
pairs <- unique(data$G)
tau_hat_vec <- sapply(pairs, function(idx) {
  Yt <- data$Y[data$G == idx & data$W == 1]
  Yc <- data$Y[data$G == idx & data$W == 0]
  return(Yt - Yc)
})

tau_hat <- mean(tau_hat_vec)
sd_tau_hat <- sd(tau_hat_vec)/sqrt(length(tau_hat_vec))
  
## 95% CI
print(c(tau_hat- 1.96 * sd_tau_hat, tau_hat + 1.96 * sd_tau_hat))

[1]  4.337779 22.512221

• What does Neyman’s CI look like if we treat the data as from completely randomized experiment

## CI from Neyman

tau_hat <- mean(data$Y[data$W == 1]) - mean(data$Y[data$W == 0])

st <- sd(data$Y[data$W == 1])
sc <- sd(data$Y[data$W == 0])
varhat_tauhat <- st^2/sum(data$W == 1) + sc^2/sum(data$W == 0)

print(c(tau_hat - 1.96 * sqrt(varhat_tauhat), tau_hat + 1.96 * sqrt(varhat_tauhat)))

[1] -1.957377 28.807377

• Linear regression results

## need to create pair-level data
pairs <- unique(data$G)
X_bar <- mean(data$X)
group_data <- sapply(pairs, function(idx) {
  Yt <- data$Y[data$G == idx & data$W == 1]
  Yc <- data$Y[data$G == idx & data$W == 0]
  Xt <- data$X[data$G == idx & data$W == 1]
  Xc <- data$X[data$G == idx & data$W == 0]
  Y_diff <- Yt - Yc
  X_diff <- Xt - Xc
  X_mean <- (Xt + Xc)/2 - X_bar
  return(c(X_mean, X_diff, Y_diff))
})
group_data

        [,1]    [,2]   [,3]    [,4]    [,5]    [,6]    [,7]    [,8]
[1,] -3.1875 -1.9375 1.3625  1.7125 -1.0375 -0.0875  0.8625  2.3125
[2,] -0.9000 -2.8000 0.4000  3.1000  1.4000  2.1000 -1.0000  4.3000
[3,]  6.0000 -1.0000 9.6000 26.3000 15.3000 19.1000 -2.8000 34.9000

group_data <- data.frame(t(group_data))
colnames(group_data) <- c("X_mean", "X_diff", "Y_diff")
group_data

lm_result1 <- lm(Y_diff ~ X_diff, data =group_data)
summary(lm_result1) ## substantially smaller sd of the estimate

Call:
lm(formula = Y_diff ~ X_diff, data = group_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4238 -1.2954 -0.2577  2.0825  5.0432 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.9943     1.4096   6.381 0.000697 ***
X_diff        5.3705     0.5992   8.964 0.000108 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.734 on 6 degrees of freedom
Multiple R-squared:  0.9305,    Adjusted R-squared:  0.9189 
F-statistic: 80.34 on 1 and 6 DF,  p-value: 0.0001077

lm_result2 <- lm(Y_diff ~ X_diff + X_mean, data =group_data)
summary(lm_result2) 

Call:
lm(formula = Y_diff ~ X_diff + X_mean, data = group_data)

Residuals:
      1       2       3       4       5       6       7       8 
-0.4555  5.1236  0.1038  1.1179 -2.6106 -1.9925 -4.4907  3.2041 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.5220     1.4745   5.780 0.002181 ** 
X_diff        5.9431     0.8138   7.303 0.000754 ***
X_mean       -1.0298     0.9969  -1.033 0.348978    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.713 on 5 degrees of freedom
Multiple R-squared:  0.9427,    Adjusted R-squared:  0.9198 
F-statistic: 41.15 on 2 and 5 DF,  p-value: 0.0007848