# A g-computation [Validation Overall]

## A.1 Parametric g-formula

Based on Hernan and Robins, 2020 Chapter 13

Zou’s modified Poisson regression code

## A.2 NEW FUNCTIONS

``````# Helpers g-computation

# standardization =============================================================

standardization <- function(data, indices, treatment, outcome, formula1) {

treatment1 <- enquo(treatment)
outcome1 <- enquo(outcome)

# create a dataset with 3 copies of each subject
# 1st copy: equal to original one`
d <- data[indices, ]
d <- d %>% mutate(interv = -1,
seqno = 1:n()) # for gee estimation
# 2nd copy: treatment set to 0, outcome to missing
d0 <- d %>% mutate(interv = 0,
!!treatment1 := first(levels(droplevels(!!treatment1))),
!!outcome1 := NA)
# 3rd copy: treatment set to 1, outcome to missing
d1 <- d %>% mutate(interv = 1,
!!treatment1 := last(levels(droplevels(!!treatment1))),
!!outcome1 := NA)
d.onesample <- rbind(d, d0, d1) # combining datasets

if (class(unlist(select(data, !!outcome1)))=="factor") {
d.onesample <- d.onesample %>%
mutate(!!outcome1 := as.numeric(!!outcome1)-1)
}

# linear model to estimate mean outcome conditional on treatment and confounders
# parameters are estimated using original observations only (interv= -1)
# parameter estimates are used to predict mean outcome for observations with set
# treatment (interv=0 and interv=1)
# fit <- glm(formula = formula1 ,data = d.onesample, family = "binomial")
# d.onesample\$predicted_meanY <- predict(fit, d.onesample, type="response")

# Zou, 2004
require(geepack)
fit <- geeglm(formula = formula1,
data    = d.onesample,
family  = poisson(link = "log"),
id      = seqno,
corstr  = "exchangeable")

d.onesample\$predicted_meanY <- predict(fit, d.onesample, type="response")

# estimate mean outcome in each of the groups interv=-1, interv=0, and interv=1
M.ob <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == -1], na.rm=T)
M.t0 <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == 0], na.rm=T)
M.t1 <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == 1], na.rm=T)
M.t10 <- M.t1 - M.t0

return(c(M.ob, M.t0, M.t1, M.t10))
}

# plot_eff =============================================================

plot_eff <- function(dat1, dat2, dat3, coefs, legend.title = "Total Effect", model.names = c("Overall", "Rural", "Periurban"), size = 1, base_size = 13, xlims = c(-.3,1.3)) {
library(ggsci)

dat1\$coeftable %>% as_tibble(rownames = "term") %>% mutate(model = model.names[1]) %>%
rbind(dat2\$coeftable %>% as_tibble(rownames = "term") %>% mutate(model = model.names[2])) %>%
rbind(dat3\$coeftable %>% as_tibble(rownames = "term") %>% mutate(model = model.names[3])) %>%
clean_names() %>%
filter(term %in% coefs) %>%
ggplot(aes(y = model, x = est, col = model)) +
geom_pointrange(aes(xmin = x2_5_percent, xmax = x97_5_percent), size=size) +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_color_npg(limits = model.names) +
scale_y_discrete(limits = rev(model.names)) +
scale_x_continuous(limits = xlims) +
theme_bw(base_size = base_size) +
guides(color = guide_legend(title.position = "top",
title.hjust = 0)) +
labs(y = "", x = "Estimate", col = legend.title) +
theme(legend.position = "top")
}

# standardization_sim =============================================================

standardization_sim <- function(data, indices, treatment, outcome, formula1, p_out) {

treatment1 <- enquo(treatment)
outcome1 <- enquo(outcome)

# create a dataset with 3 copies of each subject
# 1st copy: equal to original one`
d <- data[indices, ]
d <- d %>% mutate(interv = -1,
seqno = 1:n()) # for gee estimation
# 2nd copy: treatment set to 0, outcome to missing
d0 <- d %>% mutate(interv = 0,
!!treatment1 := first(levels(droplevels(!!treatment1))),
!!outcome1 := NA)
# 3rd copy: treatment set to 1, outcome to missing
d1 <- d %>% mutate(interv = 1,
!!treatment1 := last(levels(droplevels(!!treatment1))),
!!outcome1 := NA)
# 4th copy: treatment set to simulated prop, outcome to missing
d2 <- d %>% mutate(interv = 2,
!!treatment1 := rbinom(n=nrow(d), size=1, prob=p_out),
!!treatment1 := ifelse(!!treatment1 == 1,
"Outside", "Inside"),
!!outcome1 := NA)
d.onesample <- rbind(d, d0, d1, d2) # combining datasets

if (class(unlist(select(data, !!outcome1)))=="factor") {
d.onesample <- d.onesample %>%
mutate(!!outcome1 := as.numeric(!!outcome1)-1)
}

# linear model to estimate mean outcome conditional on treatment and confounders
# parameters are estimated using original observations only (interv= -1)
# parameter estimates are used to predict mean outcome for observations with set
# treatment (interv=0 and interv=1)
# fit <- glm(formula = formula1 ,data = d.onesample, family = "binomial")
# d.onesample\$predicted_meanY <- predict(fit, d.onesample, type="response")

# Zou, 2004
require(geepack)
fit <- geeglm(formula = formula1,
data    = d.onesample,
family  = poisson(link = "log"),
id      = seqno,
corstr  = "exchangeable")

d.onesample\$predicted_meanY <- predict(fit, d.onesample, type="response")

# estimate mean outcome in each of the groups interv=-1, interv=0, and interv=1
M.ob <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == -1], na.rm=T)
M.t0 <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == 0], na.rm=T)
M.t1 <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == 1], na.rm=T)
M.ts <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == 2], na.rm=T)
M.t10 <- M.t1 - M.t0
M.t1s <- M.t1 - M.ts
M.ts0 <- M.ts - M.t0
M.tsob <- M.ts - M.ob

return(c(M.ob, M.t0, M.t1, M.t10, M.ts, M.t1s, M.ts0, M.tsob))
}

# tab_sim =============================================================

tab_sim <- function(results) {
se <- c(sd(results\$t[, 1]),
sd(results\$t[, 2]),
sd(results\$t[, 3]),
sd(results\$t[, 4]),
sd(results\$t[, 5]),
sd(results\$t[, 6]),
sd(results\$t[, 7]),
sd(results\$t[, 8]))
mean <- results\$t0
ll <- mean - qnorm(0.975) * se
ul <- mean + qnorm(0.975) * se

# OBS = observed, NE = No Exposed, FE = Exposed, SE = Simulated Exposure

bootstrap <-data.frame(cbind(type=c("OBS",
"NE",
"FE",
"FE - NE",
"SE",
"FE - SE",
"SE - NE",
"SE - OBS"),
mean, se, ll, ul)) %>%
mutate_at(.vars = c(2:5), ~as.numeric(as.character(.)))

return(bootstrap)

}

# standardization_bin =============================================================

standardization_bin <- function(data, indices, treatment, outcome, formula1, var_sim, cat_sim) {

treatment1 <- enquo(treatment)
outcome1 <- enquo(outcome)
var_sim1 <- enquo(var_sim)

# create a dataset with 3 copies of each subject
# 1st copy: equal to original one`
d <- data[indices, ]
d <- d %>% mutate(interv = -1,
seqno = 1:n(), # for gee estimation
cat = !!var_sim1)
# # 2nd copy: treatment set to 0, outcome to missing
# d0 <- d %>% mutate(interv = 0,
#                    !!treatment1 := first(levels(droplevels(!!treatment1))),
#                    !!outcome1 := NA)
# # 3rd copy: treatment set to 1, outcome to missing
# d1 <- d %>% mutate(interv = 1,
#                    !!treatment1 := last(levels(droplevels(!!treatment1))),
#                    !!outcome1 := NA)
# 4th copy: treatment set to 0 for cat_sim, outcome to missing
d2 <- d %>% mutate(interv = 2,
!!treatment1 := as.character(!!treatment1),
!!treatment1 := ifelse(cat == cat_sim,
"Inside", !!treatment1),
!!treatment1 := as.factor(!!treatment1),
!!outcome1 := NA)
# 5th copy: treatment set to 1 for cat_sim, outcome to missing
d3 <- d %>% mutate(interv = 3,
!!treatment1 := as.character(!!treatment1),
!!treatment1 := ifelse(cat != cat_sim,
"Outside", !!treatment1),
!!treatment1 := as.factor(!!treatment1),
!!outcome1 := NA)
d.onesample <- rbind(d, d2, d3) # combining datasets

if (class(unlist(select(data, !!outcome1)))=="factor") {
d.onesample <- d.onesample %>%
mutate(!!outcome1 := as.numeric(!!outcome1)-1)
}

# linear model to estimate mean outcome conditional on treatment and confounders
# parameters are estimated using original observations only (interv= -1)
# parameter estimates are used to predict mean outcome for observations with set
# treatment (interv=0 and interv=1)
# fit <- glm(formula = formula1 ,data = d.onesample, family = "binomial")
# d.onesample\$predicted_meanY <- predict(fit, d.onesample, type="response")

# Zou, 2004
require(geepack)
fit <- geeglm(formula = formula1,
data    = d.onesample,
family  = poisson(link = "log"),
id      = seqno,
corstr  = "exchangeable")

d.onesample\$predicted_meanY <- predict(fit, d.onesample, type="response")

# estimate mean outcome in each of the groups interv=-1, interv=0, and interv=1
M.ob <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == -1], na.rm=T)
M.tc0 <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == 2], na.rm=T)
M.tc1 <- mean(d.onesample\$predicted_meanY[d.onesample\$interv == 3], na.rm=T)
M.tc10 <- M.tc1 - M.tc0
M.tc1ob <- M.tc1 - M.ob
M.tc0ob <- M.tc0 - M.ob

return(c(M.ob, M.tc0, M.tc1, M.tc10, M.tc1ob, M.tc0ob))
}

# tab_bin =============================================================

tab_bin <- function(results, cat_sim) {
se <- c(sd(results\$t[, 1]),
sd(results\$t[, 2]),
sd(results\$t[, 3]),
sd(results\$t[, 4]),
sd(results\$t[, 5]),
sd(results\$t[, 6]))
mean <- results\$t0
ll <- mean - qnorm(0.975) * se
ul <- mean + qnorm(0.975) * se

# OBS = observed, NE = No Exposed, FE = Exposed

bootstrap <-data.frame(cbind(type=c("OBS",
"NE",
"FE",
"FE - NE",
"FE - OBS",
"NE - OBS"),
mean, se, ll, ul)) %>%
mutate_at(.vars = c(2:5), ~as.numeric(as.character(.)))

return(bootstrap)

}``````

### A.2.1 ESTIMATION

``````f <- SEROPOSITIVE ~ work_out + edad + nm_sex

library(boot)

# bootstrap
results <- boot(data = dat,
statistic = standardization,
treatment = work_out,
outcome = SEROPOSITIVE,
formula1 = f, #parameters from function
R = 5)

# generating confidence intervals
se <- c(sd(results\$t[, 1]),
sd(results\$t[, 2]),
sd(results\$t[, 3]),
sd(results\$t[, 4]))
mean <- results\$t0
ll <- mean - qnorm(0.975) * se
ul <- mean + qnorm(0.975) * se

bootstrap <-data.frame(cbind(type=c("Observed",
"No Treatment",
"Treatment",
"Treatment - No Treatment"),
mean, se, ll, ul)) %>%
mutate_at(.vars = c(2:5), ~as.numeric(as.character(.)))

kable(bootstrap, digits = 5)``````
type mean se ll ul
Observed 0.48771 0.01295 0.46233 0.51309
No Treatment 0.41185 0.01626 0.37998 0.44371
Treatment 0.62928 0.01428 0.60130 0.65727
Treatment - No Treatment 0.21744 0.01519 0.18766 0.24722

## A.3 Benchmark [RISCA package]

To compute Marginal effects of the treatment (ATE) based on RISCA package

### A.3.1 Using logistic regression as the Q-model

``````library(RISCA)

dat_gf <- dat %>%
dplyr::select(travel_1m, work_out, edad, SEROPOSITIVE, nm_sex) %>%
mutate(id = seq(1:n()),
time = 0,
travel_1m = as.integer(as.numeric(travel_1m)-1),
work_out = as.numeric(work_out)-1,
SEROPOSITIVE = as.numeric(SEROPOSITIVE)-1,
filter(complete.cases(.)) %>%
as.data.frame()

glm.multi <- glm(formula = f, data=dat_gf,

gc.ate1 <- gc.logistic(glm.obj=glm.multi, data=dat_gf,
group = "work_out", effect="ATE",
var.method="bootstrap", iterations=1000,
n.cluster=1)

gc.ate1.dat <-data.frame(cbind(type=c("No Treatment",
"Treatment",
"Treatment - No Treatment"),
bind_rows(gc.ate1\$p0, gc.ate1\$p1,
gc.ate1\$delta)))

kable(gc.ate1.dat, digits = 5)``````
type estimate ci.lower ci.upper std.error
2.5%…1 No Treatment 0.42545 0.39665 0.45410 NA
2.5%…2 Treatment 0.65798 0.60958 0.70779 NA
2.5%…3 Treatment - No Treatment 0.23252 0.17215 0.28979 0.03071