Power Calculation Example
using paramtest package
Rafael F. Bressan
2021-08-20
This is an old script I used to benchmark some power test algorithms, including the usage of the paramtest package and its parallelism capabilities. You will notice I don’t use data.table and fixest here, up to you to adapt the following code, :-)
library(dplyr)
library(sandwich)
library(paramtest)
library(microbenchmark)
library(ggplot2)Load data
data <- read.csv("Data/escolas.csv")
# Clean up the data. Never forget it!!
data <- data %>%
filter(!is.na(logico)) %>%
mutate(escola = if_else(escola == "S/escola", "99", escola)) %>%
mutate(escola = as.integer(escola))We’ll consider a fine grid of treatment effects
grid = seq(0, 0.5, 0.05)pt_func <- function(simNum, eff) {
#Drawing a sample, with replacement, of schools from our data. The idea here is that our sample well approximates the population distribution.
#and that the sampling scheme is close to random sampling from the population of interest.
sample.draw = sample(list.schools, replace = TRUE)
#Construct the artificial data from a draw. Note that repeated schools should be treated as distinct so as to mimic random sampling.
list.new.data = lapply(1:length(sample.draw), function(x){
extract = data[data$escola == sample.draw[x],]
return(extract)
})
#Concatenating data on lists
# data.artificial ends up with lower number of observations that original
# data
data.artificial = do.call(rbind, list.new.data)
#Next, we select that half of the schools will be treated
treat.draw = sample(unique(data.artificial$escola),size = length(unique(data.artificial$escola))/2, replace = F)
data.artificial$treat = 1*(data.artificial$escola %in% treat.draw)
#Create outcome
data.artificial$y = data.artificial$logico + data.artificial$treat*eff
#Running models and storing whether we reject the null that effect is 0
model1 = lm(y~treat, data = data.artificial)
se = sqrt(diag(sandwich::vcovCL(model1, cluster = data.artificial$escola)))
rej.col.1 = 1*(abs(model1$coefficients["treat"]/se["treat"]) > qnorm(0.975) )
model2 = lm(y~treat+mulher, data = data.artificial)
se = sqrt(diag(sandwich::vcovCL(model2, cluster = data.artificial$escola)))
rej.col.2 = 1*(abs(model2$coefficients["treat"]/se["treat"]) > qnorm(0.975) )
model3 = lm(y~treat+mulher+idade, data = data.artificial)
se = sqrt(diag(sandwich::vcovCL(model3, cluster = data.artificial$escola)))
rej.col.3 = 1*(abs(model3$coefficients["treat"]/se["treat"]) > qnorm(0.975) )
return(c(M1 = rej.col.1, M2 = rej.col.2, M3 = rej.col.3))
}Passing data frame along
pt_func2 <- function(simNum, df, eff) {
list.schools <- unique(df$escola)
#Drawing a sample, with replacement, of schools from our data. The idea here is that our sample well approximates the population distribution.
#and that the sampling scheme is close to random sampling from the population of interest.
sample.draw = sample(list.schools, replace = T)
#Construct the artificial data from a draw. Note that repeated schools should be treated as distinct so as to mimic random sampling.
list.new.data = lapply(1:length(sample.draw), function(x){
extract = df[df$escola == sample.draw[x],]
#extract$escola = x
return(extract)
})
#Concatenating data on lists
# data.artificial ends up with lower number of observations that original
# data
data.artificial = do.call(rbind, list.new.data)
#Next, we select that half of the schools will be treated
treat.draw = sample(unique(data.artificial$escola),
size = length(unique(data.artificial$escola))/2,
replace = F)
data.artificial$treat = 1*(data.artificial$escola %in% treat.draw)
#Create outcome
data.artificial$y = data.artificial$logico + data.artificial$treat*eff
#Running models and storing whether we reject the null that effect is 0
model1 = lm(y~treat, data = data.artificial)
se = sqrt(diag(sandwich::vcovCL(model1, cluster = data.artificial$escola)))
rej.col.1 = 1*(abs(model1$coefficients["treat"]/se["treat"]) > qnorm(0.975) )
model2 = lm(y~treat+mulher, data = data.artificial)
se = sqrt(diag(sandwich::vcovCL(model2, cluster = data.artificial$escola)))
rej.col.2 = 1*(abs(model2$coefficients["treat"]/se["treat"]) > qnorm(0.975) )
model3 = lm(y~treat+mulher+idade, data = data.artificial)
se = sqrt(diag(sandwich::vcovCL(model3, cluster = data.artificial$escola)))
rej.col.3 = 1*(abs(model3$coefficients["treat"]/se["treat"]) > qnorm(0.975) )
return(c(M1 = rej.col.1, M2 = rej.col.2, M3 = rej.col.3))
}Testing paramtest
First let as global variables: data and list.schools
list.schools <- unique(data$escola)
pt_power <- run_test(pt_func, n.iter = 1000, output = "data.frame",
eff = 0.5)## Running 1,000 tests...Test power results for one effect level only
colMeans(results(pt_power))[2:4]## M1.treat M2.treat M3.treat
## 0.663 0.664 0.720Now test passing data to pt_func2
pt_power2 <- run_test(pt_func2, n.iter = 1000, output = "data.frame",
df = data, eff = 0.5)## Running 1,000 tests...Test power results
colMeans(results(pt_power2))[2:4]## M1.treat M2.treat M3.treat
## 0.683 0.679 0.724WARNING: Benchmarking may take a long time to run depending on your number of simulations and number of replication in the benchmark (parameter “times”)
Benchmarking. Looks like there is no harm in passing the data.frame to the function.
pt_bench <- microbenchmark(
pt1 = run_test(pt_func, n.iter = 1000, output = "data.frame", eff = 0.5),
pt2 = run_test(pt_func2, n.iter = 1000, output = "data.frame", df = data,
eff = 0.5),
times = 5
)## Running 1,000 tests...
## Running 1,000 tests...
## Running 1,000 tests...
## Running 1,000 tests...
## Running 1,000 tests...
## Running 1,000 tests...
## Running 1,000 tests...
## Running 1,000 tests...
## Running 1,000 tests...
## Running 1,000 tests...autoplot(pt_bench)## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Now let’s vary the eff argument, to capture the Minimum Detectable Effect
pt_pwr2 <- grid_search(pt_func2, n.iter = 1000, output = "data.frame",
df = data, params = list(eff = grid))## Running 11,000 tests...pwr2_df <- results(pt_pwr2) %>%
group_by(eff.test) %>%
summarise(across(contains("treat"), mean))
pwr2_df## # A tibble: 11 x 4
## eff.test M1.treat M2.treat M3.treat
## <dbl> <dbl> <dbl> <dbl>
## 1 0 0.11 0.111 0.102
## 2 0.05 0.13 0.127 0.117
## 3 0.1 0.143 0.143 0.138
## 4 0.15 0.181 0.18 0.188
## 5 0.2 0.22 0.221 0.232
## 6 0.25 0.275 0.28 0.306
## 7 0.3 0.318 0.329 0.373
## 8 0.35 0.435 0.441 0.487
## 9 0.4 0.512 0.511 0.564
## 10 0.45 0.564 0.568 0.613
## 11 0.5 0.65 0.648 0.709pt_par <- grid_search(pt_func2, n.iter = 1000, output = "data.frame",
df = data, params = list(eff = grid),
parallel = "multicore", ncpus = 4)## Running 11,000 tests...par_df <- results(pt_par) %>%
group_by(eff.test) %>%
summarise(across(contains("treat"), mean))
par_df## # A tibble: 11 x 4
## eff.test M1.treat M2.treat M3.treat
## <dbl> <dbl> <dbl> <dbl>
## 1 0 0.112 0.109 0.102
## 2 0.05 0.137 0.135 0.131
## 3 0.1 0.115 0.117 0.109
## 4 0.15 0.199 0.193 0.205
## 5 0.2 0.214 0.217 0.239
## 6 0.25 0.267 0.273 0.294
## 7 0.3 0.357 0.362 0.39
## 8 0.35 0.445 0.448 0.49
## 9 0.4 0.521 0.526 0.584
## 10 0.45 0.578 0.585 0.625
## 11 0.5 0.662 0.67 0.712Is it worth to use parallelization?
Benchmark it! Parallelized simulation is almost twice as fast.
par_bench <- microbenchmark(
no_par = grid_search(pt_func2, n.iter = 1000, output = "data.frame",
df = data, params = list(eff = grid)),
par = grid_search(pt_func2, n.iter = 1000, output = "data.frame",
df = data, params = list(eff = grid),
parallel = "multicore", ncpus = 4),
times = 5
)## Running 11,000 tests...
## Running 11,000 tests...
## Running 11,000 tests...
## Running 11,000 tests...
## Running 11,000 tests...
## Running 11,000 tests...
## Running 11,000 tests...
## Running 11,000 tests...
## Running 11,000 tests...
## Running 11,000 tests...par_bench## Unit: seconds
## expr min lq mean median uq max neval
## no_par 174.20182 174.2672 187.46722 174.35237 174.42401 240.0907 5
## par 86.00311 86.5560 92.83074 86.60456 87.20923 117.7808 5