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)))