]> code.communitydata.science - stats_class_2019.git/blobdiff - problem_sets/week_06/ps6-worked-solution.Rmd
updates to cover t-tests and ANOVA
[stats_class_2019.git] / problem_sets / week_06 / ps6-worked-solution.Rmd
index f5791b0a0b6b80e85abe67c24192cbeb2bbf8275..986a9a0147f27890d9b6c2925e9898f9204666e2 100644 (file)
@@ -65,7 +65,7 @@ h_plot
 
 # In this case, faceted histograms is probably better
 
-h_facet = df %>% ggplot(aes(x=weeks_alive, # What to summarize
+h_facet = df %>% ggplot(aes(x=weeks_alive # What to summarize
                       )) + geom_histogram(bins = 5) + facet_grid(~dose)
 
 h_facet
@@ -90,7 +90,7 @@ box_plot
 library(ggridges) 
 
 ridge_plot = df %>% ggplot(aes(x=weeks_alive, y = dose)) + 
-  geom_density_ridges(jittered_points = T)
+  geom_density_ridges(jittered_points = T, fill = 'orange') + theme_minimal()
 ridge_plot
 
 ```
@@ -103,7 +103,17 @@ The global mean is
 mean(df$weeks_alive)
 ```
 
-PC3. T-test between None and Any, and between None and High.
+
+PC3. Anova
+
+```{r}
+summary(aov(weeks_alive ~ dose, data = df))
+
+```
+
+This provides evidence that the group means are different.
+
+PC4. T-test between None and Any, and between None and High.
 
 ```{r}
 
@@ -129,16 +139,10 @@ t.test(weeks_alive ~ dose, data = tmp)
 
 ```
 
-The t-test supports the idea that receiving a dose of RD40 reduces lifespan
-
-PC4. Anova
-
-```{r}
-summary(aov(weeks_alive ~ dose, data = df))
+These t-tests both support the idea that receiving a dose of RD40 reduces lifespan. However, we should not completely trust these p-values, since we are doing multiple comparisons. One option is to do a Bonferroni correction, where we only cnsider things significant if $\alpha < .05/m$ where $m$ is the number of tests. In this case, we would fail to reject the null for the second test because we would set $\alpha = .025$.
 
-```
+The Bonferroni correction is more conservative than it needs to be, ane there are other approaches; for example,  the `TukeyHSD` function takes in an anova result and does post-hoc comparisons with corrections for all of the groups.
 
-This provides evidence that the group means are different.
 
 ## Statistical Questions
 
@@ -202,6 +206,9 @@ chisq.test(x = c(83,121,193,103), p = c(.18,.22,.37,.23))
 # Use the formula for chi-squared
 chisq = (83-90)^2/90 + (121-110)^2/110 + (193-185)^2/185 + (103-115)^2/115
 
+chisq
+
+# We could then look up the chi-square distribution for 3 degrees of freedom
 
 ```
 
@@ -218,3 +225,66 @@ x <- matrix(c(29,54,44,77,62,131,36,67), nrow = 4, # this makes a matrix with 4
             byrow=T) # And this says that we've entered it row by row
 chisq.test(x)
 ```
+
+## Empirical Questions
+
+
+EQ0.
+
+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 unkown gender customers are being compared. This is a two-way test.
+
+b) The null hypothesis is that the board purchased is independent of the gender of the customer. The alternative hypothesis is that if we know the gender of the customer that will tell us something about the type of board they purchased.
+
+c) A $\chi^2$ test found statistically significant evidence that board purchase behavior differs by gender. This difference is convincing, but it does directly not tell us what we really want to know, which is the difference between men and women. It could be possible that it is simply identifying 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 reality of this difference.
+
+d) Statistical tests help to give (or take away) confidence in a conclusion. People are not natively good at estimating how likely something is due to chance and tests help us to make these judgments. 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 numeric data (e.g., the amount of money spent on electronics for men and women) then we could choose a t-test to compare those distributions.
+
+
+EQ1
+
+a) These are counts for two categorical variables, so the procedure used was a $\chi^2$ test. The null hypothesis is that whether or not a blog is governed by one person is independent of whether it is on the left or the right.
+
+b) It would be surprising to see these results by chance and it make sense to believe that this difference is real. However, the main reason to be skeptical is the way that the data are grouped. 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 then we have good reason to be skeptical.
+
+c) 
+
+```{r}
+
+# First we create the dataframe
+df = data.frame(Governance=c('Individual','Multiple', 'Individual','Multiple'),
+                Ideology = c('Left','Left','Right','Right'),
+                Count = c(13,51,27,38))
+
+# We can make sure it's the same by testing the Chi-squared                
+chisq.test(matrix(df$Count, nrow=2))
+
+percentage_data = df %>% 
+  group_by(Ideology) %>% 
+  summarize(individual_ratio = sum(Count[Governance=='Individual']) / sum(Count),
+            group_count = sum(Count))
+
+shaw_benkler_plot = percentage_data %>%
+  ggplot(aes(x=Ideology, y=individual_ratio * 100)) + 
+    geom_bar(stat='identity', aes(fill=c('red','blue')), show.legend=F) + 
+    ylab('Percentage of Blogs') + theme_minimal()
+
+shaw_benkler_plot
+
+# 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 - always share your data!)
+
+# I decided to use confidence intervals. The standard error is another reasonable choice
+
+# Remember that for a binomial distribution (we can consider individual/non-individual as binomial), confidence intervals are mu +- x * sqrt((p*(1-p))/n)
+
+ci_95 = 1.96 * sqrt(percentage_data$individual_ratio * (1 - percentage_data$individual_ratio)/percentage_data$group_count)
+
+shaw_benkler_plot + geom_errorbar(aes(ymin=(individual_ratio-ci_95)*100, ymax=(individual_ratio + ci_95)*100),
+                                  alpha = .3, 
+                                  size=1.1, 
+                                  width=.4)
+
+```
+
+The error bars do overlap in this case, but that actually doesn't tell us whether they are significantly different.
+
+d) We don't need to be very worried about the base rate fallacy because the sizes of both groups are about the same.

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