Example 1

D <- c(1, 1, 0, 0)
Y <- c(40, 20, 10, 15)
cbind(D, Y)
##      D  Y
## [1,] 1 40
## [2,] 1 20
## [3,] 0 10
## [4,] 0 15
obs_diff <-  mean(Y[D==1]) - mean(Y[D==0])
obs_diff
## [1] 17.5
n <- 4 # number of obs
treatment <- c(0,1) # types of treatments

# list up all possible treatment status profile
randomization_profile <- expand.grid(rep(list(treatment),n)) 
randomization_profile
##    Var1 Var2 Var3 Var4
## 1     0    0    0    0
## 2     1    0    0    0
## 3     0    1    0    0
## 4     1    1    0    0
## 5     0    0    1    0
## 6     1    0    1    0
## 7     0    1    1    0
## 8     1    1    1    0
## 9     0    0    0    1
## 10    1    0    0    1
## 11    0    1    0    1
## 12    1    1    0    1
## 13    0    0    1    1
## 14    1    0    1    1
## 15    0    1    1    1
## 16    1    1    1    1
# select only where sum is equal to 2
randomization_profile <- randomization_profile[rowSums(randomization_profile) == 2,] 
randomization_profile
##    Var1 Var2 Var3 Var4
## 4     1    1    0    0
## 6     1    0    1    0
## 7     0    1    1    0
## 10    1    0    0    1
## 11    0    1    0    1
## 13    0    0    1    1
# permute
permute_tau <- c() 
for(i in 1:nrow(randomization_profile)){
  D_star <- unlist(randomization_profile[i,])
  permute_tau[i] <- mean(Y[D_star == 1]) - mean(Y[D_star == 0])
}
permute_tau
## [1]  17.5   7.5 -12.5  12.5  -7.5 -17.5
#Two sided test
pval <- length(permute_tau[abs(permute_tau) >= abs(obs_diff)]) / nrow(randomization_profile)


density_perm <- density(permute_tau)
plot(density_perm,
     main = 'All Possile Effects',
     xlab = 'Treatment Effects')

#Row names of levels indicating areas of coloring
x1 <- min(which(density_perm$x >= -abs(obs_diff)))
x2 <- max(which(density_perm$x <=  min(density_perm$x)))
x3 <- min(which(density_perm$x >= abs(obs_diff)))
x4 <- max(which(density_perm$x <  max(density_perm$x)))

with(density_perm, polygon(x = c(x[c(x1,x1:x2,x2)]),
                           y = c(0, y[x1:x2], 0),
                           col="gray"))
with(density_perm, polygon(x = c(x[c(x3,x3:x4,x4)]),
                           y = c(0, y[x3:x4], 0),
                           col="gray"))
abline(v = obs_diff, col = 'red', lty = 2)
text(x = -10,
     y = 0.01,
     labels = paste('Pvalue:', round(pval, 2), '\n','Observed tau = ', round(obs_diff, 2)))

Example 2

set.seed(21930)

D <- c(rep(1, 5), rep(0, 15))
Y <- rnorm(20, 10, 5) + rnorm(20, 7, 1) * D
cbind(D, Y)
##       D         Y
##  [1,] 1 11.809519
##  [2,] 1 19.018942
##  [3,] 1 18.123395
##  [4,] 1 22.474232
##  [5,] 1 14.985495
##  [6,] 0  7.460489
##  [7,] 0 14.568118
##  [8,] 0 10.939037
##  [9,] 0  6.143088
## [10,] 0 16.505258
## [11,] 0  3.092529
## [12,] 0 16.387006
## [13,] 0 13.646043
## [14,] 0  9.512907
## [15,] 0 -1.438535
## [16,] 0  9.011769
## [17,] 0  3.493713
## [18,] 0 13.173865
## [19,] 0  8.739871
## [20,] 0 16.252682
obs_diff <-  mean(Y[D==1]) - mean(Y[D==0])
obs_diff
## [1] 7.449794
n <- 20 # number of obs
treatment <- c(0,1) # types of treatments

# list up all possible treatment status profile
randomization_profile <- expand.grid(rep(list(treatment), n)) 
# select only where sum is equal to 2
randomization_profile <- randomization_profile[rowSums(randomization_profile) == 5,] 
# permute
permute_tau <- c() 
for(i in 1:nrow(randomization_profile)){
  D_star <- unlist(randomization_profile[i,])
  permute_tau[i] <- mean(Y[D_star == 1]) - mean(Y[D_star == 0])
}
#Two sided test
pval <- length(permute_tau[abs(permute_tau) >= abs(obs_diff)]) / nrow(randomization_profile)
pval
## [1] 0.0119324
density_perm <- density(permute_tau)
plot(density_perm,
     main = 'All Possile Effects',
     xlab = 'Treatment Effects')

#Row names of levels indicating areas of coloring
x1 <- min(which(density_perm$x >= -abs(obs_diff)))
x2 <- max(which(density_perm$x <=  min(density_perm$x)))
x3 <- min(which(density_perm$x >= abs(obs_diff)))
x4 <- max(which(density_perm$x <  max(density_perm$x)))

with(density_perm, polygon(x = c(x[c(x1,x1:x2,x2)]),
                           y = c(0, y[x1:x2], 0),
                           col="gray"))
with(density_perm, polygon(x = c(x[c(x3,x3:x4,x4)]),
                           y = c(0, y[x3:x4], 0),
                           col="gray"))
abline(v = obs_diff, col = 'red', lty = 2)
text(x = -10,
     y = 0.01,
     labels = paste('Pvalue:', round(pval, 2), '\n','Observed tau = ', round(obs_diff, 2)))

  1. Graduate School of Economics, Waseda University, ↩