problem_sets/week_07/ps7-worked-solutions.Rmd

   1 ---
   2 title: "Week 7 problem set: Worked solutions"
   3 author: "Jeremy Foote & Aaron Shaw"
   4 date: "May 16, 2019"
   5 output: html_document
   6 ---
   7
   8 ```{r setup, include=FALSE}
   9 knitr::opts_chunk$set(echo = TRUE)
  10 ```
  11 ## Programming Challenges
  12
  13 ### PC1-3
  14
  15 We'll import the .dta file first using an appropriate command from the `readstata13` package. Also, after looking through the dataverse files, it turns out there is a version of the data which is a TSV file, and can be imported with `read_delim()`
  16
  17 ```{r}
  18 library(readstata13)
  19
  20 df <- read.dta13('~/Documents/Teaching/2019/stats/data/week_07/Halloween2012-2014-2015_PLOS.dta')
  21
  22 ## Same result
  23 ## df.t <- read.delim('~/Documents/Teaching/2019/stats/data/week_07/Halloween2012-2014-2015_PLOS.tab')
  24
  25 head(df)
  26 summary(df)
  27 ```
  28 There are a few strange 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.
  29
  30 ### PC4
  31
  32 The `table` function is a great way to create contingency tables. We'll recode the variables as logical TRUE/FALSE values first.
  33
  34 ```{r}
  35
  36 # Change both measures into T/F
  37 df$obama = as.logical(df$obama)
  38 df$fruit = as.logical(df$fruit)
  39
  40 # create the table with nice labels
  41 obama.tbl <- table(fruit=df$fruit, flotus=df$obama)
  42 obama.tbl
  43
  44 ```
  45
  46 ### PC5
  47
  48 The simplest way to determine if the groups are independent is a $\chi^2$ test. Since it's a 2x2 comparison, we could also test for a difference in proportions using the `prop.test()` function.
  49
  50 ```{r}
  51
  52 chisq.test(obama.tbl)
  53
  54 prop.test(obama.tbl)
  55
  56 ```
  57 Notice that both functions report identical $\chi^2$ test results and p-values.
  58
  59 There are many ways you could answer the question about why these results are different from the regression results presented in the paper. We'll discuss some of them in class this week and some more next week as part of our discussion of multiple regression.
  60
  61 ### PC6
  62
  63 First we want to get the proportion and standard error for fruit in each group. These can be calculated individually or using a function (guess which one we'll document here). Also, note that I'm just going to use the `complete.cases()` function to eliminate the missing items for the sake of simplicity.
  64
  65 ```{r}
  66 df <- df[complete.cases(df),]
  67
  68 prop.se = function(values){# Takes in a vector of T/F values
  69   N = length(values)
  70   prop = mean(as.numeric(values))
  71   se = sqrt(prop * (1-prop)/N)
  72   return(c(prop, se))
  73 }
  74
  75 prop.se(df$fruit[df$obama])
  76 prop.se(df$fruit[!df$obama])
  77
  78 ```
  79 In order to graph that it will help to convert the results into a data frame:
  80 ```{r}
  81 library(ggplot2)
  82
  83 prop.and.se <- data.frame(rbind(
  84   prop.se(df$fruit[df$obama]),
  85   prop.se(df$fruit[!df$obama])
  86 ))
  87 names(prop.and.se) <- c("proportion", "se")
  88 prop.and.se$obama <- c(TRUE, FALSE)
  89
  90 ggplot(prop.and.se, aes(x=obama,y=proportion)) +
  91   geom_point(aes(color=obama), size=5) +  # Add the points for the proportions
  92   geom_errorbar(aes(ymin=proportion - 1.96 * se, # Add error bars
  93                     ymax=proportion + 1.96 * se,
  94                     width=0, # Remove the whiskers
  95                     color=obama),size=1.1) + # Make them a little bigger and color them
  96   coord_flip() + # Flip the chart
  97   theme_light() + # Change the theme (theme_minimal is also nice)
  98   scale_color_manual(values = c('gray','black'), guide=F) + # Change the colors
  99   ylim(0,.5) + # Change the y axis to go from 0 to .5
 100   ylab('Proportion choosing fruit') + # Add labels
 101   xlab('Picture shown was Michelle Obama')
 102
 103 ```
 104
 105 Another way to do that involves a slightly different version of the function we created above and then using the `group_by` and `summarize` functions in the `dplyr` library.
 106 ```{r}
 107
 108 prop.se = function(values){# Takes in a vector of T/F values
 109   N = length(values)
 110   prop = mean(values)
 111   se = sqrt(prop * (1-prop)/N)
 112   return(se)
 113 }
 114
 115 library(dplyr)
 116
 117 prop.and.se = df %>% filter(!is.na(fruit)) %>%
 118   group_by(obama) %>%
 119   summarize(
 120     proportion=mean(fruit),
 121     se = prop.se(fruit)
 122   )
 123
 124 ## Same exact plotting code
 125 ggplot(prop.and.se, aes(x=obama,y=proportion)) +
 126   geom_point(aes(color=obama), size=5) +
 127   geom_errorbar(aes(ymin=proportion - 1.96 * se,
 128                     ymax=proportion + 1.96 * se,
 129                     width=0,
 130                     color=obama),size=1.1) +
 131   coord_flip() +
 132   theme_light() +
 133   scale_color_manual(values = c('gray','black'), guide=F) +
 134   ylim(0,.5) +
 135   ylab('Proportion choosing fruit') +
 136   xlab('Picture shown was Michelle Obama')
 137
 138 ```
 139
 140 ### PC7
 141
 142 Here's one way to export our table, using write.csv
 143 ```{r}
 144 write.csv(obama.tbl, file = 'crosstabs.csv')
 145
 146 ```
 147 We can make sure it worked
 148 ```{r}
 149 read.csv('crosstabs.csv')
 150 ```
 151 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:
 152 ```{r}
 153 as.data.frame(obama.tbl)
 154 ```
 155 and then save that dataframe.
 156 ```{r}
 157 write.csv(as.data.frame(obama.tbl), file = 'crosstabs.csv')
 158 read.csv('crosstabs.csv')
 159 ```
 160 You could also use the `xtable` package to do this. The package has many functions to customize table outputs, but a relatively simple way to generate an html table looks like this:
 161 ```{r}
 162 library(xtable)
 163 print(xtable(obama.tbl), type="html")
 164 ```
 165 There is a lot of documentation and examples online to help you customize as you see fit.
 166
 167 ## Statistical Questions
 168
 169 ### SQ1—7.6
 170 a) Husband and wife ages are positively, linearly correlated, with a few possible outliers.
 171 b) Husband and wife heights also appear to be positively and linearly correlated, but with very weakly.
 172 c) The age plot shows a much more evident linear trend. The data points also align more tightly.
 173 d) No. Correlation is not affected by rescaling. Revisit the formulas for calculating correlation for more depth on this..
 174
 175 ### SQ2—7.30
 176 a) $\widehat{heart~weight} = -0.357 + 4.034 \times body_wt$
 177 b) If a cat had a body weight of 0, we would expect it to have a heart weight of -0.357 (which in itself is obviously non-sensical. The intercept just serves to anchor the regression line and often has no substantive meaning in the real world).
 178 c) For every additional kilogram in weight, heart weight is expected to increase 4.034 grams on average.
 179 d) About 64% of the variance in heart weights is explained by body weight.
 180 e) $\sqrt{0.6466} = .8041$
 181
 182 ### SQ3—7.40
 183 a) $\bar{y} = b_0 + b_1 \bar{x}$
 184 $3.9983 = 4.010 + b_1 \times -0.0883$
 185 ```{r}
 186 b_1 = (3.9983-4.010) / -0.0883
 187 b_1
 188 ```
 189 b) The null hypothesis is that $H_0: \beta_1 = 0$. With a t-score of 4.13 and a p-value of \~0.00 there is strong evidence of a positive relationship between beauty and teaching scores.
 190 c)
 191   1. Nearly normal residuals: The distriubtion of residuals looks a little bit skewed but fairly normal.
 192   2. Homoscedasticity of residuals: The residuals have approximately equal variance across the range of beauty around zero.
 193   3. Independent observations/residuals: The question does not really shed much light on this. Looking at the original paper, it is actually a bit unclear how the instructors were chosen and whether or not they are a random sample of all instructors at UT Austin.
 194   4. Linear relationship between independent/dependent variables: From the scatterplot, the relationship appears to be linear (or at least not obviously non-linear).
 195
 196 ### SQ4—7.42
 197 a) Substituting into the regression equation:
 198 $\bar{y} = 3.91 + .78 * 28$
 199 ```{r}
 200
 201 y_bar = 3.91 + .78 * 28
 202
 203 y_bar
 204
 205 ```
 206 b) If our null hypothesis is that $\beta_{gestational\_age} = 0$, we calculate the t-score as
 207 $$T = \frac{0.78 - 0}{.35} = 2.23$$
 208 We can look up a p-value for this in a t-table for $df = 23$ or we can figure it out in R. The `pt(q,df)`  function gives the proportion of the T-distribution which is less than `q` with `df` degrees of freedom.
 209
 210 ```{r}
 211 (1 - pt(q = 2.23, df = 23)) * 2 # To get 2-tailed p-value, we multiply this by 2
 212
 213 pt(q = 2.23, df = 23, lower.tail = FALSE) * 2 # Alternatively, we can get the area of the upper tail and multiply that by 2 (multiplying by 2 works because t-distributions are symmetrical)
 214 ```
 215 The p-value is less than 0.05, so with a traditional $\alpha = 0.05$ we would reject the null hypothesis and conclude that there is substantial evidence of an association between gestational age and head circumference for low birth-weight babies.
 216
 217 ## Empirical Questions
 218
 219 ### EQ1.
 220 Final score had a strong postive correlation with "Modded", "Starting score", and "Karma", and  strong negative relationship with "Anonymous user". Other correlations were fairly weak.
 221
 222 ### EQ2.
 223 The authors have presented a histogram of the dependent variable and while it's not quite normal, it does not seem so skewed that we would be too worried about fitting a linear model. It might be helpful to see more information about the linearity of relationships between the independent variables and the outcome. In addition, the fact that this is data collected over time suggests that there may be temporal dependencies in the data that undermine the assumptions of a linear model. It might be helpful to see a residual plot and QQ-plot to check for homogeneity of variance and normal distribution of errors.
 224
 225 ## EQ3.
 226 The coefficients in the table represent the expected amount that a final score would change for a one-unit change in a given measure (in this sense, multiple regression is the same as regression with a single predictor). For example, if we focus on the first row in the table, for every 1-point increase in starting score, we would expect the final score to increase by 1.08 on average.
 227
 228 The t-statistic is the (coefficient - 0) / standard error, and the P value is the proportion of the t-distribution which is as extreme or more extreme than that t-statistic. For the purposes of the first coefficient, the extremely high t-statistic and extremely low p-value indicate that the probability of observing a relationship this strong under the null hypothesis of no association between starting score and final score is very, very small. This indicates that we can reject the null hypothesis.
 229
 230 The $R^2$ value is the amount of variation in final scores which is explained by these measures. In this case, the $R^2$ value of 0.52 indicates that the variables in the model explain a substantial amount of the variation in the final scores, but that other explanatory factors likely exist that are not captured by the model.