PH724_Assignment3_VDC
Vanessa Davila
2025-04-11
Load R packages
#install.packages("devtools")
#devtools::install_version("epicalc",version="2.15.1.0")
#install.packages("epiDisplay")
library(tidyverse)
library(MatchIt)
library(jsmodule)
library(lattice)
library(epicalc)
library(cobalt)
library(Hmisc)
library(ipw)
library(survey)
library(epicalc)
library(epiDisplay)
library(plyr)
library(modelsummary)
library(knitr)
setwd('/Users/vanessa/Library/Mobile Documents/com~apple~CloudDocs/Courses/PH724-01 Advanced Methods in Epidemiology/Assignments/Assignment3_IPW_MSM')#PART 1 MOSQUITO DATASET ### Q1.1 Create a DAG to represent the causal relationship between mosquito net use and malaria risk, in the presence of other variables.
Create a Directed Acyclic Graph (DAG) to determine the associations of smoking intensity and smoking duration with systolic and diastolic blood pressure, as well as the diagnosis of high blood pressure.
Q1.2
To address confounding in this analysis, utilize the propensity score approach. Construct propensity scores for mosquito net use by developing an exposure model based on relevant variables. Justify the selection of variables for the exposure model and describe the resulting distribution of propensity scores by mosquito net use status.
## read the dataset
mosquito <- read.csv("mosquito_nets.csv")
## Calculate the propensity score for exposure
ps1 <- glm(net_num~ income + health + household + eligible + temperature, data=mosquito, family = binomial)
mosquito$psvalue1 <- predict(ps1,type = "response")
## Histogram of distribution of propensity scores by mosquito net use status
hist(mosquito$psvalue1[mosquito$net_num==1], density = 10, angle = 45, main="Propensity Scores",
xlab="Shaded = smoker | Gray = nonsmoker")
hist(mosquito$psvalue1[mosquito$net_num==0], col=gray(0.4,0.25), add=T)
## Boxplot distribution of propensity scores by mosquito net use status
lattice::bwplot(psvalue1~as.factor(net_num), data = mosquito, ylab = "Propensity Scores", auto.key = TRUE)
Potential confounders and precision variables were included in the
exposure model (as shown in the DAG). There’s good overlapping of
exposed and unexposed so we can assume a balanced distribution of PS in
both groups.
Q1.3
Perform propensity score matching and examine the balance of confounders before and after matching. In addition, calculate inverse probability weights on top of propensity score and assess the distribution of the weights. Report your assumptions, methods and results in detail.
#PS matching
#table(mosquito$net_num) 681 exposed and 1071 unexposed
out.matchit <- MatchIt::matchit(## Give formula for propensity score model)
formula = net_num~ income + health + household + eligible + temperature,
data = mosquito, # Cannot have missing values
method = "nearest", # Match using near neighbor
distance = "logit", # Distance defined by usual propensity score from logistic model
ratio = 1, # 1:1 match is the default
caliper=.25, replace=FALSE)
## Diagnostics
summary(out.matchit)##
## Call:
## MatchIt::matchit(formula = net_num ~ income + health + household +
## eligible + temperature, data = mosquito, method = "nearest",
## distance = "logit", replace = FALSE, caliper = 0.25, ratio = 1)
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.4581 0.3446 0.6170 1.9857 0.1876
## income 955.1938 872.7526 0.4089 1.3633 0.1044
## health 54.9090 48.0570 0.3619 1.2083 0.0714
## household 3.1791 2.8674 0.2097 1.2127 0.0346
## eligibleTRUE 0.0529 0.0028 0.2237 . 0.0501
## temperature 23.3809 24.0880 -0.1685 1.0846 0.0424
## eCDF Max
## distance 0.2920
## income 0.1983
## health 0.1683
## household 0.0795
## eligibleTRUE 0.0501
## temperature 0.0972
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.4197 0.4046 0.0823 1.2208 0.0210
## income 957.1843 941.0147 0.0802 1.1993 0.0214
## health 54.9755 53.6949 0.0676 1.1340 0.0151
## household 3.0277 2.9935 0.0230 1.0055 0.0067
## eligibleTRUE 0.0049 0.0049 0.0000 . 0.0000
## temperature 23.5478 23.6558 -0.0257 1.1045 0.0155
## eCDF Max Std. Pair Dist.
## distance 0.0848 0.0830
## income 0.0930 0.4476
## health 0.0522 0.6761
## household 0.0114 1.0482
## eligibleTRUE 0.0000 0.0000
## temperature 0.0473 1.0313
##
## Sample Sizes:
## Control Treated
## All 1071 681
## Matched 613 613
## Unmatched 458 68
## Discarded 0 0## density plot for covariates in the exposure model before and after matching
cobalt::bal.plot(out.matchit, var.name = 'income', which = 'both', grid=TRUE) 
## love plot for standardized mean difference
love.plot(bal.tab(out.matchit, m.threshold=0.1), stat = "mean.diffs", grid=TRUE,
stars="raw", abs = F)
# we set the m.threshold for mean difference to identify the variable(s) that don't match very well.
bal.tab(out.matchit, m.threshold=0.1) ## Balance Measures
## Type Diff.Adj M.Threshold
## distance Distance 0.0823 Balanced, <0.1
## income Contin. 0.0802 Balanced, <0.1
## health Contin. 0.0676 Balanced, <0.1
## household Contin. 0.0230 Balanced, <0.1
## eligible Binary 0.0000 Balanced, <0.1
## temperature Contin. -0.0257 Balanced, <0.1
##
## Balance tally for mean differences
## count
## Balanced, <0.1 6
## Not Balanced, >0.1 0
##
## Variable with the greatest mean difference
## Variable Diff.Adj M.Threshold
## income 0.0802 Balanced, <0.1
##
## Sample sizes
## Control Treated
## All 1071 681
## Matched 613 613
## Unmatched 458 68## Output Matched Data Sets for further analysis
dat.matchit <- match.data(out.matchit)
head(dat.matchit) # weight == 1 for matched sets and 0 for unmatched sets under the method of "nearest"#Truncated stabilized IPW
IPW <- ipwpoint(exposure=net_num, family="gaussian",
numerator= ~ 1, ## Stabilized weight / Remove numerator= ~ 1 get unstabilized weight
denominator= ~ income + health + household + eligible + temperature,
trunc = 0.01, data=mosquito)
#Compare truncated vs untruncated weights
summary(IPW$weights.trunc)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4348 0.7768 0.9008 0.9908 1.1075 2.5205# Visually assess weights
ipwplot(weights=IPW$weights.trunc,logscale=FALSE, main="Truncated stabilized weights",
xlim=c(0,6))
# Paste inverse probability weights
mosquito$weights.trunc <- IPW$weights.trunc
boxplot(weights.trunc ~ net_num, data = mosquito, main = "IPW Weights by Mosquito net use status")
The 1:1 propensity score matching summary shows small standardized mean differences. The plots comparing covariate distributions before and after matching indicate improved overlap. The love plot shows balanced covariates post-matching, with all standardized mean differences falling within ±0.1. Based on these diagnostics, the propensity score matching appears successful, and we assume conditional exchangeability of exposed and unexposed.
Also, the estimated stabilized, truncated inverse probability weights are approximately centered around 1, indicating proper scaling relative to the probability of exposure. The distribution shows limited extreme values indicating a lack of influence from a small number of observations. A lack of heavy right skewness suggests that no individuals are given disproportionately large weights, which would indicate model misspecification or poor covariate overlap. Finally, sufficient overlap in covariate distributions between exposure groups supports the validity of the weighting model, as it implies that comparisons are being made across comparable individuals.
Q1.4
Compare the effect estimates (i.e., the outcome model) from the crude linear model, the traditional linear model with confounding adjustments, the linear model with propensity score matching, and the linear model with inverse probability weighting.
library(gtsummary)
# Crude linear model
M_crude <- lm(malaria_risk ~ net_num, data = mosquito) %>%
tbl_regression(
label= list(net_num ~ "Crude model")
)%>%
add_global_p()
M_crude| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| Crude model | -16 | -18, -15 | <0.001 |
| 1 CI = Confidence Interval | |||
# Traditional linear model with confounding adjustments
M_traditional <- lm(malaria_risk ~ net_num + income + health + household + eligible + temperature, data = mosquito) %>%
tbl_regression(
label= list(net_num ~ "Traditional model")
) %>%
modify_table_body(~ .x %>% filter(variable == "net_num")) %>%
add_global_p()
M_traditional| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| Traditional model | -10 | -11, -9.9 | <0.001 |
| 1 CI = Confidence Interval | |||
# Linear regression using matched dataset
M_PS_match <- lm(malaria_risk ~ net_num + strata(subclass), data = dat.matchit) %>%
tbl_regression(
label= list(net_num ~ "PS match")
) %>%
modify_table_body(~ .x %>% filter(variable == "net_num")) %>%
add_global_p()
M_PS_match| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| PS match | -11 | -11, -10 | <0.001 |
| 1 CI = Confidence Interval | |||
# Linear model with IPW
# Estimate the outcome effect
M_IPW <- svyglm(malaria_risk ~ net_num, family = gaussian(), design=svydesign(~1,weights=~weights.trunc, data=mosquito)) %>%
tbl_regression(
label= list(net_num ~ "IPW")
) %>%
modify_table_body(~ .x %>% filter(variable == "net_num")) %>%
add_global_p()
M_IPW | Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| IPW | -11 | -12, -9.2 | <0.001 |
| 1 CI = Confidence Interval | |||
#Stacked tables
compare <- tbl_stack(
list(M_crude, M_traditional, M_PS_match, M_IPW )
) %>%
modify_header(label = "**Effect of net use on malaria risk comparing different models**")
compare| Effect of net use on malaria risk comparing different models | Beta | 95% CI1 | p-value |
|---|---|---|---|
| Crude model | -16 | -18, -15 | <0.001 |
| Traditional model | -10 | -11, -9.9 | <0.001 |
| PS match | -11 | -11, -10 | <0.001 |
| IPW | -11 | -12, -9.2 | <0.001 |
| 1 CI = Confidence Interval | |||
#PART 2 HAPPINESS DATASET ### Q2.1 All the data are longitudinal. The current exposure status may be impacted by previous exposure status and outcome, as well as time-varying confounders. Create a DAG to represent the causal relationships of the policy implementation of a 6-hour workday and the number of mandated vacation days with people’s happiness, in the presence of other variables
Q2.2
We decide to use marginal structural models with inverse probability weighting in this analysis. You may consider building different models for the two hypotheses. Please assess whether the weights are balanced along time. Report your methods and results in detail.
#Create lag variables
happy <- happy %>%
group_by(ID) %>%
dplyr::mutate(
lag_corrupt = lag(corruption, n = 1, default = 0),
lag_democr = lag(democracy, n = 1, default = 0),
lag_corrupt = replace(lag_corrupt, 1, corruption[1]),
lag_democr = replace(lag_democr, 1, democracy[1]),
lag_policy = lag(policy, n = 1, default = 0),
lag_vacdays = lag(vacation_days, n = 1, default = 0),
lag_happolicy = lag(happiness_policy, n = 1, default = NA),
lag_happolicy = replace(lag_happolicy, 1, happiness_policy[1]),
lag_hapvac = lag(happiness_vacation, n = 1, default = NA),
lag_hapvac = replace(lag_hapvac, 1, happiness_vacation[1])
)
#Exposure model weight with IPW
# use GEE models for the numerator and denominator, and account for autoregressive time structures in the data
#Continuous exposure
w1 <- ipwtm(
exposure = vacation_days,
family = "gaussian",
corstr = "ar1",
numerator = ~ lag_vacdays + as.factor(year),
denominator = ~ lag_vacdays + lag_hapvac + lag_corrupt + lag_democr + as.factor(year),
id = ID,
timevar=year,
type="all",
data = as.data.frame(happy))
## 1.Weights distribution
happy$ipw1 <- w1$ipw.weights
ipwplot(weights = happy$ipw1,
timevar = happy$year,
binwidth = 0.5,
main = "Stabilized weights vacation days")
#Dichotomous exposure
w2 <- ipwtm(
exposure = policy,
family = "binomial",
link = "logit",
numerator = ~ lag_policy + as.factor(year),
denominator = ~ lag_policy + lag_happolicy + lag_corrupt + lag_democr + as.factor(year),
id = ID,
timevar=year,
type="all",
data = as.data.frame(happy))
## 2.Weights distribution
happy$ipw2 <- w2$ipw.weights
ipwplot(weights = happy$ipw2,
timevar = happy$year,
binwidth = 0.5,
main = "Stabilized weights policy of working hours")
The plots show the distribution of stabilized inverse probability
weights from 2010 to 2019, estimated using marginal structural models to
evaluate the effect of vacation days (continuous exposure) and policy
(binary exposure) on happiness. The numerator models include exposure
and time-invariant covariates, while the denominator models additionally
incorporate time-varying covariates. Across all years, the weights are
centered, indicating stability over time, with moderate variability and
few extreme values. This supports the appropriateness of the weights for
use in outcome modeling. For binary exposures, stabilized weights are
optional but often preferred for interpretability. For continuous
exposures, unstabilized weights are not usable, as they result in
infinite variance; stabilized weights are required. With the weights
validated, the next step is to build the outcome model.
Q2.3
Outcome models. Describe your assumptions and interpret the causal relationships of interest.
#install.packages("geepack")
library(geepack)
library(gtsummary)
#Model happines-vacation days
happy_vac <- geeglm(happiness_vacation ~ vacation_days + as.factor(year), id=ID, data=happy,
family=gaussian("identity"), corstr="ar1",
weights=ipw1) %>%
tbl_regression()%>%
modify_caption("Causal effect of vacation days on people's happiness")
happy_vac| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| vacation_days | 1.3 | 1.0, 1.6 | <0.001 |
| as.factor(year) | |||
| 2010 | — | — | |
| 2011 | 0.70 | 0.41, 0.98 | <0.001 |
| 2012 | 2.5 | 1.9, 3.1 | <0.001 |
| 2013 | 2.8 | 2.1, 3.6 | <0.001 |
| 2014 | 3.5 | 2.5, 4.5 | <0.001 |
| 2015 | 4.6 | 3.4, 5.8 | <0.001 |
| 2016 | 5.1 | 3.7, 6.4 | <0.001 |
| 2017 | 6.4 | 4.8, 8.0 | <0.001 |
| 2018 | 7.6 | 5.5, 9.7 | <0.001 |
| 2019 | 8.3 | 5.9, 11 | <0.001 |
| 1 CI = Confidence Interval | |||
#Model happines-policy of working hours
happy_policy <- geeglm(happiness_policy ~ policy + as.factor(year), id = ID, family = gaussian,
corstr = "ar1",
weights=ipw2, data = happy) %>%
tbl_regression()%>%
modify_caption("Causal effect of policy of working hours on people's happiness")
happy_policy| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| policy | 6.8 | 5.9, 7.7 | <0.001 |
| as.factor(year) | |||
| 2010 | — | — | |
| 2011 | 1.2 | 0.58, 1.8 | <0.001 |
| 2012 | 3.6 | 2.9, 4.4 | <0.001 |
| 2013 | 3.7 | 2.9, 4.6 | <0.001 |
| 2014 | 4.7 | 3.6, 5.8 | <0.001 |
| 2015 | 6.1 | 4.9, 7.3 | <0.001 |
| 2016 | 7.1 | 5.8, 8.3 | <0.001 |
| 2017 | 8.8 | 7.5, 10 | <0.001 |
| 2018 | 11 | 9.1, 12 | <0.001 |
| 2019 | 11 | 9.7, 13 | <0.001 |
| 1 CI = Confidence Interval | |||
The marginal structural models using Generalized Estimating Equations (GEE) estimate the causal effects of vacation days and work-hour policies on happiness over time. In Model 1, each additional vacation day is associated with a 1.3-point increase in happiness (95% CI: 1.0 to 1.6, p < 0.001), after adjusting for time-varying confounding. The year coefficients suggest a consistent upward trend in happiness from 2011 to 2019.
In Model 2, the presence of a work-hour policy is associated with a 6.8-point higher happiness score (95% CI: 5.9 to 7.7, p < 0.001). Similar to Model 1, happiness increases steadily over time. Both models show strong temporal effects, but the effect size for policy is notably larger than that for vacation days, suggesting a more substantial impact of structural policy on happiness.
Assuming no unmeasured confounding and no selection bias, we can intepret the results from these models as the causal effect of policy and vacation days on people’s happiness.
#PART 3 NATURAL EXPERIMENTAL STUDIES ### Q3.1 In the methods section, the authors described the baseline temporal trends of the outcome in states that did or did not implement same-sex marriage policies by 2015. They also tested this assumption with a linear regression analysis of suicide attempts before implementation of same-sex marriage policy by states. The figure below shows the baseline temporal trends among “never legal”, “Wave 1 before registration”, and “Wave 2 before registration”. In the results section, they also reported that “Before the implementation of same-sex marriage policies, temporal trends in suicide attempts were not significantly different in states with same-sex marriage policies relative to temporal trends in suicide attempts in states without same-sex marriage policies (0.001 percentage points; 95%CI, more than –0.001 to 0.001)”. What assumption did the authors try to evaluate before the difference-and difference analysis? Was this assumption fulfilled based on their reports? Briefly explain why.
The authors aimed to evaluate the parallel trends assumption between exposed states (those with same-sex marriage policy) and unexposed states (those without). This assumption implies a stable difference in outcomes between groups over time. The observed trends suggest otherwise: in both waves 1 and 2, suicide attempts were higher around 2000, but after 2003, the trends diverge—wave 1 shows a decline, while wave 2 shows an increase relative to states where same-sex marriage was never legal. This suggests a potential violation of the assumption, even if the authors note that the differences were not statistically significant. Furthermore, although the figure distinguishes between waves, this separation is not reflected in the interpretation or conclusions.
Q3.2
In this model, same sex marriage state is a binary indicator for with or without implementation of same sex marriage policies (i.e., Wave 1 and Wave 2 are combined in the implementation category). Policy enactment is a binary indicator for before and after the implementation of same sex marriage policies. What are the interpretations for 𝛽1, 𝛽2, and 𝛽3 in this study? Does it make sense to control for annual state unemployment rates and state-level policies preventing employment discrimination on the basis of sexual orientation in this analysis? Briefly explain why.
B1: difference in suicide rates before and after the intervention in the unexposed group B2: difference in suicide rates among the 2 groups before and after the policy enactment B3: difference in suicide rates pre- and postintervention, comparing states with policy enactment (exposed) and states with policy enactment (unexposed). It only makes sense to control for annual state unemployment rates if this time-variant characteristic remains unvariant within groups. It only makes sense to control for state-level policies preventing employment discrimination if those were implemented more or less simultaneosly in all states but not at the same time that the same-sex marriage policies.
#PART 4 STUDY DESIGN Since its establishment in 1979, the Department of Education has overseen various levels of education from elementary schools to universities, primarily focusing on public schools. The Department’s primary responsibilities include providing financial aid, ensuring equal access, preventing discrimination, fostering national partnerships, and more. One notable initiative supported by the Department of Education is the Early College High School (ECHS) model, which integrates high school and college education. This model allows students to earn both a high school diploma and college credits simultaneously, offering rigorous college-preparatory courses, dual enrollment opportunities, and comprehensive academic and social-emotional support at little to no cost to students. The ECHS model, formally started in 2002, focuses particularly on students from low-income backgrounds and underrepresented racial and ethnic groups, addressing barriers to college attainment. Please use your skills in natural experimental design to evaluate whether the ECHS model improves college credit accumulation.
Q4.1
Define your target population. Assuming you have access to all historical data and resources to collect additional prospective data, how would you measure exposure, outcomes, and other relevant variables?
target population: students from 15-21 years old exposure: SES (low vs high) intervention: ECHS model implementation (2002) outcome(s): college credit accumulation (continuos), formal employment 1year after ECHS completion, proportion of students in ECHS that pursue a graduate degree other relevant variables: major, race/ethnicity, financial aid status, migration status, gender, marital status, housing stability, drug use, US region or state, recession periods, inflation
Q4.2
What specific research question would you like to answer? Based on your research question, which natural experimental study design would you choose? Justify your decision, provide details of your methods, and discuss the interpretation of potential results and how you would draw conclusions.
I think based on the exposure and outcomes, several research questions can be formulated and responded using difference-in-differences and interrupted time series. I think the interpretation is not as complicated as other natural experiments, and both methodologies are appropiated to address the propposed research questions:
##Difference-in-Differences Research Questions
- Did ECHS implementation reduce disparities in college credit accumulation between low- and high-SES students aged 15–21, compared to the period before implementation?
Rationale:Low-SES students historically accumulate fewer college credits during high school due to limited access to dual enrollment programs, lack of academic advising, and financial constraints. ECHS programs are designed to provide structured access to college coursework, tuition waivers, and support services, which may disproportionately benefit low-SES students.
Expected Results:Following ECHS implementation, the gap in college credit accumulation between low- and high-SES students is expected to narrow, indicating that the intervention helped improve postsecondary readiness among disadvantaged youth.
- Following the introduction of the ECHS model, did the gap in formal employment one year after high school (or college) completion between low- and high-SES students change significantly over time?
Rationale:ECHS programs often include career-oriented coursework and credentials, which may improve employability. Since low-SES students face more barriers to labor market entry, ECHS could reduce employment disparities by providing early exposure to workforce pathways and reducing time to postsecondary credentialing. Expected results: We would expect a bigger difference in employment status 1 year after program completion, with the prop. of employment to remain relatively constant among high SES students and have a greater increase after the intervention among low SES students.
- Did ECHS implementation increase the likelihood of pursuing graduate education among low-SES students relative to high-SES students, compared to pre-intervention trends?
Rationale: Early exposure to college coursework and credentials may raise academic aspirations among low-SES students and reduce time and cost barriers associated with extended educational trajectories. While high-SES students are already more likely to pursue graduate study, ECHS may have a greater marginal impact on low-SES students’ access to graduate education.
Expected Results: An increase in graduate degree pursuit among low-SES students post-ECHS implementation, narrowing the SES gap, though the effect may be smaller or more delayed than for earlier outcomes (e.g., credit accumulation or employment).
##Interrupted Time Series Research Questions
- Does the introduction of ECHS in a given district/state/school produce a significant shift in the trend of average college credit accumulation among 15–21-year-old students in general? Among low SES students?
Rationale:ECHS programs are designed to increase access to and completion of college coursework during high school, particularly through structured pathways, academic advising, and financial support. These programs lower institutional barriers to dual enrollment and may change the trajectory of average credit accumulation, especially for students with limited prior access.
Expected Results:A positive level change (immediate jump) and/or upward shift in trend (sustained increase over time) in average college credit accumulation after ECHS implementation. The effect may vary depending on the intensity and fidelity of implementation at the district or state level. It would be of special interest to analyze the trend among States/districts with low and high proportions of schools following ECHS model.
- Was there a detectable level or trend change in formal employment one year after ECHS completion following the implementation of the ECHS model?
Rationale:By accelerating postsecondary credential attainment and exposing students to workforce-relevant courses earlier, ECHS may improve employment outcomes. ITS can detect whether implementation coincides with structural changes in labor market entry patterns for recent graduates.
Expected Results:A moderate level increase in employment rates one year after high school for ECHS completers, with possible continued upward trend if the labor market advantages persist. The change may appear with a lag, depending on time to degree or credential.
- After the 2002 introduction of ECHS, does the proportion of students pursuing graduate education change significantly over time within implementing regions?
Rationale:ECHS may influence long-term educational trajectories by building early academic momentum and increasing confidence and aspiration among students, particularly first-generation or low-income youth. ITS can capture delayed structural changes in graduate education enrollment.
Expected Results:A delayed but positive trend change in the proportion of students pursuing graduate education in regions with sustained ECHS adoption. A level shift may not be immediate, but a gradual increase is expected if ECHS improves degree completion and continuation rates.