### Section 5.1 - Single-stage problem ### R code from vignette source "example1.Rnw" ################################################### ### code chunk number 1: preliminaries ################################################### library("knitr") options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE) opts_chunk$set(fig.path="figure/fig_") options(tikzDefaultEngine = "xetex", tikzMetricPackages = c("\\usepackage[utf8{inputenc}", "\\usepackage[T1]{fontenc}", "\\usetikzlibrary{calc}", "\\usepackage{amssymb}")) library("polle") library("mets") library("SuperLearner") library("lava") library("ggplot2") library("data.table") ################################################### ### code chunk number 2: ex1_sim ################################################### library("polle") par0 <- c(p = .3, k = .1, d = .5, a = 1, b = -2.5, c = 3) d <- sim_single_stage(n = 5e2, seed = 1, par = par0) head(d, 4) ################################################### ### code chunk number 3: ex1_data ################################################### pd <- policy_data(d, action = "A", covariates = list("Z", "B", "L"), utility = "U") pd ################################################### ### code chunk number 4: ex1_static1 ################################################### p1 <- policy_def(1, name = "(A=1)") p1 ################################################### ### code chunk number 5: ex1_pol_eval ################################################### p1(pd) |> head(3) ################################################### ### code chunk number 6: ex1_eval1 ################################################### (pe1 <- policy_eval(pd, policy = p1)) ################################################### ### code chunk number 7: ex1_static3 ################################################### (pe0 <- policy_eval(pd, policy = policy_def(0, name = "(A=0)"))) ################################################### ### code chunk number 8: ex1_estimateATE ################################################### estimate(merge(pe0, pe1), function(x) x[2] - x[1], labels = "ATE") ################################################### ### code chunk number 9: ex1_merge ################################################### IC(merge(pe0, pe1)) |> head(3) ################################################### ### code chunk number 10: ex1_opt_pol ################################################### p_opt <- policy_def( function(Z, L) 1 * ((par0["c"] * Z + par0["a"] * L + par0["b"]) > 0), name = "optimal") ################################################### ### code chunk number 11: ex1_opt_pol_value ################################################### set.seed(1) policy_eval(pd, policy = p_opt, g_models = g_sl(SL.library = c("SL.glm", "SL.ranger", "SL.gam")), q_models = q_sl(SL.library = c("SL.glm", "SL.ranger", "SL.gam")), M = 5) ################################################### ### code chunk number 12: ex1_drql ################################################### pl <- policy_learn(type = "blip", L = 5, control = control_blip(blip_models = q_glm(formula = ~ Z + L))) ################################################### ### code chunk number 13: ex1_drql_pol ################################################### set.seed(1) po <- pl(pd, g_models = g_sl(SL.library = c("SL.glm", "SL.ranger", "SL.gam")), q_models = q_sl(SL.library = c("SL.glm", "SL.ranger", "SL.gam"))) po ################################################### ### code chunk number 14: ex1_drql_d ################################################### get_policy(po)(pd) |> head(3) ################################################### ### code chunk number 15: getpolfunex ################################################### pf <- get_policy_functions(po, stage = 1) pf(H = data.table(Z = c(1, -2), L = c(0.5, -1))) ################################################### ### code chunk number 16: ex1_drql_eval ################################################### set.seed(1) pe <- policy_eval(pd, policy_learn = pl, g_models = g_sl(SL.library = c("SL.glm", "SL.ranger", "SL.gam")), q_models = q_sl(SL.library = c("SL.glm", "SL.ranger", "SL.gam")), M = 5) pe ################################################### ### code chunk number 17: ex1_plot ################################################### plot_data <- get_history(pd)$H plot_data <- merge(plot_data, pe$policy_actions) ggplot(plot_data) + geom_point(aes(x = Z, y = L, color = d, shape = d)) + geom_abline(intercept = -par0["b"] / par0["a"], slope = -par0["c"] / par0["a"]) + theme_bw() ### Section 5.2 - Two-stage problem ### R code from vignette source "example2.Rnw" ################################################### ### code chunk number 1: preliminaries ################################################### library("knitr") options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE) opts_chunk$set(fig.path="figure/fig_") options(tikzDefaultEngine = "xetex", tikzMetricPackages = c("\\usepackage[utf8{inputenc}", "\\usepackage[T1]{fontenc}", "\\usetikzlibrary{calc}", "\\usepackage{amssymb}")) library("polle") library("lava") library("data.table") library("mets") library("SuperLearner") library("ggplot2") library("gridExtra") library("policytree") ################################################### ### code chunk number 2: ex2_sim ################################################### par0 <- c(gamma = 0.5, beta = 1) d <- sim_two_stage(2e3, seed = 1, par = par0) head(d[, -(1:2)], 3) ################################################### ### code chunk number 3: ex2_pd ################################################### pd <- policy_data(d, action = c("A_1", "A_2"), covariates = list(L = c("L_1", "L_2"), C = c("C_1", "C_2")), utility = c("U_1", "U_2", "U_3")) pd ################################################### ### code chunk number 4: ex2_his ################################################### get_history(pd, stage = 1, full_history = TRUE)$H |> head(3) get_history(pd, stage = 2, full_history = TRUE)$H |> head(3) ################################################### ### code chunk number 5: ex2_p_opt ################################################### kappa <- function(mu){ pnorm(q = -mu, lower.tail = FALSE) * (mu + dnorm(-mu) / pnorm(-mu, lower.tail = FALSE)) } p_opt <- policy_def( list(function(C_1, L_1){ 1 * ((C_1 + kappa(par0[["gamma"]] * L_1 + 1) - kappa(par0[["gamma"]] * L_1)) > 0) }, function(C_2){ 1 * (C_2 > 0) }), full_history = TRUE, name = "optimal") p_opt ################################################### ### code chunk number 6: ex2_p_opt_head ################################################### p_opt(pd) |> head(3) ################################################### ### code chunk number 7: ex2_p_opt_head_full ################################################### get_history(pd, full_history = FALSE)$H |> head(3) get_history(pd, full_history = FALSE)$A |> head(3) ################################################### ### code chunk number 8: ex2_p_opt_eval ################################################### pe_opt <- policy_eval(pd, policy = p_opt, g_models = g_glm(), g_full_history = FALSE, q_models = list(q_glm(), q_glm()), q_full_history = TRUE) pe_opt ################################################### ### code chunk number 9: ex2_p_opt_g ################################################### get_g_functions(pe_opt) ################################################### ### code chunk number 10: ex2_p_opt_q ################################################### get_q_functions(pe_opt) ################################################### ### code chunk number 11: ex2_predict ################################################### predict(get_g_functions(pe_opt), pd) |> head(3) predict(get_q_functions(pe_opt), pd) |> head(3) ################################################### ### code chunk number 12: ex2_ptl ################################################### pl <- policy_learn(type = "ptl", control = control_ptl(policy_vars = c("C", "L")), full_history = FALSE, L = 5) ################################################### ### code chunk number 13: ex2_ptl_object ################################################### po <- pl(pd, g_models = g_glm(), g_full_history = FALSE, q_models = q_glm(), q_full_history = TRUE) get_policy(po)(pd) |> head(3) ################################################### ### code chunk number 14: ex2_ptl_eval ################################################### set.seed(1) pe <- policy_eval(pd, policy_learn = pl, g_models = g_glm(), g_full_history = FALSE, q_models = q_glm()) pe ################################################### ### code chunk number 15: ex2_ptl_insp ################################################### po <- get_policy_object(pe) po$ptl_objects ################################################### ### code chunk number 16: ex2_plot_data ################################################### library("policytree") plot_data <- merge(get_history(pd)$H, get_policy(po)(pd)) plot_data$d <- as.factor(plot_data$d) plot_data_2 <- plot_data[stage == 2, ] plot_data_2$`d: stage 2` <- plot_data_2$d plot_data_1 <- plot_data[stage == 1, ] plot_data_1$`d: stage 1` <- plot_data_1$d ################################################### ### code chunk number 17: ex2_plot2 ################################################### plot2 <- ggplot(plot_data_2) + geom_point(aes(x = L, y = C, shape = `d: stage 2`, color = `d: stage 2`)) + geom_hline(yintercept = 0, linewidth = 1) + theme_bw() ################################################### ### code chunk number 18: ex2_plot1 ################################################### fun <- function(l){ -(kappa(par0[["gamma"]] * l + 1) - kappa(par0[["gamma"]] * l)) } plot1 <- ggplot(plot_data_1) + geom_point(aes(x = L, y = C, shape = `d: stage 1`, color = `d: stage 1`)) + geom_function(fun = fun, size = 1) + theme_bw() ################################################### ### code chunk number 19: ex2_plot ################################################### grid.arrange(plot1, plot2) ### Section 6 - K-2 literacy intervention ### R code from vignette source "harvard_example.Rnw" ################################################### ### code chunk number 1: preliminaries ################################################### library("knitr") options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE) opts_chunk$set(fig.path="figure/fig_") options(tikzDefaultEngine = "xetex", tikzMetricPackages = c("\\usepackage[utf8{inputenc}", "\\usepackage[T1]{fontenc}", "\\usetikzlibrary{calc}", "\\usepackage{amssymb}")) library("polle") library("future.apply") ################################################### ### code chunk number 2: harvard_data_load ################################################### library("readstata13") d <- readstata13::read.dta13("k2smart_public.dta") d <- transform(d, utility = as.numeric(fa_maprit) - as.numeric(sp_maprit), cct = as.logical(cct), responder = as.logical(responder), text = as.logical(text), maprit = as.numeric(sp_maprit), teacher = as.character(t_id_public), attend = as.numeric(sp_pctattend_yr), trc = as.numeric(sp_trcbook_num), dib = as.numeric(sp_dib_score), ell = as.logical(ell), iep = as.logical(iep), grade = as.character(grade_final), male = as.logical(male), familynight = as.logical(familynight) ) ################################################### ### code chunk number 3: harvard_subset ################################################### d <- subset(d, !is.na(utility)) ################################################### ### code chunk number 4: harvard_transform ################################################### d <- transform(d, A_1 = ifelse(cct, "cct", "lt"), A_2 = ifelse(responder, "continue", ifelse(text, "text", "notext")) ) ################################################### ### code chunk number 5: harvard_pd ################################################### pd <- policy_data(d, action = c("A_1", "A_2"), utility = "utility", baseline=c("maprit", "male", "ell", "iep", "attend", "trc", "dib", "grade", "familynight"), covariates = list(responder = c(NA, "responder"))) print(pd) ################################################### ### code chunk number 6: harvard_stat_pol ################################################### actions <- list(c("cct", "text"), c("cct", "notext"), c("lt", "text"), c("lt", "notext")) static_policies <- lapply(actions, function(a){ policy_def(list( function(...) a[1], function(responder) ifelse(responder, "continue", a[2])), name = paste(a, collapse = "_")) } ) head(static_policies[[1]](pd), 4) ################################################### ### code chunk number 7: harvard_pe_stat_pol ################################################### gm <- list(g_empir(~ 1), g_empir(~ responder)) pe_static_policies_ipw <- lapply(static_policies, function(p){ policy_eval(pd, policy = p, g_models = gm, type = "ipw") } ) do.call("merge", pe_static_policies_ipw) ################################################### ### code chunk number 8: harvard_print_g_stat ################################################### print(get_g_functions(pe_static_policies_ipw[[1]])) predict(get_g_functions(pe_static_policies_ipw[[1]]), pd)[c(1, 2, 43, 44), ] ################################################### ### code chunk number 9: harvard_get_his_names ################################################### get_history_names(pd, stage = 1) get_history_names(pd, stage = 2) ################################################### ### code chunk number 10: havard_q_sl ################################################### sl_lib <- c("SL.mean", "SL.glm", "SL.gam", "SL.ranger", "SL.nnet") qm <- list( q_sl(formula = ~ . - responder_1, SL.library = sl_lib), q_sl(formula = ~ . - responder_1 - responder_2, SL.library = sl_lib) ) ################################################### ### code chunk number 11: harvard_static_dr ################################################### library("future.apply") plan(list(tweak("multisession", workers = 4))) set.seed(1) pe_static_policies_dr <- lapply(static_policies, function(p){ set.seed(1) policy_eval(pd, policy = p, g_models = gm, q_models = qm, q_full_history = TRUE, type = "dr", M = 25) } ) plan("sequential") print(do.call("merge", pe_static_policies_dr)) ################################################### ### code chunk number 12: havard_est_id_print ################################################### library("lava") (est <- estimate(do.call("merge", pe_static_policies_dr), id = d$teacher)) ################################################### ### code chunk number 13: harvard_example.Rnw:265-267 ################################################### pdiff <- function(n) lava::contr(lapply(seq(n - 1), \(x) seq(x, n))) estimate(est, f = pdiff(4)) ################################################### ### code chunk number 14: harvard_example.Rnw:275-279 ################################################### estimate(conditional(pe_static_policies_dr[[1]], pd, "male"), id = d$teacher ) ################################################### ### code chunk number 15: harvard_example.Rnw:290-303 ################################################### qvm_formulas <- list( qvm_1 = list(~ 1, ~ A_1), qvm_2 = list(~ maprit, ~ A_1 + maprit), qvm_3 = list(~ male, ~ A_1 + male), qvm_4 = list(~ grade, ~ A_1 + grade), qvm_5 = list(~ male + maprit, ~ A_1 + male + maprit), qvm_6 = list(~ male + grade + maprit, ~ A_1 + grade + male + maprit)) qvm <- lapply(qvm_formulas, function(form){ list(q_glm(form[[1]]), q_glm(form[[2]])) } ) ################################################### ### code chunk number 16: harvard_pl_drql ################################################### pl_drql <- mapply(qvm, names(qvm), FUN = function(qv, name){ policy_learn(type = "drql", control = control_drql(qv_models = qv), full_history = TRUE, alpha = 0.01, L = 25, cross_fit_g_models = FALSE, name = name) } ) ################################################### ### code chunk number 17: harvard_pe_drql ################################################### plan(list(tweak("multisession", workers = 4))) set.seed(1) pe_drql <- lapply(pl_drql, function(pl){ set.seed(1) policy_eval(pd, policy_learn = pl, g_models = gm, q_models = qm, q_full_history = TRUE, type = "dr", M = 25) } ) plan("sequential") estimate(do.call(what = "merge", unname(pe_drql)), id = d$teacher) ################################################### ### code chunk number 18: harvard_qvm_male ################################################### set.seed(1) po_drql_male <- pl_drql[["qvm_3"]](pd, g_models = gm, q_models = qm, q_full_history = TRUE) pa_drql_male <- get_policy(po_drql_male)(pd) head(pa_drql_male, 4) ################################################### ### code chunk number 19: harvard_qvm_male_2 ################################################### pa_drql_male <- merge(pa_drql_male, get_history(pd)$H) pa_drql_male[, .N, list(stage, male, d)][order(stage, male, d)] ################################################### ### code chunk number 20: harvard_pol_6 ################################################### po_drql_6 <- pl_drql[["qvm_6"]](pd, g_models = gm, q_models = qm, q_full_history = TRUE) ################################################### ### code chunk number 21: harvard_plot ################################################### plot_data <- get_history(pd)$H plot_data <- merge(plot_data, get_policy(po_drql_6)(pd), by = c("id", "stage")) library("ggplot2") ggplot(plot_data) + geom_point(aes(x = grade, y = maprit, color = d)) + facet_wrap(~ stage + male, labeller = "label_both") + theme_bw() ################################################### ### code chunk number 22: harvard_qvm_6 ################################################### plan(list(tweak("multisession", workers = 4))) set.seed(1) pe_plugin <- policy_eval(pd, policy = get_policy(po_drql_6), g_models = gm, q_models = qm, q_full_history = TRUE, type = "dr", M = 25, name = "qvm_6_plugin") estimate(pe_plugin + pe_drql[["qvm_6"]], id = d$teacher) ################################################### ### code chunk number 23: harvard_example.Rnw:414-415 ################################################### plan("sequential")