Lecture 8

class: center, middle, inverse, title-slide

# Lecture 8
### Jung-Jin Lee
### Mar 5, 2019

---

## One-way ANOVA, review

Recall `iris` data set:

```r
head(iris)
```

```
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
```

Can we say that `Sepal.Width` is associated with `Species`?

---

## One-way ANOVA, review

```r
ggplot(iris, aes(Species, Sepal.Width)) +
  geom_boxplot()
```

---

## One-way ANOVA, review

```r
oneway.test(iris$Sepal.Width ~ iris$Species, var.equal = T)
```

```
## 
## 	One-way analysis of means
## 
## data:  iris$Sepal.Width and iris$Species
## F = 49.16, num df = 2, denom df = 147, p-value < 2.2e-16
```

```r
summary(aov(iris$Sepal.Width ~ iris$Species))
```

```
##               Df Sum Sq Mean Sq F value Pr(>F)    
## iris$Species   2  11.35   5.672   49.16 <2e-16 ***
## Residuals    147  16.96   0.115                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

- Small `$p$`-value indicates that we should _reject_ the null hypothesis, i.e., population means are not the same for all species.
- In other words, at least one mean is different from another.
- Which pair(s) show difference?

---

## Pairwise _t_-tests

- Compare all possible pairs

- Done by applying two-sample `$t$`-tests for each pair

```r
setosa.SW <- iris %>% 
  filter(Species == "setosa") %>%
  pull(Sepal.Width)

versicolor.SW <- iris %>% 
  filter(Species == "versicolor") %>%
  pull(Sepal.Width)

virginica.SW <- iris %>% 
  filter(Species == "virginica") %>%
  pull(Sepal.Width)
```

---

## Pairwise _t_-tests

```r
t.test(setosa.SW, versicolor.SW)
```

```
## 
## 	Welch Two Sample t-test
## 
## data:  setosa.SW and versicolor.SW
## t = 9.455, df = 94.698, p-value = 2.484e-15
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.5198348 0.7961652
## sample estimates:
## mean of x mean of y 
##     3.428     2.770
```

---

## Pairwise _t_-tests

```r
t.test(setosa.SW, virginica.SW)
```

```
## 
## 	Welch Two Sample t-test
## 
## data:  setosa.SW and virginica.SW
## t = 6.4503, df = 95.547, p-value = 4.571e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.3142808 0.5937192
## sample estimates:
## mean of x mean of y 
##     3.428     2.974
```

---

## Pairwise _t_-tests

```r
t.test(versicolor.SW, virginica.SW)
```

```
## 
## 	Welch Two Sample t-test
## 
## data:  versicolor.SW and virginica.SW
## t = -3.2058, df = 97.927, p-value = 0.001819
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.33028364 -0.07771636
## sample estimates:
## mean of x mean of y 
##     2.770     2.974
```

---

## Pairwise _t_-tests

Short cut for all pairwise comparisons using `pairwise.t.test()`:

```r
*pairwise.t.test(iris$Sepal.Width, iris$Species,
*               pool.sd = F, p.adjust.method = "none")
```

```
## 
## 	Pairwise comparisons using t tests with non-pooled SD 
## 
## data:  iris$Sepal.Width and iris$Species 
## 
##            setosa  versicolor
## versicolor 2.5e-15 -         
## virginica  4.6e-09 0.0018    
## 
## P value adjustment method: none
```

---

## Multiple comparisons

- Suppose we want to test if a coin is fair (i.e. probability of getting `$H$` is `$\frac{1}{2}$`)

- As a test:

- Toss a coin five times, and say the outcome is `$\{H,H,H,H,H\}$`
  - If the coin is fair, the probability of getting `$\{H,H,H,H,H\}$` is given by 
  `$$\frac{1}{2} \times \frac{1}{2} \times \frac{1}{2} \times \frac{1}{2} \times \frac{1}{2} = \frac{1}{32} = 0.03125,$$` or 3.125%
  - It is a rare event and so we have an evidence to claim that the coin IS biased in favor of `$H$`.

- Still it is possible this happened by chance

- About 1 out of 20 fair coins (notice that `$\frac{1}{20}$` = 5%) can show this behavior.

---

## Multiple comparisons

Simulation of 20 unbiased coins, five tosses each

```r
set.seed(2019)  
sam <- sample(c("H", "T"), size = 100, replace = T)  
df <- matrix(sam, ncol = 5) %>%  
  as.data.frame() %>%  
  setNames(paste0("Toss", 1:5))

dim(df)
```

```
## [1] 20  5
```

```r
head(df)
```

```
##   Toss1 Toss2 Toss3 Toss4 Toss5
## 1     T     H     T     H     H
## 2     T     T     H     H     T
## 3     H     T     T     T     T
## 4     T     H     T     H     H
## 5     H     H     T     T     T
## 6     H     H     H     T     H
```

---

## Multiple comparisons

<div id="htmlwidget-7aca68d55ea43992310d" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-7aca68d55ea43992310d">{"x":{"filter":"none","fillContainer":false,"data":[["1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20"],["T","T","H","T","H","H","T","H","H","T","T","T","H","H","T","T","H","T","H","H"],["H","T","T","H","H","H","T","T","T","H","T","H","T","H","H","H","H","H","H","T"],["T","H","T","T","T","H","T","T","H","H","H","H","H","H","T","H","T","T","H","H"],["H","H","T","H","T","T","H","H","T","T","H","H","H","H","T","T","H","T","H","T"],["H","T","T","H","T","H","T","T","T","H","T","T","T","H","T","T","H","T","H","T"]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th> <\/th>\n      <th>Toss1<\/th>\n      <th>Toss2<\/th>\n      <th>Toss3<\/th>\n      <th>Toss4<\/th>\n      <th>Toss5<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"pageLength":7,"order":[],"autoWidth":false,"orderClasses":false,"columnDefs":[{"orderable":false,"targets":0}],"lengthMenu":[7,10,25,50,100]}},"evals":[],"jsHooks":[]}</script>

---

## Multiple comparisons

- Coin 14 (and 19) should not be viewed as a biased coin!

- https://xkcd.com/882/

- False discovery (or Type I error): rejecting the null hypothesis when it should not be

- `$p$`-values from multiple comparisons should be penalized to reduce false discoveries

---

## One-way ANOVA, post hoc analysis

```r
anova_res <- aov(iris$Sepal.Width ~ iris$Species)
*TukeyHSD(anova_res)
```

```
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = iris$Sepal.Width ~ iris$Species)
## 
## $`iris$Species`
##                        diff         lwr        upr     p adj
## versicolor-setosa    -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa     -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor  0.204  0.04314472  0.3648553 0.0087802
```

All pairs show a significant difference in means.

---

## One-way ANOVA, post hoc analysis

**Example**: recall `crabs.tsv` which has data involving features of female horseshoe crabs.

```r
crabs <- read.table("crabs.tsv", header = T, sep = "\t") 
crabs <- crabs %>%
  mutate(color = factor(color)) %>%
  mutate(spine = factor(spine)) %>%
  mutate(id = paste0("obs", 1:nrow(crabs)))
head(crabs)
```

```
##   color spine width satell weight   id
## 1     3     3  28.3      8   3050 obs1
## 2     4     3  22.5      0   1550 obs2
## 3     2     1  26.0      9   2300 obs3
## 4     4     3  24.8      0   2100 obs4
## 5     4     3  26.0      4   2600 obs5
## 6     3     3  23.8      0   2100 obs6
```

---

## One-way ANOVA, post hoc analysis

```r
ggplot(crabs, aes(color, weight)) +
  geom_boxplot()
```

---

## One-way ANOVA, post hoc analysis

```r
anova_res <- aov(crabs$weight ~ crabs$color)
summary(anova_res)
```

```
##              Df   Sum Sq Mean Sq F value  Pr(>F)   
## crabs$color   3  3766792 1255597   3.966 0.00917 **
## Residuals   169 53502001  316580                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## One-way ANOVA, post hoc analysis

```r
hsd <- TukeyHSD(anova_res)
print(hsd)
```

```
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = crabs$weight ~ crabs$color)
## 
## $`crabs$color`
##           diff       lwr       upr     p adj
## 3-2  -91.35614 -538.6378 355.92548 0.9516931
## 4-2 -329.91667 -805.3818 145.54843 0.2767204
## 5-2 -455.30303 -979.2401  68.63399 0.1129741
## 4-3 -238.56053 -504.7929  27.67184 0.0964136
## 5-3 -363.94689 -709.3779 -18.51583 0.0346278
## 5-4 -125.38636 -506.6065 255.83381 0.8286679
```

```r
names(hsd)
```

```
## [1] "crabs$color"
```

---

## One-way ANOVA, post hoc analysis

```r
df <- hsd$`crabs$color` %>%
  as.data.frame() %>%
  rownames_to_column("Pair") %>%
  mutate(Sig = ifelse(`p adj` < 0.05, "Sig.", "Not Sig."))
```

```r
ggplot(df, aes(Pair, color = Sig)) +
  geom_hline(yintercept = 0) +
* geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.2) +
  labs(color = "") 
```

.pull-left[
<img src="Drexel_2019_Lecture_8_files/figure-html/unnamed-chunk-20-1.png" width="90%" style="display: block; margin: auto;" />
]

.pull-right[
- `lwr` is the lower limit of 95% confidence interval of mean difference.
- `upr` is the upper limit of 95% confidence interval of mean difference.
- The pair color 5 and color 3 shows a significant difference in weights.
]

---

## False discoveries

```r
pairwise.t.test(crabs$weight, crabs$color, 
                pool.sd = F, p.adjust = "none")
```

```
## 
## 	Pairwise comparisons using t tests with non-pooled SD 
## 
## data:  crabs$weight and crabs$color 
## 
##   2      3      4     
## 3 0.4845 -      -     
## 4 0.0259 0.0233 -     
## 5 0.0049 0.0029 0.3284
## 
## P value adjustment method: none
```

Without penalties on `$p$`-values, we could have false discoveries that the color pairs (2, 4), (3, 4), (2, 5) are also showing differences in weight.

---

## Chi-squared distribution

---

## Chi-squared test

**Example**: 40 mice were randomly selected from 100 male mice and were castrated one day after birth. They were compared with the remaining 60 sham-operated mice. The mice were maintained for 140 days and blood samples were collected biweekly starting at 42 days old. The frequencies of mice with and without diabetes are recorded in the `$2\times 2$` contingency table below:

`$$\begin{array}{c|ccc} & \text{Castrated Mice} & \text{Control Mice} & \text{Total} \\ \hline   
\text{With Diabetes} & 22 & 15 & 37 \\  
\text{Without Diabetes} & 18 & 45 & 63 \\
\text{Total} & 40 & 60 & 100
\end{array}$$`
It was shown that neonatal castration more than doubled the incidence of diabetes (22/40=55%) when compared with controls (15/60=25%). Is this difference significant? That is, can we say that two population proportions of diabetes are different?

This example is different in that the response variable is not continuous, but categorical (With/Without Diabetes).

---

## Chi-squared test, continued

- Observed contingency table

- Overall proportion of diabetes: 37%

- If two treatments (Castrated/Control) do not affect the incidence of diabetes, then the proportion diabetes in both treatments should be the same 37%.

- In other words, the expected number of diabetic mice should be about 40 `$\times$` 0.37 = 14.8 in castrated group and 60 `$\times$` 0.37 = 22.2 in control group.

- Similarly, the expected number of non-diabetic mice should be about 40 `$\times$` 0.63 = 25.2 in castrated group and 60 `$\times$` 0.63 = 37.8 in control group.

---

## Chi-squared test, continued

- Observed contingency table

- Expected contingency table under the assumption that treatments do not affect the incidence of diabetes

`$$\begin{array}{c|ccc} & \text{Castrated Mice} & \text{Control Mice} & \text{Total} \\ \hline   
\text{With Diabetes} & 14.8 & 22.2 & 37 \\  
\text{Without Diabetes} & 25.2 & 37.8 & 63 \\
\text{Total} & 40 & 60 & 100
\end{array}$$`

---

## Chi-squared test, continued

If treatments do not affect the incidence of diabetes, then observed values should not be too different from expected. The following is true in general:

`$$\begin{array}{c|ccc} & \text{X} & \text{Y} & \text{Total} \\ \hline   
\text{Success} & a & b & a+b \\  
\text{Failure} & c & d & c+d \\
\text{Total} & n=a+c & m=b+d & n+m=a+b+c+d
\end{array}$$`

Let `$p_0=(a+b)/(a+b+c+d)=(a+b)/(n+m)$`. Under the assumption that two treatments `$X$` and `$Y$` have the same success rate, 
`$$\chi^2 = \frac{(a-np_0)^2}{np_0}+\frac{(b-mp_0)^2}{mp_0}+\frac{(c-n(1-p_0))^2}{n(1-p_0)}+\frac{(d-m(1-p_0))^2}{m(1-p_0)}$$`
approximately follows a `$\chi^2$`-distribution with `$1$` degree of freedom.

---

## Chi-squared test, continued

- Observed contingency table

- Expected contingency table

```r
X2 <- (22-14.8)^2/14.8 + (15-22.2)^2/22.2 +
  (18-25.2)^2/25.2 + (45-37.8)^2/37.8; print(X2)
```

```
## [1] 9.266409
```

---

## Chi-squared test, continued

```r
pval <- 1-pchisq(X2, df = 1); print(pval)
```

```
## [1] 0.002333948
```

With small `$p$`-value, we conclude that the treatments differently affect on diabetes.

---

## Chi-squared test, Yates' continuity correction

Instead of `$$\chi^2 = \frac{(a-np_0)^2}{np_0}+\frac{(b-mp_0)^2}{mp_0}+\frac{(c-n(1-p_0))^2}{n(1-p_0)}+\frac{(d-m(1-p_0))^2}{m(1-p_0)},$$`
Yates suggested the use of 
`$$\chi^2 = \frac{(|a-np_0|-0.5)^2}{np_0}+\frac{(|b-mp_0|-0.5)^2}{mp_0} \\+\frac{(|c-n(1-p_0)|-0.5)^2}{n(1-p_0)}+\frac{(|d-m(1-p_0)|-0.5)^2}{m(1-p_0)}.$$`

```r
X2_Yates <- 
  (abs(22-14.8)-0.5)^2/14.8 + (abs(15-22.2)-0.5)^2/22.2 +
  (abs(18-25.2)-0.5)^2/25.2 + (abs(45-37.8)-0.5)^2/37.8
print(X2_Yates)
```

```
## [1] 8.024096
```

---

## Chi-squared test, Yates' continuity correction

```r
pval <- 1-pchisq(X2_Yates, df = 1); print(pval)
```

```
## [1] 0.004615907
```

Again with small `$p$`-value, we conclude that the treatments differently affect on diabetes.

---

## Chi-squared test, continued

R time!

```r
contingency <- matrix(c(22, 18, 15, 45), byrow = F, nrow = 2)
print(contingency)
```

```
##      [,1] [,2]
## [1,]   22   15
## [2,]   18   45
```

---

## Chi-squared test, continued

```r
chisq.test(contingency, correct = F) # No Yates' correction
```

```
## 
## 	Pearson's Chi-squared test
## 
## data:  contingency
## X-squared = 9.2664, df = 1, p-value = 0.002334
```

```r
chisq.test(contingency, correct = T) # With Yates' correction
```

```
## 
## 	Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency
## X-squared = 8.0241, df = 1, p-value = 0.004616
```

Default: With Yates' correction

---

## Chi-squared test, continued

**Example**: 993 cows were randomly assigned to one of three insemination methods. Each cow was inseminated once and the proportion of cows that became pregnant in each group is given in the following table:
`$$\begin{array}{c|cccc} & \text{Method I} & \text{Method II} & \text{Method III} & \text{Total} \\ \hline   
\text{Pregnant} & 275 & 192 & 261 & 728 \\  
\text{Not Pregnant} & 78 & 64 & 123 & 265 \\
\text{Total} & 353 & 256 & 384 & 993
\end{array}$$`
Is there any evidence for believing that the insemination methods show different proportions of pregnant animals?

Note that we have more than two treatments.

---

## Chi-squared test, continued

Suppose that the following `$r\times c$` contingency table is given:
`$$\begin{array}{c|ccccccc} & G_1 & G_2 & \cdots & G_j & \cdots & G_c & \text{Total} \\ \hline   
\text{Outcome 1} & n_{11} & n_{12} & \cdots & n_{1j} & \cdots & n_{1c} & n_{1+} \\  
\text{Outcome 2} & n_{21} & n_{22} & \cdots & n_{2j} & \cdots & n_{2c} & n_{2+} \\  
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\text{Outcome i} & n_{i1} & n_{i2} & \cdots & n_{ij} & \cdots & n_{ic} & n_{i+} \\ 
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\text{Outcome r} & n_{r1} & n_{r2} & \cdots & n_{rj} & \cdots & n_{rc} & n_{r+} \\
\text{Total} & n_{+1} & n_{+2} & \cdots & n_{+j} & \cdots & n_{+c} & n
\end{array}$$`
Under the assumption that population outcome distributions are the same for all groups, 
`$\sum_{i=1}^r\sum_{j=1}^c \frac{(n_{ij}-n_{+j}\frac{n_{i+}}{n})^2}{n_{+j}\frac{n_{i+}}{n}} \quad \left(\sum \frac{(Observed-Expected)^2}{Expected}\right)$`
approximately follows a `$\chi^2$`-distribution with `$(r-1)(c-1)$` degrees of freedom.

---

## Chi-squared test, continued

```r
obs <- matrix(c(275, 78, 192, 64, 261, 123), 
               nrow = 2, byrow = F)
n_1_plus <- 728; n_2_plus <- 265;
n_plus_1 <- 353; n_plus_2 <- 256; n_plus_3 <- 384
n <- 993
expect <- matrix(1/n*c(n_plus_1*n_1_plus, n_plus_1*n_2_plus, 
                       n_plus_2*n_1_plus, n_plus_2*n_2_plus,
                       n_plus_3*n_1_plus, n_plus_3*n_2_plus), 
                 nrow = 2, ncol = 3)
print(obs)
```

```
##      [,1] [,2] [,3]
## [1,]  275  192  261
## [2,]   78   64  123
```

```r
print(expect)
```

```
##           [,1]      [,2]     [,3]
## [1,] 258.79557 187.68177 281.5227
## [2,]  94.20443  68.31823 102.4773
```

---

## Chi-squared test, continued

```r
X2 <- sum((obs-expect)^2/expect) 
# Yates' correction is only for 2 x 2
pval <- 1-pchisq(X2, df = 2)
print(c(X2, pval))
```

```
## [1] 9.780370691 0.007520029
```

---

## Chi-squared test, continued

```r
chisq.test(obs)
```

```
## 
## 	Pearson's Chi-squared test
## 
## data:  obs
## X-squared = 9.7804, df = 2, p-value = 0.00752
```

Conclusion: the insemination methods show different proportions of pregnant animals.