| Title | IPW models |
|---|---|
| author | TingFang Lee, URI Avenir Team |
Load packages and function script...
library(lme4)
library(plyr)
library(igraph)
library(numDeriv)
library(dplyr)
source("https://github.com/uri-ncipher/Nearest-Neighbor-estimators/blob/main/IPW.R?raw=TRUE")
Read in the synthetic data and network structure
sim.data <- read.csv("https://github.com/uri-ncipher/Nearest-Neighbor-estimators/blob/main/sim.data.csv?raw=TRUE")
net0 <- readRDS(file = url("https://github.com/uri-ncipher/Nearest-Neighbor-estimators/blob/main/net0.100.rds?raw=TRUE", open = "rb"))
The synthetic dataset sim.data and network net0 include observations with missing outcomes. These data can be analyzed using the IPCW estimators, which account for outcome missingness.
Because the IPW estimators require complete outcome data, we first calculate each individual's number of neighbors and number of treated neighbors in the original network. We then remove individuals with missing outcomes from the network, along with any resulting isolated nodes, to create the complete-case network net1. The number of neighbors and number of treated neighbors are then recalculated based on this reduced network for use in the IPW analysis.
Finally, for both the original and complete-case networks, we calculate the average covariate value among each individual's neighbors (avg_L), which is used in the subsequent IPCW2 and IPW2.
V(net0)$name <- as.character(V(net0))
remove.ids <- as.character(sim.data$id[sim.data$censored.cor == 1])
net1 <- delete_vertices(net0, remove.ids)
net1 <- delete_vertices(net1, V(net1)[degree(net1) == 0])
keep.ids <- as.numeric(V(net1)$name)
sim.data1 <- sim.data[match(keep.ids, sim.data$id), ]
sim.data1$n.neighbor <- degree(net1)
sim.data1$n.trt.neighbor <- sapply(sim.data1$id, function(i) {
nb <- neighbors(net1, as.character(i))
nb.ids <- as.numeric(V(net1)[nb]$name)
sum(sim.data1$A[match(nb.ids, sim.data1$id)])
})
nn.sim.data <- data.frame()
for (i in sim.data[["id"]]) {
nb <- neighborhood(net0, order = 1, nodes = as.character(i))[[1]]
ids <- as.numeric(V(net0)[nb]$name)
ids <- setdiff(ids, i)
temp.data <- sim.data[sim.data[["id"]] %in% ids, ]
temp.data$ego <- i
nn.sim.data <- rbind(nn.sim.data, temp.data)
}
nn.sim.data <- nn.sim.data %>%
group_by(ego) %>%
dplyr::summarise(avg_L = mean(L), .groups = "drop")
sim.data <- merge(sim.data, nn.sim.data, by.x = "id", by.y = "ego")
nn.sim.data.1 <- data.frame()
for (i in sim.data1[["id"]]) {
nb <- neighborhood(net1, order = 1, nodes = as.character(i))[[1]]
ids <- as.numeric(V(net1)[nb]$name)
ids <- setdiff(ids, i)
temp.data <- sim.data1[sim.data1[["id"]] %in% ids, ]
temp.data$ego <- i
nn.sim.data.1 <- rbind(nn.sim.data.1, temp.data)
}
nn.sim.data.1 <- nn.sim.data.1 %>%
group_by(ego) %>%
dplyr::summarise(avg_L = mean(L), .groups = "drop")
sim.data1 <- merge(sim.data1, nn.sim.data.1, by.x = "id", by.y = "ego")
In this section, we apply the IPW1, IPW2, IPCW1, and IPCW2 estimators to estimate average potential outcomes and the corresponding causal effects. The average potential outcome is defined by an individual's exposure/treatment status and the allocation strategy (\alpha), where (\alpha) represents the probability that individuals in the nearest-neighbor set receive the exposure/treatment under the counterfactual allocation strategy. Thus, (\alpha) takes a value between 0 and 1.
For each estimator, the R functions return the point estimates and estimated variances of the average potential outcomes, (\boldsymbol{\widehat{Y}(1,\alpha), \widehat{Y}(0,\alpha), \widehat{Y}(\alpha)}), under three allocation strategies, (\alpha = 0.25, 0.50,) and (0.75). The functions also return the point estimates and estimated variances of four causal effects: the direct effect (DE), indirect effect (IE), total effect (TE), and overall effect (OE).
The IPW1 and IPW2 estimators are applied to complete outcome data, whereas the IPCW1 and IPCW2 estimators accommodate missing outcomes through inverse probability of censoring weighting.
Equations for calculating the four causal effects are:
${\widehat{DE} = \widehat{Y}(1, \alpha) - \widehat{Y}(0, \alpha)}$ ${\widehat{IE} = \widehat{Y}(0, \alpha_0) - \widehat{Y}(0, \alpha_1)}$ ${\widehat{TE} = \widehat{Y}(1, \alpha_0) - \widehat{Y}(0, \alpha_1)}$ -
${\widehat{OE} = \widehat{Y}(\alpha_0) - \widehat{Y}(\alpha_1)}$ ,
where
- (i)
$\alpha_0 = 0.75$ versus$\alpha_1 = 0.25$ ; - (ii)
$\alpha_0 = 0.75$ versus$\alpha_1 = 0.50$ ; - (iii)
$\alpha_0 = 0.25$ versus$\alpha_1 = 0.50$ .
point.est.IPCW1 = Y_IPCW1(data = sim.data, net = net0, alphas = c(0.25, 0.5, 0.75), id.var = "id", trt.var = "A", n.neighbor = "n.neighbor", n.trt.neighbor = "n.trt.neighbor", propensity.covars = c("L"), censor.indicator = "censored.cor", censor.covars = c("L"), censor.correlation = TRUE, component.var = "component", obs.y = "obs.y.cor")
var.est.IPCW1 = Var_Y_IPCW1(alphas = c(0.25, 0.5, 0.75), net = net0, data = sim.data, id.var = "id", obs.y= "obs.y.cor", n.neighbor = "n.neighbor", n.trt.neighbor = "n.trt.neighbor", trt.var = "A", propensity.covars = c("L"), censor.indicator = "censored.cor", censor.covars = c("L"), censor.correlation = TRUE, component.var = "component")
point.est.IPCW2 = Y_IPCW2(data = sim.data, net = net0, alphas = c(0.25, 0.5, 0.75), id.var = "id", trt.var = "A", n.neighbor = "n.neighbor", n.trt.neighbor = "n.trt.neighbor", propensity.covars = c("L"), censor.covars = c("L"), censor.correlation = TRUE, component.var = "component", obs.y = "obs.y.cor")
var.est.IPCW2 = Var_Y_IPCW2(alphas = c(0.25, 0.5, 0.75), net = net0, data = sim.data, id.var = "id", obs.y= "obs.y.cor", n.neighbor = "n.neighbor", n.trt.neighbor = "n.trt.neighbor", trt.var = "A", propensity.covars = c("L"), propensity.nn.covars = c("avg_L"), censor.indicator = "censored.cor", censor.covars = c("L"), censor.correlation = TRUE, component.var = "component")
point.est.IPW1 = Y_IPW1(data = sim.data1, net = net1, alphas = c(0.25, 0.5, 0.75), id.var = "id", trt.var = "A", n.neighbor = "n.neighbor", n.trt.neighbor = "n.trt.neighbor", propensity.covars = c("L"), obs.y = "obs.y.cor")
var.est.IPW1 = Var_Y_IPW1(alphas = c(0.25, 0.5, 0.75), net = net1, data = sim.data1, id.var = "id", obs.y= "obs.y.cor", n.neighbor = "n.neighbor", n.trt.neighbor = "n.trt.neighbor", trt.var = "A", propensity.covars = c("L"), component.var = "component")
point.est.IPW2 = Y_IPW2(data = sim.data1, net = net1, alphas = c(0.25, 0.5, 0.75), id.var = "id", trt.var = "A", n.neighbor = "n.neighbor", n.trt.neighbor = "n.trt.neighbor", propensity.nn.covars = c("avg_L"), propensity.covars = c("L"), obs.y = "obs.y.cor")
var.est.IPW2 = Var_Y_IPW2(alphas = c(0.25, 0.5, 0.75), net = net1, data = sim.data1, id.var = "id", obs.y= "obs.y.cor", n.neighbor = "n.neighbor", n.trt.neighbor = "n.trt.neighbor", trt.var = "A", propensity.covars = c("L"), propensity.nn.covars = c("avg_L"), component.var = "component")
- Lee T, Buchanan AL, Katenka NV, Forastiere L, Halloran ME, Friedman SR, Nikolopoulos G. Estimating causal effects of HIV prevention interventions with interference in network-based studies among people who inject drugs. Annals of Applied Statistics. 2023;17(3):2165–2191. doi:10.1214/22-AOAS1713.
- Lee T, Buchanan AL, Katenka N, Forastiere L, Halloran ME, Nikolopoulos G. Assessing spillover effects: Handling missing outcomes in network-based studies. Statistical Methods in Medical Research. 2025;34(12):2284–2301. doi:10.1177/09622802251382586.