-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path2_run_sims.R
154 lines (146 loc) · 6.85 KB
/
2_run_sims.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
library(tidyverse)
library(ltmle)
library(future)
library(doFuture)
library(foreach)
library(sandwich)
library(lme4)
library(gee)
source("0_DGP_and_helpers.R")
source("1_fit_models.R")
one_set <- function(specification, nsims = 4, cores = 8, trial = T, seed = NA,
tmle_misspec_test = F,
estimate_pscore = T){
sizes <- specification[[1,"sizes"]]
sizes <- sizes[[1]]
#print(sizes)
if(!is.na(seed)){
set.seed(seed, kind = "L'Ecuyer-CMRG")
}
if(specification[[1,"link_function"]] == "identity"){
PATE <- cbind.data.frame(target_parameter = specification[[1,"t_eff"]])
} else { # Get true PATE values by averaging...
pop <- bind_rows(lapply(as.list(1:500), function(x) gen_obs_dat(LPS_only = T,
trial = trial,
sizes = sizes,#specification[[1,"sizes"]],
t_eff = specification[[1,"t_eff"]],
#control_arm_no_ranef = specification[[1,"control_arm_no_ranef"]],
ranef_sd = specification[[1,"ranef_sd"]],
link_function = specification[[1,"link_function"]],
sim_type = specification[[1,"sim_type"]])))
UY <- runif(nrow(pop))
P0 <- mean(as.numeric(plogis(pop$LP0) > UY))
P1 <- mean(as.numeric(plogis(pop$LP1) > UY))
PATE <- cbind.data.frame(target_parameter = (P1/(1-P1)) / (P0/(1-P0)))
}
#print(sizes)
#print(specification)
plan(multisession, workers = cores)
if(!is.na(seed)){
set.seed(seed, kind = "L'Ecuyer-CMRG")
}
results <- data.frame(foreach(j = 1:nsims, .combine = rbind,
.options.future = list(seed = T, globals = structure(T)),
.errorhandling = "remove") %dofuture% {
fit_models(dat = gen_obs_dat(sizes = sizes,
trial = trial,
#control_arm_no_ranef = sizes[[1,"control_arm_no_ranef"]],
t_eff = specification[[1,"t_eff"]],
ranef_sd = specification[[1,"ranef_sd"]],
link_function = specification[[1,"link_function"]],
sim_type = specification[[1,"sim_type"]]),
sizes = sizes,
tmle_misspec_test = tmle_misspec_test,
#control_arm_irgtt = sizes[[1]][["control_arm_irgtt"]],
link_function = specification[[1,"link_function"]],
#v = specification[[1,"v"]],
sim_type = specification[[1,"sim_type"]],
estimate_pscore = estimate_pscore)
})
plan(sequential)
return(cbind.data.frame(nsims = nsims,
specification,
control_arm_irgtt = sizes[["control_arm_irgtt"]],
control_arm_no_ranef = sizes[["control_arm_no_ranef"]],
control_arm_clusters = sizes[["clusters_per_arm_c"]],
treatment_arm_clusters = sizes[["clusters_per_arm_t"]],
control_arm_cluster_size = sizes[["cluster_size_c"]],
treatment_arm_cluster_size = sizes[["cluster_size_t"]],
results, PATE))
}
run_full_set <- function(sim_type = c("complex_covars", "main_terms", "overfit"),
link_function = "identity",
t_eff = c(.25, 0),
ranef_sd = c(icc_to_sd(.001),
icc_to_sd(.01),
icc_to_sd(.05),
icc_to_sd(.10)),
#v = 10,
tmle_misspec_test = F,
estimate_pscore = T, nsims = 4, cores = 8, seed = NA, trial = T,
sizes = list(c(cluster_size_t = 2, clusters_per_arm_t = 10,
cluster_size_c = 2, clusters_per_arm_c = 10,
control_arm_irgtt = T, control_arm_no_ranef = T))){
specification <- expand_grid(sizes,
sim_type,
link_function,
t_eff,
ranef_sd,
#v,
seed)
#print(specification)
#return(specification)
return(data.frame(foreach(j = 1:nrow(specification), .combine = rbind) %do%
{one_set(specification = specification[j,],
seed = seed,
trial = trial,
cores = cores,
nsims = nsims,
tmle_misspec_test = tmle_misspec_test,
estimate_pscore = estimate_pscore)}))
}
# run_full_set(link_function = "logit", tmle_misspec_test = T)
# sim_type = "nonadditive"
# link_function = "identity"
# t_eff = c(.25, 0)
# ranef_sd = c(.05)
# control_arm_no_ranef = T
# v = 10
# control_arm_irgtt = T
# estimate_pscore = T
# nsims = 2
# cores = 10
# seed = NA
#
# sizes <- list(c(cluster_size_t = 5, clusters_per_arm_t = 10,
# cluster_size_c = 5, clusters_per_arm_c = 10,
# control_arm_irgtt = T, control_arm_no_ranef = T))
#
# specification <- expand_grid(sizes,
# sim_type,
# link_function,
# t_eff,
# ranef_sd,
# control_arm_no_ranef,
# v,
# control_arm_irgtt,
# seed)
# specification[1,]
#
# one_set(specification = specification[1,],
# cores = cores,
# nsims = nsims,
# estimate_pscore = estimate_pscore)
#
#
# fit_models(dat = gen_obs_dat(seed = specification[[1,"seed"]],
# sizes = specification[[1,"sizes"]],
# control_arm_no_ranef = specification[[1,"control_arm_no_ranef"]],
# t_eff = specification[[1,"t_eff"]],
# ranef_sd = specification[[1,"ranef_sd"]],
# link_function = specification[[1,"link_function"]],
# sim_type = specification[[1,"sim_type"]]),
# control_arm_irgtt = specification[[1,"control_arm_irgtt"]],
# link_function = specification[[1,"link_function"]],
# v = specification[[1,"v"]],
# estimate_pscore = T)