]> code.communitydata.science - stats_class_2020.git/blob - psets/pset4-worked_solution.rmd
62b88f9390e23b57d5f35dca063f512805455850
[stats_class_2020.git] / psets / pset4-worked_solution.rmd
1 ---
2 title: "Problem set 4: Worked solutions"
3 subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
4 author: "Aaron Shaw"
5 date: "October 19, 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 Load the dataset. As usual, you'll need to edit this block of code to point at the appropriate file location to make it work on your machine.
29 ```{r}
30 pop <- read.delim(file=url("https://communitydata.science/~ads/teaching/2020/stats/data/week_06/population.tsv"))
31
32 # Same thing using `read_tsv()`. Notice that it can handle URLs directly.
33 library(tidyverse)
34
35 pop <- read_tsv("https://communitydata.science/~ads/teaching/2020/stats/data/week_06/population.tsv")
36
37 head(pop)
38
39 ```
40
41 While we're at it, let's import the data from Problem Set 2:
42
43 ```{r}
44 ## Insert an appropriate group number as needed:
45 ps2 <- read.csv(file=url("https://communitydata.science/~ads/teaching/2020/stats/data/week_04/group_03.csv"), row.names=NULL)
46
47 head(ps2)
48 ```
49
50
51 ## PC2. Compare means  
52 Straightforward enough:  
53
54 ```{r}
55 mean(pop$x, na.rm=T)
56
57 mean(ps2$x, na.rm=T)
58 ``` 
59
60 Given the structure of the dataset for this week, it's also no big deal to calculate all of the group means (and you can compare your answers against these if you like). Notice that I can define and call a function within my call to `tapply()`:
61 ```{r}
62 s.means <- tapply(pop$x, pop$group, function(x){
63                     round(mean(x, na.rm=T), 2)})
64
65 s.means
66 ```
67
68 For anyone who decided to try this using a Tidyverse approach, here's a way you might do that with a call to `group_by()` and `summarize()`:
69 ```{r}
70 pop %>%
71   group_by(group) %>%
72   summarize(
73     mean(x, na.rm=T)
74   )
75 ```
76
77 Knowing that each group was a random sample from the entire population, we might think of the individual group (sample) means as constructing a *sampling distribution of the (population) mean*.
78
79
80 ## PC3. CI of the mean  
81
82 I'll do this two ways. First, I can plug in the values from one group sample into the formula for the standard error and then add/subtract twice the standard error from the mean to find the 95% CI.
83 ```{r}
84 se <- sd(ps2$x, na.rm=T) / sqrt(length(ps2$x[!is.na(ps2$x)]))
85
86 mean(ps2$x, na.rm=T)-(2*se) ## lower
87 mean(ps2$x, na.rm=T)+(2*se) ## upper
88 ```
89
90 Now, I'll write a more general function to calculate confidence intervals. Note that my function here takes an argument for `x` as well as an `alpha` argument with a default value of 0.05. The function then goes on to divide alpha (and it's complementary probability) by 2. Can you explain why this division step is necessary?
91
92 ```{r}
93 ci <- function (x, alpha=0.05) {
94     x <- x[!is.na(x)]
95     probs <- c(alpha/2, 1-alpha/2)
96     critical.values <- qnorm(probs, mean=0, sd=1)
97     se <- sd(x)/sqrt(length(x))
98     return(
99       round( mean(x) + critical.values * se, 2) )
100 }
101
102 ## Here I run the function on the group 3 data
103 ci(ps2$x)
104 ```
105
106 Again, it's possible to use `tapply()` to estimate this for every group:
107
108 ```{r}
109 group.confints <- tapply(pop$x, pop$group, ci)
110 group.confints
111 ```
112
113 Recall that the population mean for `x` was 2.60. Since the group samples are random samples, we should not be surprised that the group means are different from the population mean. We should also not be surprised that the 95% CI for the population mean estimated from at least one of the samples does *not* include the true population mean. Recall that our confidence interval is 95%, so we can expect to be wrong about 1/20 times (on average)! In this case, we got unlucky since 2 of our confidence intervals around the sample means do not include the population mean.  
114
115 ## PC4. Compare distributions  
116
117 I'll start with the single comparison between the population and my Problem Set 2 sample using base-R functions:
118
119 ```{r}
120 hist(pop$x)
121 hist(ps2$x)
122
123 summary(pop$x)
124 sd(pop$x, na.rm=T)
125
126 summary(ps2$x)
127 sd(ps2$x, na.rm=T)
128 ```
129
130 Notice the differences between the shapes of the histograms as well as the summary statistics. In particular, you might consider that the sample has a *larger* standard deviation than the population. Given what you know about this statistic and the conditions under which this data was generated, could it have worked out otherwise (i.e., could the standard deviation of the population have been larger than that of the sample?). How?
131
132
133 Moving right along, here's `tapply()` code to construct the same comparisons for each group. 
134
135 ```{r}
136 tapply(pop$x, pop$group, summary)
137
138 tapply(pop$x, pop$group, sd, na.rm=T)
139
140 ```
141
142 And here's a visual comparison. A ggplot2 "faceted" set of histograms will produce an easy visual comparison (and also makes for concise code). First I'll do it in a slightly simplified way:
143 ```{r}
144 library(ggplot2)
145 ggplot(pop, aes(x)) + geom_histogram() + labs(title="Population distribution of X")
146
147 ggplot(pop, aes(x)) + geom_histogram() + facet_wrap(. ~ group) + labs(title="Conditional distributions of X")
148 ```
149
150 It's possible to see some of the differences there, but it might be helpful to overlay a brightly colored density plot. In the process, I can also convert the histograms themselves into density plots to make it easier to directly contrast each group against the population without having to worry about the divergent y-axis scales. The code below does that by providing a call to `after_stat(density)` inside the initial plotting aesthetics and a `color` argument to `geom_density()`. For more examples like this, you might search the phrase "density histogram" online.
151
152 ```{r}
153 library(ggplot2)
154 ggplot(pop, aes(x, after_stat(density))) + geom_histogram() + geom_density(color="red") + labs(title="Population distribution of X")
155
156 ggplot(pop, aes(x, after_stat(density))) + geom_histogram(aes(x,)) + facet_wrap(. ~ group) + geom_density(color="red") + labs(title="Group distributions of X")
157
158 ```
159
160 The conditional distributions all look a little bit different from each other and from the population distribution. Again, none of this should be shocking given the relationship of the samples to the population.  
161
162 ## PC5. Std. dev. of conditional means  
163
164 I can do this by constructing my vector of group means again and calling `sd()` on that.
165 ```{r}
166 s.means <- tapply(pop$x, pop$group, mean, na.rm=T)
167 s.means
168
169 sd(s.means)
170
171 ## This was the standard error from one of the groups that I calculated earlier:
172 se
173 ```
174 As mentioned earlier, the distribution of sample means drawn from the population is the *sampling distribution*. The standard error of the mean estimated from any of the individual groups/samples should be a good approximation of (but not necessarily equal to!) the standard deviation of the sampling distribution of the means. 
175
176 ## PC6. A simulation  
177
178 Since there's going to be some randomness in the next few commands, I'll set a seed value to ensure the results are reproducible and consistent on other machines.
179
180 ```{r}
181 set.seed(20201019)
182 ```
183
184 **(a) Simulate draws from a known distribution**  
185
186 There are a few ways to do this, but one of the most intuitive is to define my vector of possible values and then sample it repeatedly *with replacement*.
187 ```{r}
188 my.vals <- seq(0, 9)
189 pop.unif <- sample(my.vals, 10000, replace=TRUE)
190 ```
191
192 **(b) Take the mean**  
193
194 ```{r}
195 mean(pop.unif)
196 hist(pop.unif)
197 ```
198
199 **(c) Draw samples and describe them**  
200
201 Again, many ways to go about this. A good first step might be to do the sampling part once:
202
203 ```{r}
204 sample(pop.unif, 2, replace=T)
205 ```
206 Now I can call the mean of that:
207
208 ```{r}
209 mean(sample(pop.unif, 2, replace=T))
210 ```
211 Now, I want to run that 100 times. I might do that with a for-loop and store the values in a new vector:
212 ```{r}
213 sample.means <- 0
214
215 for(i in 1:100){
216   sample.means[i] <- mean(sample(pop.unif, 2, replace=T))
217 }
218 ```
219
220 You could also do the same thing by nesting the sampling step inside of a call to `sapply()` that runs over an arbitrary index (`rep(1, 100)` here).
221
222 ```{r}
223 sample.means.2 <- sapply( rep(1, 100), function (x) { 
224   mean(sample(pop.unif, 2))}
225   )
226 ```
227
228 Then I can plot them: 
229 ```{r}
230 hist(sample.means)
231
232 hist(sample.means.2)
233 ```
234
235 Note that you certainly didn't need to do this twice, but having done so provides a nice illustration of sampling variability!
236
237 **(d) Draw more samples and describe them too**  
238
239 Same `tapply()` command as (c) just changing the sample size
240
241 ```{r}
242 hist(sapply(rep(1, 100), function (x) { mean(sample(pop.unif, 10))}))
243
244 hist(sapply(rep(1, 100), function (x) { mean(sample(pop.unif, 100))}))
245 ```
246
247 Those axis labels are just terrible, but let's try to focus on the substance of the plots. Some notable things you might observe include that the sampling distribution of the means approaches normality as each sample gets larger in size (and this is true whether the population we draw from is uniform, log-normal, or really just about any other smooth distribution). In this simulation, the number of samples is constant (100), so the changes in the distribution are *solely* due to changes in sample size. This is an illustration of some aspects of the *central limit theorem* applied to sampling distributions. It is also an illustration of the *t-distribution* (the basis for "t-tests" that you'll read about soon).  
248
249 # Reading Questions  
250 ## RQ1. CIs vs. P-values?  
251
252 We'll discuss this one as a group and I'm eager to hear others' thoughts. Personally, I am inclined to agree with Reinhart as I find the focus on p-values somewhat thought-stopping and would prefer to see researchers and publications embrace reporting standards that yield more intuitive, flexible, uncertain, and meaningful interpretations of findings in terms of the original measurements/units. Confidence intervals usually accomplish these goals more effectively than p-values and, in that respect, I like confidence intervals quite a lot. 
253
254 ## RQ2 Emotional contagion revisited  
255 **(a) Hypotheses**    
256
257 In my words (or rather formulas since I think that's less ambiguous), the key pairs of null/alternative hypotheses look something like the following:
258
259 Let $\Delta$ be the parameter estimate for the difference in mean percentage of positive ($\mu_{pos}$) and negative ($\mu_{neg}$) words between the experimental and control conditions for the treatments of reduced negative content ($R_{neg}$ and reduced positive content ($R_{pos}$).
260
261 For the reduced negative content conditions (the left-hand side of Figure 1), the paper tests:
262
263 $$HR_{neg}1_0: \Delta_{\mu_{pos}} = 0$$
264 $$HR_{neg}1_a: \Delta{\mu_{pos}} \neq 0$$
265 And:
266 $$HR_{neg}2_0: \Delta_{\mu_{neg}} = 0$$
267 $$HR_{neg}2_a: \Delta_{\mu_{neg}} \neq 0$$
268 Then, for the reduced positive content conditions (the right-hand side of Figure 1), the paper tests:
269
270 $$HR_{pos}1_0:~~ \Delta_{\mu_{pos}} = 0$$
271 $$HR_{pos}1_a:~~ \Delta{\mu_{pos}} \neq 0$$
272
273 And:
274
275 $$HR_{pos}2_0:~~ \Delta_{\mu_{neg}} = 0$$
276 $$HR_{pos}2_a:~~ \Delta_{\mu_{neg}} \neq 0$$
277 Note that the theories the authors used to motivate the study imply directions for the alternative hypotheses, but nothing in the description of the analysis suggests that they used one-tailed tests. I've written these all in terms of undirected two-tailed tests to correspond to the analysis conducted in the paper. That said, given that the theories correspond to specific directions you might (arguably more accurately) have written the hypotheses in terms of directed inequalities (e.g., "$\gt$" or "$\lt$").
278
279
280 **(b) Describing the effect**s  
281
282 The authors' estimates suggest that reduced negative News Feed content causes an increase in the percentage of positive words and a decrease in the percentage of negative words in subsequent News Feed posts by study participants (supporting $HR_{neg}1_a$ and $HR_{neg}2_a$ respectively).
283
284 They also find that reduced positive News Feed content causes a decrease in the percentage of negative words and an increase in the percentage of positive words in susbequent News Feed posts (supporting $HR_{pos}1_a$ and $HR_{pos}2_a$)
285
286 **(c) Statistical vs. practical significance**  
287
288 Cohen's $d$ puts estimates of experimental effects in standardized units (much like a Z-score!) in order to help understand their size relative to the underlying distribution of the dependent variable(s). The $d$ values for each of the effects estimated in the paper are 0.02, 0.001, 0.02, and 0.008 respectively (in the order presented in the paper, not in order of the hypotheses above). These are miniscule by the standards of most approaches to Cohen's $d$! However, as the authors' argue, the treatment itself is quite narrow in scope, suggesting that the presence of any treatment effect at all is an indication of the underlying phenomenon (emotional contagion). Personally, I find it difficult to attribute much substantive significance to the results because I'm not even convinced that tiny shifts in the percentage of positive/negative words used in News Feed updates accurately index meaningful emotional shifts (I might call it linguistic contagion instead?). That said, I have a hard time thinking about micro-level psychological processes and I'm probably being overly narrow/skeptical in my response. Despite these concerns and the ethical considerations that attracted so much public attention, I consider this a clever, well-executed study and I think it's quite compelling. I expect many of you will have different opinions of various kinds and I'm eager to hear about them.  

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