updating explanation of part I to address reasons for dropped observations
[stats_class_2020.git] / psets / pset8-worked_solution.rmd
1 ---
2 title: "Problem set 8: Worked solutions"
3 subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
4 author: "Aaron Shaw"
5 date: "November 23, 2020"
6 output:
7   html_document:
8     toc: yes
9     toc_depth: 3
10     toc_float:
11       collapsed: true
12       smooth_scroll: true
13     theme: readable
14   pdf_document:
15     toc: yes
16     toc_depth: '3'
17   header-includes:
18   - \newcommand{\lt}{<}
19   - \newcommand{\gt}{>}
20 ---
21
22 ```{r setup, include=FALSE}
23 knitr::opts_chunk$set(echo = TRUE, message=FALSE, tidy='styler')
24 ```
25
26 # Import libraries
27
28 I'll start by loading some useful libraries...
29 ```{r dependencies}
30 library(openintro)
31 library(tidyverse)
32 library(ggfortify)
33 library(haven)
34 ```
35
36 # Part I: Mario kart replication/extension  
37
38 ## Import data and setup
39 ```{r}
40 data(mariokart)
41
42 mariokart
43 ```
44
45 To make things a bit easier to manage, I'll select the variables I want to use in the analysis and do some cleanup. Note that I convert the `cond_new` and `stock_photo` variables to logical `TRUE/FALSE` values first (using boolean comparisons) and then coerce them to be numeric values (using `as.numeric()`). This results in 1/0 values corresponding to the observations shown/described in Figure 9.13 and 9.14 on p. 365 of the textbook.
46
47 ```{r}
48 mariokart <- mariokart %>% 
49   select(
50     price = total_pr, cond_new = cond, stock_photo, duration, wheels
51   ) %>% mutate(
52     cond_new = as.numeric(cond_new == "new"),
53     stock_photo = as.numeric(stock_photo == "yes")
54   ) 
55
56 mariokart
57 ```
58
59 Now let's summarize and explore the variables in our model. I'll calculate summary statistics for each, univariate density plots for the continuous measures, boxplots for the dichotomous indicators, and some bivariate plots with each predictor and the outcome variable.
60
61 ```{r}
62 summary(mariokart)
63
64 sd(mariokart$price)
65 sd(mariokart$duration)
66 sd(mariokart$wheels)
67
68 qplot(data=mariokart, x=price, geom="density")
69 qplot(data=mariokart, y=price, geom="boxplot") # check out the outliers!
70 qplot(data=mariokart, x=duration, geom="density")
71 qplot(data=mariokart, x=wheels, geom="density")
72
73 ggplot(data=mariokart, aes(cond_new, price)) + geom_boxplot()
74 ggplot(data=mariokart, aes(stock_photo, price)) + geom_boxplot()
75 ggplot(data=mariokart, aes(duration, price)) + geom_point()
76 ggplot(data=mariokart, aes(wheels, price)) + geom_point()
77
78 ```
79
80 I'm also going to calculate correlation coefficients for all of the variables. We can discuss what to make of this in class:
81
82 ```{r}
83 cor(mariokart)
84 ```
85
86 ## Replicate model results
87
88 Based on the information in the textbook, I'm assuming that the model is an ordinary least squares regression.
89 ```{r}
90 model <- lm(price ~ cond_new + stock_photo + duration + wheels, data=mariokart)
91 summary(model)
92 ```
93
94 Huh, that doesn't quite look like what's in the textbook Figure 9.15...
95
96 I'll go ahead and calculate a confidence interval around the only parameter for which the model rejects the null hypothesis (`wheels`):
97 ```{r}
98 confint(model, "wheels")
99 ```
100
101 Without more information from the authors of the textbook it's very hard to determine exactly where or why the differences emerge. My guess at this point is that it might have something to do with those outlying `price` values (take a look at that boxplot and density plot again). Maybe they did something to transform the variable? Remove the outliers? I have no idea.
102
103 ## Assess model fit/assumptions  
104
105 I've already generated a bunch of univariate and bivariate summaries and plots. Let's inspect the residuals more closely to see what else we can learn. I'll use the `autoplot()` function from the `ggfortify` package.
106
107 ```{r}
108 autoplot(model)
109 ```
110
111 There are those outliers again. At this point, there are a number of issues with this model fit that I'd want to mention/consider:
112 * The distribution of the dependent variable (`price`) is skewed, with two extreme outliers. I'd recommend trying some transformations to see if it would look more appropriate for a linear regression and/or inspecting the cases that produced the outlying values more closely to understand what's happening there and identify reasons why they're so different.
113 * The plots of the residuals reveal those same two outlying points are also outliers with respect to the line of best fit. That said, they are not exerting huge amounts of leverage on the estimates, so it's possible that the estimates from the fitted model wouldn't change *too much* without those two points. Indeed, based on the degrees of freedom reported in Figure 9.15 (136) vs. the number reported in our version of the model (138) my best guess at this point is that the textbook authors silently dropped those two outlying observations from their model.
114
115 More out of curiosity than anything, I'll create a version of the model that drops the two largest values of price. From the plots, I can see that those two are the only ones above \$100, so I'll use that information here:
116
117 ```{r}
118 summary(
119   lm(price ~ cond_new + stock_photo + duration + wheels, 
120      data = mariokart[mariokart$price < 100,]
121      )
122   )
123 ```
124
125 What do you know. That was it. The difference in $R^2$ is huge! 
126
127 A little further digging (by Nick Vincent) revealed that these two outliers come from auctions where the Mario kart game was being sold as part of a bundle along with other games. You can look this up in the `title` field from the original dataset using the following block of code:
128
129 ```{r}
130 data(mariokart)
131
132 mariokart %>% 
133   filter(total_pr > 100) %>% 
134   select(id, total_pr, title)
135
136 ```
137
138
139 What do you make of the textbook authors' decision to drop the observations? Can you make a case for/against doing so? What seems like the right decision and the best way to handle this kind of situation?
140
141 ## Interpret some results  
142
143 The issues above notwithstanding, we can march ahead and interpret the results of the original model that I fit. Here are some general comments and some specifically focused on the `cond_new` and `stock_photo` variables:
144 * Overall, the linear model regressing total auction price on condition, stock photo, duration, and number of Wii wheels shows evidence of a positive, significant relationship between number of wheels and price. According to this model fit, an increase of 1 wheel in a listing is associated with a total auction price increase of \$10 on average (95% confindence interval: \$4.57-\$15.32). 
145 * The point estimate for selling a new condition game is positive, but with a large standard error. The model fails to reject the null of no association and provides no evidence of any relationship between the game condition and auction price. 
146 * The point estimate for including a stock photo is negative, but again, the standard error is very large and the model fails to reject the null hypothesis. There is no evidence of any relationship between including a stock photo and the final auction price.
147
148 ## Recommendations  
149
150 Based on this model result, I'd recommend the prospective vendor of a **used** copy of the game not worry about it too much unless they can get their hands on some extra Wii wheels, since it seems like the number of Wii wheels explains variations in the auction price outcome more than anything else they can control (such as whether or not a stock photo is included). 
151
152 # Part II: Hypothetical study 
153
154 ## Import, explore, summarize  
155
156 I'll start off by just importing things and summarizing the different variables we care about here:
157 ```{r}
158 grads <- readRDS(url("https://communitydata.science/~ads/teaching/2020/stats/data/week_11/grads.rds"))
159
160 summary(grads)
161
162 table(grads$cohort)
163 sd(grads$gpa)
164 sd(grads$income)
165
166 qplot(data=grads, x=gpa, geom="density")
167 qplot(data=grads, x=income, geom="density")
168
169 ggplot(data=grads, aes(x=gpa, y=income)) + geom_point()
170 ```
171
172 I'll also calculate some summary statistics and visual comparisons within districts (cohorts):
173
174 ```{r}
175
176 tapply(grads$income, grads$cohort, summary)
177 tapply(grads$income, grads$cohort, sd)
178
179 tapply(grads$gpa, grads$cohort, summary)
180 tapply(grads$gpa, grads$cohort, sd)
181 ```
182
183 Note that you could also do this pretty easily with a call to `group_by` and `summarize` in the tidyverse:
184
185 ```{r}
186 grads %>% 
187   group_by(cohort) %>% 
188   summarize(
189     n = n(),
190     min = min(income),
191     mean = mean(income),
192     max = max(income),
193     sd = sd(income)
194   ) 
195
196 grads %>% 
197   group_by(cohort) %>% 
198   summarize(
199     n = n(),
200     min = min(gpa),
201     mean = mean(gpa),
202     max = max(gpa),
203     sd = sd(gpa)
204   )
205 ```
206
207
208 Huh. Those are remarkably similar values for the group means and the group standard deviations...weird.
209
210 Onwards to plotting:
211 ```{r}
212 ggplot(data=grads, aes(x=cohort, y=income)) + geom_boxplot()
213
214 ggplot(data=grads, aes(x=gpa, y=income, color=cohort)) + geom_point()
215 ```
216
217 Those plots are also a little strange. Even though this is just a simulated analysis, it still seems weird that overall the scatterplot just looks like a big mass of points, but when I color the points by district, I can see some regularities within groups. At this point, I might want to facet the scatterplots by district to see any patterns more clearly.
218
219 ```{r}
220 ggplot(data=grads, aes(x=gpa, y=income, color=cohort)) + geom_point() + facet_wrap(vars(cohort))
221 ```
222
223 Okay, that's...a joke (in particular, cohort 8 looks like a sideways dinosaur). At this point, if I were really working as a consultant on this project, I would write to the client and start asking some probing questions about data quality (who collected this data? how did it get recorded/stored/etc.? what quality controls were in place?). I would also feel obligated to tell them that I suspect there's just no way the data correspond to the variables they think are here. If you did that and the client was honest, they might tell you [where the data actually came from](https://www.autodesk.com/research/publications/same-stats-different-graphs).
224
225 In the event that you marched ahead with the analysis and are curious about what that could have looked like, I've provided some example code below. That said, **this is a situation where the assumptions and conditions necessary to identify ANOVA, t-tests, or regression are all pretty broken** because the data was generated programmatically in ways that undermine the kinds of interpretation you've been asked to make. The best response here (IMHO) is to abandon these kinds of analysis once you discover that there's something systematically weird going on like this. While the experience of discovering a scatterplot dinosaur in your data is...unlikely outside of the context of a problem set, there are many situations in which careful data exploration will bring a realization that you just don't understand some important things about the sources or qualities of your data. You have to learn to identify these moments and develop strategies for dealing with them! Often, the statistical procedures will "work" in the sense that they will return a result without any apparent errors, but because those results aren't even close to meaningful, any relationships you do observe in the data reflect something different than the sorts of relationships the statistical procedures were designed to identify.
226
227 ## Fake analysis for fake data  
228
229 Okay, if you wanted example code to look at for this, here it is. Please just keep in mind that the results are not informative!
230 ```{r}
231 summary(aov(income ~ cohort, data=grads))  # no global differences of means across groups
232
233 summary(grads.model <-  lm(income ~ gpa + cohort, data=grads)) # gpa percentile has a small, negative association with income. 
234
235 confint(grads.model, "gpa") # 95% confidence interval
236
237 ```
238
239 Note that the failure to reject the null of any association between district and income in the ANOVA would not (even in the event of more realistic data) provide very compelling evidence that the relationship between GPA and income does not vary by cohort. There were several things you might have done here. One is to calculate correlation coefficients within groups. Here's some tidyverse code that does that:
240
241 ```{r}
242 grads %>% 
243   group_by(cohort) %>% 
244   summarize(
245     correlation = cor(income, gpa)
246   )
247 ```
248 Because these correlation coefficients are nearly identical, I would likely end an analysis here and conclude that the correlation between gpa and income appears to be consistently small and negative. If you wanted to go further, you could theoretically calculate an interaction term in the model (by including `I(gpa*cohort)` in the model formula), but the analysis up to this point gives no indication that you'd be likely to find much of anything (and we haven't really talked about interactions yet). On top of that, there's a literal dinosaur lurking in your data...just give up!  
249
250 # Part III: Trick or treating again  
251 ## Import and update data
252
253 Revisit the text and worked solutions for problem set 7 for more details about the study design, data collection and more.
254
255 ```{r import}
256 ## reminder that the "read_dta()" function requires the "haven" library
257
258 df <- read_dta(url('https://communitydata.science/~ads/teaching/2020/stats/data/week_07/Halloween2012-2014-2015_PLOS.dta'))
259
260 df <- df %>%
261   mutate(
262     obama = as.logical(obama),
263     fruit = as.logical(fruit),
264     year = as.factor(year),
265     male = as.logical(male),
266     age = age-mean(age, na.rm=T)
267   )
268
269 df
270 ```
271
272 That looks consistent with what we want here. Let's fit and summarize the model:  
273 ```{r}
274 f <- formula(fruit ~ obama + age + male + year)
275
276 fit <- glm(f, data=df, family=binomial("logit"))
277
278 summary(fit)
279 ```
280
281 Interesting. Looks like adjusting for these other variables in a regression setting allows us to uncover some different the results. 
282
283 Onwards to generating more interpretable results. You might recall that the big problem with interpreting logistic regression is that the results are given to you in "log-odds." Not only is it difficult to have intuitions about odds, but intuitions about the natural logarithms of odds are just intractable (for most of us). 
284
285 To make things easier, the typical first step is to calculare odds-ratios instead of log-odds. This is done by exponentiating the coefficients (as well as the corresponding 95\% confidence intervals):
286
287 ```{r}
288 ## Odds ratios (exponentiated log-odds!)
289 exp(coef(fit))
290 exp(confint(fit))
291 ```
292
293 You can use these to construct statements about the change in odds of the dependent variable flipping from 0 to 1 (or `FALSE` to `TRUE`) predicted by a 1-unit change in the corresponding predictor (where an odds ratio of 1 corresponds to unchanged odds). We'll interpret the `obamaTRUE` odds ratio below.
294
295 Now, model-predicted probabilities for prototypical observations. Recall that it's necessary to create synthetic ("fake"), hypothetical individuals to generate predicted probabilities like these. In this case, I'll create two versions of each fake kid: one assigned to the treatment condition and one assigned to the control. Then I'll use the `predict()` function to generate fitted values for each of the fake kids.
296
297 ```{r}
298 fake.kids = data.frame(
299   obama = c(FALSE, FALSE, TRUE, TRUE),
300   year = factor(rep(c("2015", "2012"), 2)),
301   age = rep(c(9, 7), 2),
302   male= rep(c(FALSE, TRUE), 2)
303 )
304
305 fake.kids.pred <- cbind(fake.kids, pred.prob = predict(fit, fake.kids, type="response"))
306
307 fake.kids.pred
308 ```
309
310 Note that [this UCLA logit regression tutorial](https://stats.idre.ucla.edu/r/dae/logit-regression/) also contains example code to help extract standard errors and confidence intervals around these predicted probabilities. You were not asked to create them here, but if you'd like an example here you go. The workhorse is the `plogis()` function, which essentially does the inverse logit transformation detailed in the textbook chapter 8:
311
312 ```{r}
313 fake.kids.more.pred <- cbind(fake.kids,
314                       predict(fit, fake.kids, type= "link", se=TRUE))
315
316 within(fake.kids.more.pred, {
317   pred.prob = plogis(fit);
318   lower = plogis(fit - (1.96*se.fit));
319   upper = plogis(fit + (1.96*se.fit))
320 } )
321 ```
322
323 ## Sub-group analysis
324
325 To do this, we'll need to create a slightly new model formula that drops the `year` term since we're going to restrict each subset of the data along that dimension.
326
327 Once I fit the models, I'll use the `stargazer` package to create a reasonably pretty regression table that incorporates all three summaries.  
328 ```{r}
329 f2 <- formula(fruit ~ obama + age + male)
330
331 m2012 <- glm(f2, data=df[df$year == "2012",], family=binomial("logit"))
332 m2014 <- glm(f2, data=df[df$year == "2014",], family=binomial("logit"))
333 m2015 <- glm(f2, data=df[df$year == "2015",], family=binomial("logit"))
334
335  
336 ## I can make a pretty table out of that:
337 library(stargazer)  
338 stargazer(m2012, m2014, m2015, column.labels = c("2012", "2014", "2015"), type="text")
339
340 ```
341
342 ## Interpret and discuss
343
344 Well, for starters, the model providing a "pooled" estimate of treatment effects while adjusting for age, gender, and study year suggests that the point estimate is "marginally" statistically significant ($p <0.1$) indicating some evidence that the data support the alternative hypothesis (being shown a picture of Michelle Obama causes trick-or-treaters to be more likely to pick up fruit than the control condition). In more concrete terms, the trick-or-treaters shown the Obama picture were, on average, about 26\% more likely to pick up fruit than those exposed to the control (95\% CI:  $-4\%~-~+66\%$).[^1] In even more concrete terms, the estimated probability that a 9 year-old girl in 2015 and a 7 year-old boy in 2012 would take fruit increase about 17\% and 19\% respectively on average (from 29\% to 34\% in the case of the 9 year-old and from 21\% to 25\% in the case of the 7 year-old). These findings are sort of remarkable given the simplicity of the intervention and the fairly strong norm that Halloween is all about candy.
345
346 [^1]: Remember when I said we would use those odds ratios to interpret the parameter on `obamaTRUE`? Here we are. The parameter value is approximately 1.26, which means that the odds of picking fruit are, on average, 1.26 times as large for a trick-or-treater exposed to the picture of Michelle Obama versus a trick-or-treater in the control condition. In other words, the odds go up by about 26\% ($= \frac{1.26-1}{1}$).
347
348 All of that said, the t-test results from Problem set 5 and the "unpooled" results reported in the sub-group analysis point to some potential concerns and limitations. For starters, the fact that the experiment was run iteratively over multiple years and that the sample size grew each year raises some concerns that the study design may not have anticipated the small effect sizes eventually observed and/or was adapted on the fly. This would undermine confidence in some of the test statistics and procedures. Furthermore, because the experiment occurred in sequential years, there's a very real possibility that the significance of a picture of Michelle Obama shifted during that time period and/or the house in question developed a reputation for being "that weird place where they show you pictures of Michelle Obama and offer you fruit." Whatever the case, my confidence in the findings here is not so great and I have some lingering suspicions that the results might not replicate.
349
350 On a more nuanced/advanced statistical note, I also have some concerns about the standard errors. This goes beyond the content of our course, but basically, a randomized controlled trial introduces clustering into the data by-design (you can think of it as analogous to the observations coming from the treatment "cluster" and the control "cluster"). In this regard, the normal standard error formulas can be biased. Luckily, there's a fix for this: compute "robust" standard errors as a result and re-calculate the corresponding confidence intervals. Indeed, robust standard errors are often considered to be the best choice even when you don't know about potential latent clustering or heteroskedastic error structures in your data. [This a short pdf](https://oes.gsa.gov/assets/files/calculating-standard-errors-guidance.pdf) provides a little more explanation, citations, as well as example R code for how you might calculate robust standard errors.

Community Data Science Collective || Want to submit a patch?