]> code.communitydata.science - stats_class_2020.git/blob - psets/pset5-worked_solution.rmd
initial commit. missing text still
[stats_class_2020.git] / psets / pset5-worked_solution.rmd
1 ---
2 title: "Problem set 5: Worked solutions"
3 subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
4 author: "Aaron Shaw (with Jeremy Foote)"
5 date: "October 26, 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 # Programming Challenges  
27 ## PC1. Import  
28
29 You had the option of downloading the dataset in several formats, including a Stata 13 ".dta" file. For demonstration purposes, I'll use that proprietary format here. I also uploaded a copy of the dataset to the course website for pedagogical purposes (so you can replicate my code here on your own machine). 
30
31 There are a few ways to import data from Stata 13 files. One involves using the appropriately named `readstata13` package. Another uses the `haven` package (more of a general purpose tool for importing data from binary/proprietary formats). I'll go with the `read_dta()` command from `haven` here because it works better with importing from a URL. 
32
33 ```{r import}
34 library(haven)
35
36 df <- read_dta(url('https://communitydata.science/~ads/teaching/2020/stats/data/week_07/Halloween2012-2014-2015_PLOS.dta'))
37
38 ```
39 ## PC2. Explore and cleanup  
40
41 ```{r explore}
42 head(df)
43 summary(df)
44 ```
45
46 There are a few noteable things about the dataset. One is the `neob`, which the codebook says means "not equal to obama"; in other words, it's the converse of the `obama` column. The `treat_year` column is a unique index of the `obama` column and the `year` column. See the codebook for more information. Happily, we only need to use the first two columns for now.
47
48 I'll drop everything else and convert those two columns into logical vectors here using a couple of tidyverse `dplyr` commands:
49
50 ```{r recodes}
51 library(tidyverse)
52
53 df <- df %>%
54   select(obama, fruit) %>%
55   mutate(
56     obama = as.logical(obama),
57     fruit = as.logical(fruit)
58   )
59
60 head(df)
61
62 ```
63
64 ## PC3. Summarize key variables  
65 I'll run summary again and then make my contingency table.
66
67 ```{r summarize}
68 summary(df)
69
70 obama.tbl <- table(took_fruit=df$fruit, saw_flotus=df$obama)
71 obama.tbl
72
73 obama.prop.tbl <- proportions(table(took_fruit=df$fruit, saw_flotus=df$obama), margin=2)
74 obama.prop.tbl
75 ```
76
77 ## PC4. Test for differences between groups
78 So, the test we want to conduct here focuses on the difference between the proportion of the two groups (those shown a picture of Michelle Obama vs those not shown a picture of Michelle Obama) who took fruit vs. candy. Labeling the difference in proportions as $\Delta_{fruit}$, the comparison can be constructed around the following (two-sided) hypothesis test
79
80 $$ H_0:~\Delta_{fruit} = 0$$
81 $$ H_A:~\Delta_{fruit} \neq 0$$
82 As discussed at length in *OpenIntro* Chapter 6, a great way to determine if two groups are independent (in terms of proportions or counts) is a $\chi^2$ test. 
83
84 Are the conditions for a valid test met? It seems so, since there are many observations in each cell of the table being compared and the observations appear to have been collected in a way that ensures independence. 
85
86 The $\chi^2$ test is implemented as `chisq.test()` in R. Since it's a 2x2 comparison, we can also test for a difference in proportions using the `prop.test()` function. Let's do both.
87
88 ```{r hypothesis tests}
89 chisq.test(obama.tbl)
90
91 prop.test(obama.tbl)
92 ```
93
94 Notice that both functions report identical $\chi^2$ test results and p-values. Did you expect this based on the *OpenIntro* reading?
95
96 Recall that if you want to double-check that p-value, you can also calculate it "by-hand" using the `pchisq()` function:
97 ```{r p-value by hand}
98 pchisq(1.8637, df=1, lower.tail=FALSE)
99 ```
100
101 We'll discuss the interpretation of these results in our class session this week.
102
103 ## PC5. Replicate a figure  
104
105 In order to replicate the top panel of Figure 1, we'll first want to calculate the proportion and standard error for fruit-takers in the treatment and control groups. These can be calculated individually or using a function (guess which one we'll document here). Also, note that I'm going to use the `complete.cases()` function to eliminate the missing items for the sake of simplicity.
106
107 ```{r calc proportions and se}
108 df <- df[complete.cases(df),]
109
110 prop.se = function(values){# Takes in a vector of T/F values
111   N = length(values)
112   prop = mean(as.numeric(values))
113   se = sqrt(prop * (1-prop)/N)  ## textbook formula for SE of a proportion
114   return(c(prop, se))
115 }
116
117 prop.se(df$fruit[df$obama])
118 prop.se(df$fruit[!df$obama])
119
120
121 ```
122 In order to graph that it will help to convert the results into a data frame with clearly-labeled variable names and values:
123 ```{r bind up proportions and se}
124 prop.and.se <- data.frame(
125   rbind(
126     prop.se(df$fruit[df$obama]),
127     prop.se(df$fruit[!df$obama])
128     )
129   )
130
131 names(prop.and.se) <- c("proportion", "se")
132 prop.and.se$obama <- c(TRUE, FALSE)
133
134 prop.and.se
135 ```
136 Now we can start to build a visualization
137 ```{r base plot}
138 ## library(ggplot2)  ## already imported w the tidyverse
139
140 p <- ggplot(prop.and.se, aes(x=obama,y=proportion)) + 
141   geom_point(aes(color=obama), size=5)
142
143 p
144
145 ```
146
147 Looking pretty good. Let's go ahead and add some error bars:
148 ```{r err bars}
149 p + geom_errorbar(aes(ymin=proportion - 1.96 * se, # Add error bars
150                     ymax=proportion + 1.96 * se))
151 ```
152
153 Now let's clean up those error bars a bit...
154 ```{r better error bars}
155 p1 <- p + geom_errorbar(aes(ymin=proportion - 1.96 * se,
156                     ymax=proportion + 1.96 * se,
157                     width=0, # Remove the whiskers
158                     color=obama),size=1.1) # Make bars thicker and apply color
159
160 p1
161 ```
162
163 Great! Now I can style it a bit more by flipping the plot on it's side with `coord_flip()`, adding a theme, converting it to grayscale, and fixing up my axes a bit:
164
165 ```{r style plot}
166 p1 + coord_flip() + # Flip the chart
167   theme_light() + # Change the theme (theme_minimal is also nice)
168   scale_color_manual(values = c('gray','black'), guide=F) + # Change the colors
169   ylim(0,.5) + # Change the y axis to go from 0 to .5
170   ylab('Proportion choosing fruit') + # Add labels
171   xlab('Picture shown was Michelle Obama')
172
173 ```
174
175 Pretty good!  
176
177 ## PC6. Export a table  
178
179 Here's one way to export our table, using `write.csv()` (and have commented it out, so you can run locally and determine where the exported file goes).
180 ```{r simple csv export}
181 ## uncomment to generate exported file
182 ## write.csv(obama.tbl, file = 'crosstabs.csv')
183 ```
184 We can make sure it worked by importing it (again, commented out here):
185 ```{r test simple export}
186 ## read.csv('crosstabs.csv')
187 ```
188 We lost some information, because the `table()` function doesn't save column names. Another way to do this would be to change it into a dataframe first, like this:
189 ```{r convert to dataframe}
190 data.frame(obama.tbl)
191 ```
192 and then save that dataframe. Note that you can drop the rownames when importing (or exporting)
193 ```{r export dataframe}
194 ## write.csv(data.frame(obama.tbl), file = 'crosstabs.csv', row.names = FALSE)
195
196 ## read.csv('crosstabs.csv')
197 ```
198 That's formatted a little bit funny, but it's still usable. 
199
200 You could also use the `xtable` package to do this. The package has many functions to customize table outputs in several formats. A relatively simple way to generate an html table looks like this:
201 ```{r xtable example}
202 library(xtable)
203 print(xtable(obama.tbl), type="html")
204 ```
205
206 You can export by assigning the html to an object and saving. I've commented it out here so you can choose whether/where to create the file:
207 ```{r export xtable}
208 ## uncomment to generate file output
209
210 ## print(xtable(obama.tbl), type="html", file="example_table.html")
211 ```
212
213 You should be able to open that file in a web-browser. 
214
215 There is a lot of documentation and examples online to help you customize as you see fit. If you're really excited about exporting tables, you might also take a look at the `tables` package, which has some nice export options.
216
217 # Empirical paper questions  
218 ## EQ1. LilyPad Arduino users  
219
220 a) The unit of analysis is the customer. The dependent variable is the type of board purchased and the independent variable is gender. Males, females, and unknown gender customers are being compared. This is a two-way test.  
221
222 b) For this type of comparison statistical tests help to give (or take away) confidence in any observed differences in counts or proportions across categories. Choosing a statistical test is based on the question that you want to answer and the type of data that you have available to answer it. For example, if this were continuous numeric data (e.g., the amount of money spent on electronics for men and women) then we would want a different to compare those distributions. Given that the test compares counts (or proportions) a $\chi^2$ test for independence is appropriate. 
223
224 c) The null hypothesis ($H_0$) is that the board purchased is independent of the gender of the customer. The alternative hypothesis ($H_A$) is that board purchase choice is dependent on gender.  
225
226 d) A $\chi^2$ test results indicate that board purchase behavior differs by gender ($p~\leq~0.05$). This difference is convincing, but it does directly not tell us what the authors set out to understand, which is the difference between men and women (the test as-implemented might have identified a significant difference in the number of unknown gender customers across board types!). Many of these concerns are addressed in the text and with additional tests, giving increased confidence in the observed differences.  
227
228
229 ## EQ2. Two blogospheres  
230
231 a) The data are counts for two categorical variables and the procedure used was a $\chi^2$ test. The null hypothesis is that blog governance (by one person or more than one person) is independent of whether the blog was on the left or the right ideologically.  
232
233 b) One way to approach this focuses on the credibility of the null hypothesis. A null of equality/independence in this case seems like it might be a bit of a stretch since there are potentially so many factors involved in determining site governance. If we accept that the null hypothesis of no difference across the two groups is compelling, it seems like it could be surprising to see these results in a world where ideological orientation and blog governance have no relationship. In this respect, it makes sense to believe that there is some relationship of dependence. A closer reading of the paper suggests a different reason to be skeptical: the way that the measure of blog governance groups the data into categories. The authors could have grouped them differently (e.g., 1-2 people, 3-4 people, and 5+ people). If the decision on how to group was made after seeing the data or if the observed result depends on the choice of grouping, then we have good reason to be skeptical.  
234
235 c) We can do this in R.
236
237 ```{r replicate Shaw and Benkler tabular}
238
239 ## First we create the dataframe
240 df = data.frame(Governance=c('Individual','Multiple', 'Individual','Multiple'),
241                 Ideology = c('Left','Left','Right','Right'),
242                 Count = c(13,51,27,38))
243
244 df
245
246 ## We can make sure it's the same by testing the Chi-squared                
247 chisq.test(matrix(df$Count, nrow=2))
248
249 ## We can convert that into proportions several ways. 
250 proportions(matrix(df$Count, nrow=2), margin =2)
251
252 ## Here's one using tidyverse code that yields more readable results:
253 percentage_data <- df %>% 
254   group_by(Ideology) %>% 
255   summarize(individual_ratio = sum(Count[Governance=='Individual']) / sum(Count),
256             ideology_count = sum(Count))
257
258 percentage_data
259 ```
260
261 And here we go with a figure using the `geom_bar()` layer in `ggplot2`. Note that by calling `stat=\'identity\'', we tell ggplot to vizualize the provided counts (instead of trying to summarize the data further):
262 ```{r replicate Shaw and Benkler figure 2}
263 shaw_benkler_plot = percentage_data %>%
264   ggplot(aes(x=Ideology, y=individual_ratio * 100)) + 
265     geom_bar(stat='identity', aes(fill=c('red','blue')), show.legend=F) + 
266     ylab('Percentage of Blogs') + theme_minimal()
267
268 shaw_benkler_plot
269 ```
270
271 If we want to add error bars, we need to calculate them (Note that ggplot could do this for us if we had raw data - what a helpful reminder to always share your data!). 
272
273 For our purposes here, you might decide to use confidence intervals or standard errors (both seem like reasonable choices as long as you label them). Either way, ggplot has a `geom_errorbar` layer that is very useful.
274
275 Remember that for a binomial distribution (we can consider individual/non-individual as binomial), confidence intervals are
276 $\hat{p} \pm z^* \sqrt{\frac{\hat{p}~(1-\hat{p})}{n}}$ and we can approximate $z^*~=~1.96$ for a 95% confidence interval.
277
278 ```{r add CIs to Shaw and Benkler Fig 2}
279 ci_95 = 1.96 * sqrt(percentage_data$individual_ratio * (1 - percentage_data$individual_ratio)/percentage_data$ideology_count)
280
281 shaw_benkler_plot + geom_errorbar(aes(ymin=(individual_ratio-ci_95)*100, ymax=(individual_ratio + ci_95)*100),
282                                   alpha = .3, 
283                                   size=1.1, 
284                                   width=.4)
285
286 ```
287
288 The 95% confidence intervals overlap in this case, indicating that the true population proportions may not be as far apart as our point estimates suggest. Note that this is not the same as the hypothesis test (illustrating one of Reinhart's points).  
289
290 d) On the one hand, we don't need to worry about the base rate fallacy because the sizes of both groups are about the same and the paper does not abuse the evidence too egregiously. The base rate fallacy would likely come into play, however, in the ways that the results are (mis)represented. For example, you might imagine some news coverage looking at these results and claiming something (totally wrong!) like "study finds right wing blogs more than twice as likely to be solo affairs." This is taking a relationship between the sample proportions ($\hat{p}$ in the language of our textbook) and converting that into a statement about the relationship between population proportions ($p$). That would be a bit of a mess (and absolutely characterizes the news coverage that did follow the study). 
291
292     Another way in which the base rate fallacy could play a role in this paper, however, concerns the presence of multiple comparisons. The authors conducted numerous statistical tests (indeed, one of the authors seems to recall that some of the tests were not even reported in the paper <gasp!>) and they make no effort to address the baseline probability of false positives. 
293
294     In any case, the point here is that the statistical tests reported in the paper may not mean exactly what the authors said they did in the context of the publication. That may or may not change the validity of the results, but it should inspire us all to do better statistical analysis.

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