updates to cover t-tests and ANOVA

[stats_class_2019.git] / problem_sets / week_03 / ps3-worked_solution.Rmd

---
title: 'Week 3 Problem set: Worked solutions'
subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
author: "Aaron Shaw"
date: "April 18, 2019"
output: html_document
---

## Programming challenges

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
### PC1 & PC2: Download and read the data into R

I like to load libraries at the beginning of my R scripts. Doing so helps me find them later. Since I know we're graphing stuff later in this problem set I'll call the ggplot2 package here.

```{r}
library(ggplot2)
```

Now, I'll go ahead and load the CSV file into R. As with last week, I'll do this directly with a `url()` command. However, I included a couple of lines (commented out) that you might adapt to set the working directory on your machine, ask R to show you what's there, and read the data. I'll walk through this in the screencast and/or in class, but **please note that you'll need to edit this code to get it to work on your own file system!**

```{r}

### Uncomment, edit, and try running this on your machine:
###
### setwd("~/Documents/Teaching/2019/stats/")
### list.files("data/week_03") # just take a look around
### w3.data <- read.csv("data/week_03/group_01.csv")

w3.dtata <- read.csv(url("https://communitydata.cc/~ads/teaching/2019/stats/data/week_03/group_02.csv"))
```

### PC3. Get to know your data! 

First, I like to see the first few rows of the dataset and the dimensions of it:

```{r}
head(w3.data)
nrow(w3.data)
ncol(w3.data)
```

Now I'll use the `lapply()` command to ask for some summary information about each variable (column) in the data frame.

```{r}
lapply(w3.data, class)
lapply(w3.data, summary)
lapply(w3.data, sd, na.rm=TRUE)

```

This code below is not ideal since it makes it hard to figure out which variable is being graphed and outputs a bunch of extra crud at the end (the figures are in order, but are all labeled as x[[l]]).  That said, the code sure is expedient and short! We'll learn more intelligible ways to do this later.


```{r}
names(w3.data)
lapply(w3.data, hist)

### You might also have done this variable by variable with code like the following:
### hist(w3.data$x)

```


### PC4. Roll your own functions to find the mean and median of "x"

First, I have to define the `my.mean()` function. Then I can run it.

```{r}

### Run this first...
my.mean <- function(z){
  z <- z[!is.na(z)] 
  sigma <- sum(z)
  n <- length(z)
  out.value <- sigma/n
  return(out.value)
}

### Now you can call the function we just defined
my.mean(w3.data$x)

```

Now, for the median bit. Medians can be complicated because the calculation depends on whether there's a midpoint of the data or not. Luckily in this case, there should be 95 non-missing values in the x vector, so a midpoint exists (the 48th value) and the following should work:

```{r}

my.median <- function(z){
  z <- z[!is.na(z)]
  midpoint <- (length(z) / 2) + .5
  sorted <- sort(z)
  output <- sorted[midpoint]
  return(output)
}

my.median(w3.data$x)

```
Since R has built-in functions that do both of these things, I can even check my answers:
```{r}
mean(w3.data$x, na.rm=T)
median(w3.data$x, na.rm=T)
```



### PC5. Load the Week 2 data and clean it up again.

```{r}
load(url("https://communitydata.cc/~ads/teaching/2019/stats/data/week_02/group_02.RData"))

### I'll rename it for clarity:
w2.data <- d
rm(d)

### and recode:
w2.data[w2.data < 0] <- NA
w2.data <- log1p(w2.data)
```

### PC6. Compare Week 2 data against the x variable from Week 3

First, some summary statistics. These look quite similar...
```{r}
summary(w2.data)
summary(w3.data$x)
```

But they don't quite match up...

```{r}
table(w2.data == w3.data$x)

head(w2.data)
head(w3.data$x)
```

Inspecting the first few values returned by `head()` gave you a clue. Rounded to six decimal places, the vectors match!

I can create a table comparing the sorted rounded values to check this.
```{r}

table(round(w2.data,6) == round(w3.data$x,6))
```

Can you explain what each piece of that last line of code is doing?


### PC7. Visualize the data using ggplot2

First create a basic plotting object. Then add the call to `geom_point()` to show just the x and y: 
```{r}
p <- ggplot(w3.data, mapping=aes(x=x, y=y))

p + geom_point() 
```

Now uncomment this line to try to add the color, size, and shape to the point layer:

```{r}
# p + geom_point(aes(color=j, size=l, shape=k))
```

Hmm. That doesn't work, but the error message is actually pretty helpful. If you search the text of that message online you might discover that we should try turning the aesthetic mappings of the call to `geom_point` into factors.

```{r}
p + geom_point(aes(color=as.factor(j), size=as.factor(l), shape=as.factor(k)))
```

Now that's more like it!


### PC8 Recoding (again)

```{r}
w3.data$j <- as.logical(w3.data$j)
w3.data$l <- as.logical(w3.data$l)

### Create a new column just so I can double check this bit:
w3.data$k.factor <- factor(w3.data$k,
                                 levels=c(0,1,2,3),
                                 labels=c("none", "some", "lots", "all"))

### Spotcheck to make sure it looks good
head(w3.data, 10)

### Now cleanup my extra column:
w3.data$k <- w3.data$k.factor
w3.data$k.factor <- NULL

head(w3.data)
```

### PC9. Summarize again

```{r}
lapply(w3.data, summary) 

### Run this line again to assign the new dataframe to p
p <- ggplot(data=w3.data, mapping=aes(x=x, y=y))

p + geom_point(aes(color=j, size=l, shape=k))
```


Check out how much more readable the legends look now!


## Statistical questions

### SQ1 — 3.4  

(a) Let $X$ denote the finishing times of *Men, Ages 30 - 34* and $Y$ denote the finishing times of *Women, Ages 25 - 29*. Then,

$$X ∼ N (μ = 4313, σ = 583)$$ 
$$Y ∼ N (μ = 5261, σ = 807)$$


(b) Z-scores are a standardization measure: to calculate for a given value you subtract the mean of the corresponding distribution from the value and divide by the standard deviation. The formula notation is given in the *OpenIntro* textbook.

    Let's let R calculate it for us:
```{r}
## Mary
(5513 - 5261) / 807

## Leo:
(4948 - 4313) / 583
```
Since the Z score tells you how many standard deviation units above/below the mean each value is, we can see that Mary finished 0.31 standard deviations above the mean in her category while Leo finished 1.09 standard deviations above the mean in his. 

(c) Mary finished in a much faster time with respect to her category.

(d) Using the Z-score table (Table 3.8) on p. 132 of the book, Leo finished *faster* than approximately $1-0.86 = .14$ or $14\%$ of his category. This corresponds the probability $P(Z \gt 1.09)$ for a normal distribution.

(e) Mary finished *faster* than approximately $1-0.62 = .38$ or $38\%$ of her category. This corresponds to the probability $P(Z \gt 0.31)$ for a normal distribution.

(f) The answer for part b would not change as standardized values (Z-scores) can be computed for any distribution. However, the interpretation and percentile calculations (parts c-e) *would* be different because they all presume a normal distribution. 


### SQ2 — 3.6  

(a) The fastest $5\%$ are in the $5^{th}$ percentile of the distribution. Using Table 3.8 again, the Z score corresponding to the $5^{th}$ percentile of the normal distribution is approximately -1.64. Then, 
$$Z = −1.65 = \frac{x − 4313}{583} → x = −1.65 × 583 + 4313 = 3351~seconds$$
    Divide that by 60 and it looks like the fastest $5\%$ of males in this age group finished in less than 56 minutes.

(b) The slowest $10\%$ are in the $90^{th}$ percentile of the distribution. The Z score corresponding to the $90^{th}$ percentile of the normal distribution is approximately 1.28. Then,
$$Z = 1.27 = \frac{y-5261}{807} → y = 1.28 \times 807 + 5261 = 6294 ~seconds$$
    Divide that by 60 and it looks like the slowest $10\%$ of females in this age group finished in about 1 hour 45 minutes (or longer).  
  
### SQ3 — 3.18  

(a) I'll do the math for this one in R. First, I can enter the data, then can use the information provided to calculate some ranges for the middle 68, 95, and 99 percent of the data (one, two, and three standard deviations from the mean in a normal distribution). Finally, I'll see what percentage of the empirical distribution falls within the ranges of those quantiles.

```{r}

d <- c(54, 55, 56, 56, 57, 58, 58, 59, 60, 60, 60, 61, 61, 62, 62, 63, 63, 63, 64, 65, 65, 67, 67, 69, 73)
m <- 61.52
sdev <- 4.58

## Here are the ranges for +/- 1, 2, and 3 SDs from the mean
r68 <- c(m-sdev, m+sdev)
r95 <- c(m-(2*sdev), m+(2*sdev))
r99 <- c(m-(3*sdev), m+(3*sdev))

## 
l68 <- length(d[d>r68[1] & d<r68[2]])
l95 <- length(d[d>r95[1] & d<r95[2]])
l99 <- length(d[d>r99[1] & d<r99[2]])

l68/length(d)
l95/length(d)
l99/length(d)
```
Those proportions look consistent with the 68-95-99 rule to me. 

(b) The data look pretty consistent with the normal distribution from the visualizations as well. The historgram is pretty symmetric around the theoretical mean of the density plot (red line). The quantile-quantile plot also looks pretty consistent with a normal distribution, altough most of the points that deviate seem to do so in the same direction and there's a slight trend away from the diagonal towards the upper end of the distribution  (right side of the plot).

### SQ4 — 3.32

This question focuses on applying the knowledge from section 3.4 of the textbook on binomial distributions. Review this section and especially equation 3.40 if this gets confusing


(a) $P(at~least~1~arachnophobe)=1-P(none)$  
    $1-P(none)=1-{10 \choose 0}0.07^{0}(1-0.07)^{10-0}$  
```{r}
1-(choose(10,0)*1*(.93^10))
```

(b) $P(2~arachnophobes)={10 \choose 2}0.07^2(1-0.07)^{(10-2)}$  
```{r}
choose(10,2)*0.07^2*0.93^8
```
(c) $P(\leq1~arachnophobes)=P(none)+P(one)$  
    ${10 \choose 0}0.07^00.93^{10}+{10 \choose 1}0.07^1 0.93^9$  

```{r}    
(choose(10,0)*1*(.93^10))+(choose(10,1)*0.07*(0.93^9))
```

(d) Is random assignment to tents likely to ensure $\leq1~arachnophobe$ per tent?
   
 Random assignment and the independence assumption means that the answer to part c is the complement of the outcome we're looking to avoid: $P(\gt1~arachnophobes) = 1-P(\leq1~arachnophobe)$. So, $P(\gt1~arachnophobes) = 1-0.84 = 0.16 = 16\%$. Those are the probabilities, but the interpretation really depends on how confident the camp counselor feels about a $16\%$ chance of having multiple arachnophobic campers in one of the tents.