-
Notifications
You must be signed in to change notification settings - Fork 1
/
13-applied-ipw.Rmd
211 lines (138 loc) · 8.61 KB
/
13-applied-ipw.Rmd
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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
```{r 13_setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, eval=FALSE, fig.align="center")
```
# Applied Analysis: IPW
## Learning Goals {-}
- Conduct and interpret results from an appropriate IPW analysis to estimate causal effects and effect modification of causal effects.
<br>
Slides from today are available [here](https://docs.google.com/presentation/d/1gclqUwHE67Q1dBv8l7TfB_-ePfYBq0v143ycJoKohRc/edit?usp=sharing).
<br><br><br>
## Analysis {-}
**You can download a template RMarkdown file to start from [here](template_rmds/13-applied-ipw.Rmd).**
### Data context and research questions {-}
**Research questions:** Does a policy that delays the tenure process (for new parents) actually help tenure outcomes? Does it help men and women equally?
The data and full codebook are available [on Moodle](https://moodle.macalester.edu/mod/folder/view.php?id=36846) (`aer_primarysample.csv` and `ReadMe.pdf`). An abbreviated codebook is below:
**Treatment and primary outcome:**
- `gncs`: Is a gender-neutral clock-stopping policy in place at this person's school? (treatment)
- `tenure_policy_school`: Did this person get tenure at the university with the clock-stopping policy? (primary outcome)
<br>
**Additional outcomes (mediators):**
- `top_pubsX`: cumulative number of peer-reviewed publications in top-5 journals by year X (3, 5, 7, and 9) since PhD completion
- `PUBSX`: cumulative number of non-top-5 peer-reviewed publications by year X (3, 5 7, and 9) since PhD completion
<br>
**Additional covariates:**
- `female`: Is the faculty member a female?
- `pol_u`: policy university identifier (values are randomized for privacy reasons)
- `pol_job_start`: identifier for year the job at the policy university started (values are randomized for privacy reasons)
- `phd_rank`: 1-5 categorical variable of PhD program tier based on placements into the top-50 departments in our sample
- `phd_rank_miss`: indicator for missing PhD information
- `post_doc`: indicator for doing a post-doc before the first tenure-track position
- `ug_students`: number of undergraduate students at the university, in thousands
- `grad_students`: number of graduate students at the university, in thousands
- `faculty`: number of faculty at the university (all disciplines), in hundreds
- `full_av_salary`: average salary of full professors at the university, in thousands
- `assist_av_salary`: average salary of assistant professors at the university, in thousands,
- `revenue`: annual revenue of the university, in 10,000,000s
- `female_ratio`: fraction of the faculty at the university who are female
- `full_ratio`: fraction of the faculty at the university who are full professors
- `RANK`: equal to 1 for top-10 departments and 2 for all other departments
- `max_csstops`: number of children born within five years of PhD completion
### Loading data and packages {-}
You'll need to install the `survey` package first.
```{r}
library(survey) # install.packages("survey")
library(tidyverse)
library(splines) # Make sure that the splines package is installed before running
tenure <- read_csv("aer_primarysample.csv")
tenure <- tenure %>%
mutate(
faculty = ifelse(faculty_miss==1, NA, faculty),
revenue = ifelse(revenue_miss==1, NA, revenue),
female_ratio = ifelse(female_ratio_miss==1, NA, female_ratio),
full_ratio = ifelse(full_ratio_miss==1, NA, full_ratio),
phd_rank = ifelse(phd_rank_miss==1, NA, phd_rank),
max_csstops = ifelse(max_csstops_miss==1, NA, max_csstops)
)
```
<br>
### Part 1: Review outcome regression results {-}
Previously we used outcome regression to estimate the overall effect of tenure clock stopping policies and the effects for males and females. Models from one of our analyses are displayed here.
Review the results of this analysis by interpreting the relevant coefficients and confidence intervals on the odds ratio scale.
```{r}
# Overall ACE
mod_overall <- glm(tenure_policy_school ~ gncs + ns(female_ratio, 3) + ns(full_ratio, 3) + factor(max_csstops)*female, data = tenure, family = "binomial")
summary(mod_overall)
confint(mod_overall)
# Subgroup ACEs
mod_subgroup <- glm(tenure_policy_school ~ gncs*female + ns(female_ratio, 3) + ns(full_ratio, 3) + factor(max_csstops)*female, data = tenure, family = "binomial")
summary(mod_subgroup)
confint(mod_subgroup)
```
<br><br>
### Part 2: Propensity score modeling {-}
We will assume that the `revenue`, `female_ratio`, `full_ratio`, and `max_csstops` are sufficient (proxies for) variables needed to achieve conditional exchangeability. (`female` was included in the outcome regresion models for assessing subgroup effects but wasn't determined to be part of a noncausal path in our causal graph analysis.)
- Fit an appropriate propensity score model called `ps_mod`. Use visualizations to inform the construction of your model. (Just focus on appropriately handling nonlinearity.)
- Use your model to compute appropriate weights, and add these weights to the dataset.
- Use "before and after" visualizations to compare balance of key variables before and after weighting. You'll need to use `factor(gncs)` to make your visualizations--this forces R to view `gncs` as a categorical variable.
- Note: If any distributions look dissimilar after weighting, this is an indication that the propensity score model may be misspecified. We won't address this in our analysis today for time reasons, but in practice, more flexible modeling techniques would be explored.
<br>
Coding notes:
The code below will be useful for exploring nonlinearity in covariate-treatment relationships. The blue smooth reflects observed data trends, and the red smooth shows the predictions from a logistic regression model with a specific model formula.
- Formula `y ~ x`: Covariate is included as a linear term
- Formula `y ~ ns(x,3)`: Covariate is modeled with a flexible nonlinear function (a "spline" with 3 degrees of freedom--don't worry about the details of this)
The model predictions should generally align with observed data trends. If they don't line up well using formula `y ~ x`, update the formula to `y ~ ns(x, 3)` to see if there is an improvement.
```{r}
ggplot(DATA, aes(x = COVARIATE, y = OUTCOME)) +
geom_point() +
geom_smooth(se = FALSE, color = "blue", method = "loess") +
geom_smooth(formula = INSERT_MODEL_FORMULA, method = "glm",
method.args = list(family="binomial"),
se = FALSE, color = "red"
)
```
<br>
The (incomplete) code below is useful for computing appropriate weights:
```{r}
DATA <- DATA %>%
mutate(
ps = predict(ps_mod, newdata = DATA, type = "response"),
ipw = case_when(
TREAT_VAR==1 ~ ???,
TREAT_VAR==0 ~ ???
)
)
```
<br>
You can incorporate weights into most `ggplot2` figures by adding `weight` to the `aes`thetics:
```{r}
aes(..., weight = ip_weight)
```
<br><br>
### Part 3: Modeling with IP weights {-}
We can use the IP weights constructed above in a weighted regression model to estimate causal effects.
- Note: The `svydesign()` and `svyglm()` functions from the `survey` package are used to ensure that the weights are treated appropriately. (An updated ("robust") standard error calculation method is used to appropriately account for the fact that the IP weights are estimated and have uncertainty.)
- Note that we are removing cases with missing values for the weights. (These cases had missing values in the covariates used in the propensity score model.) You'll think about whether this biases the analysis in Part 4.
- The code for estimating the overall ACE is complete. Adapt this code to estimate subgroup ACEs in males and females.
- Using both the confidence intervals and effect magnitudes, discuss the results of your analysis in a contextually meaningful way. (Tie these results back to the research questions.)
- How do confidence interval widths compare to those from outcome regression?
```{r}
# Remove cases with missing weights
tenure_subs <- tenure %>%
filter(!is.na(ipw))
# Set up information about weights
design <- svydesign(ids = ~0, weights = tenure_subs$ipw, data = tenure_subs)
# Fit a marginal structural model to estimate overall ACE
overall_fit <- svyglm(
tenure_policy_school ~ gncs,
data = tenure_subs,
design = design,
family = "quasibinomial"
)
summary(overall_fit)
confint(overall_fit)
# Adapt the svyglm() code above to include an interaction between treatment and female
subgroup_fit <-
```
<br><br>
### Part 4: Causal graphs and missing data {-}
Challenge! How might we use causal graphs to analyze the impact of dropping missing data (a complete case analysis)? Can we use ideas from representing selection bias with causal graphs?