---
title: "ADA1: Class 23, ANOVA, Pairwise comparisons"
author: Your Name
date: last-modified
description: |
[Advanced Data Analysis 1](https://StatAcumen.com/teach/ada1),
Stat 427/527, Fall 2023, Prof. Erik Erhardt, UNM
format:
html:
theme: litera
highlight-style: atom-one
page-layout: full # article, full # https://quarto.org/docs/output-formats/page-layout.html
toc: true
toc-location: body # body, left, right
number-sections: false
self-contained: false # !!! this can cause a render error
code-overflow: scroll # scroll, wrap
code-block-bg: true
code-block-border-left: "#30B0E0"
code-copy: false # true, false, hover a copy buttom in top-right of code block
fig-width: 6
fig-height: 4
---
# Rubric
```{r load-packages, message=FALSE}
library(erikmisc)
library(tidyverse)
ggplot2::theme_set(ggplot2::theme_bw()) # set theme_bw for all plots
```
# Example: Child Health and Development Study ([CHDS](https://dash.nichd.nih.gov/study/8))
We consider data from the birth records of 680 live-born white
male infants. The infants were born to mothers who reported for
pre-natal care to three clinics of the Kaiser hospitals in
northern California. As an initial analysis, we will examine
whether maternal smoking has an effect on the birth weights of
these children. To answer this question, we define 3 groups based
on mother's smoking history:
1. mother does not currently smoke or never smoked,
2. mother smoked less than one pack of cigarettes a day during pregnancy, and
3. mother smoked at least one pack of cigarettes a day during pregnancy.
Let $\mu_i$ = pop mean birth weight (lb) for children in group $i$, $(i=1,2,3)$.
We wish to test $H_0: \mu_1 = \mu_2 = \mu_3$ against $H_A: \textrm{not } H_0$.
We read in the data, create a `smoke` factor variable,
and plot the data by smoking group.
```{R}
#### Example from the Child Health and Development Study (CHDS)
# description at http://statacumen.com/teach/ADA1/ADA1_notes_05-CHDS_desc.txt
dat_chds <-
read_csv("ADA1_notes_05-CHDS.csv")
dat_chds <-
dat_chds %>%
mutate(
# create a factor variable based on number of cigarettes smoked
smoke = case_when(
# no cigs
m_smok == 0 ~ "0 cigs"
# less than 1 pack (20 cigs = 1 pack)
, (m_smok > 0) & (dat_chds$m_smok < 20) ~ "1-19 cigs"
# at least 1 pack (20 cigs = 1 pack)
, m_smok >= 20 ~ "20+ cigs"
)
, smoke = factor(smoke)
) %>%
select(
c_bwt
, smoke
)
summary(dat_chds)
```
Plot the data in a way that compares the means.
Error bars are 95% confidence intervals of the mean.
```{R}
#| fig-width: 3
#| fig-height: 5
# Plot the data using ggplot
library(ggplot2)
p <- ggplot(dat_chds, aes(x = smoke, y = c_bwt))
# plot a reference line for the global mean (assuming no groups)
p <- p + geom_hline(yintercept = mean(dat_chds$c_bwt),
colour = "black", linetype = "dashed", size = 0.3, alpha = 0.5)
# boxplot, size=.75 to stand out behind CI
p <- p + geom_violin(width = 0.5, alpha = 0.25)
p <- p + geom_boxplot(width = 0.25, alpha = 0.25)
# points for observed data
p <- p + geom_point(position = position_jitter(w = 0.05, h = 0), alpha = 0.2)
# diamond at mean for each group
p <- p + stat_summary(fun = mean, geom = "point", shape = 18, size = 4,
colour = "red", alpha = 0.8)
# confidence limits based on normal distribution
p <- p + stat_summary(fun.data = "mean_cl_normal", geom = "errorbar",
width = .2, colour = "red", alpha = 0.8)
p <- p + labs(x = "Maternal smoking (per day)"
, y = "Child birthweight (lb)"
, title = "Child birthweight vs\nmaternal smoking"
)
print(p)
```
__Hypothesis test__
## (0 p) 1. Set up the __null and alternative hypotheses__ in words and notation.
* In words: ``The population mean birthweight is different between smoking groups.''
* In notation: $H_0: \mu_1=\mu_2=\mu_3$ versus $H_A: \textrm{not } H_0$ (at least one pair of means differ).
## (0 p) 2. Let the significance level of the test be $\alpha=0.05$.
## (0 p) 3. Compute the __test statistic__.
```{R}
fit_c <-
aov(
c_bwt ~ smoke
, data = dat_chds
)
summary(fit_c)
```
The $F$-statistic for the ANOVA is $F = `r signif(unlist(summary(fit_c))["F value1"], 3)`$.
## (0 p) 4. Compute the __$p$-value__ from the test statistic.
The p-value for testing the null hypothesis is
$p = `r signif(unlist(summary(fit_c))["Pr(>F)1"], 3)`$.
## (2 p) 5. State the __conclusion__ in terms of the problem.
[ANSWER HERE]
## 6. __Check assumptions__ of the test.
### (1 p) Plot the residuals and assess whether they appear normal.
a. Residuals are normal
b. Populations have equal variances.
* Check whether residuals are normal.
```{R}
#| fig-width: 5
#| fig-height: 4
# Plot the data using ggplot
df_res <- data.frame(res = fit_c$residuals)
library(ggplot2)
p <- ggplot(df_res, aes(x = res))
p <- p + geom_histogram(aes(y = ..density..), binwidth = 0.2)
p <- p + geom_density(colour = "blue")
p <- p + geom_rug()
p <- p + stat_function(fun = dnorm, colour = "red", args = list(mean = mean(df_res$res), sd = sd(df_res$res)))
p <- p + labs(title = "ANOVA Residuals"
, caption = "Blue = Kernal density curve\nRed = Normal distribution")
print(p)
```
Describe the plot of residuals as it relates to model assumptions.
[ANSWER HERE]
### (1 p) Describe the plot of residuals as it relates to model assumptions.
- Plot the residuals versus the normal quantiles.
If the residuals are normal, then the will fall on the center line and
very few will be outside the error bands.
```{R}
#| fig-width: 5
#| fig-height: 5
# QQ plot
par(mfrow=c(1,1))
#library(car)
car::qqPlot(
fit_c$residuals
, las = 1
, id = list(n = 4, cex = 1)
, lwd = 1
, main="QQ Plot"
)
```
[ANSWER HERE]
- A formal test of normality on the residuals tests the hypothesis
$H_0:$ The distribution is Normal vs
$H_1:$ The distribution is not Normal.
We can test the distribution of the residuals.
Three tests for normality are reported below.
I tend to like the Anderson-Darling test.
Different tests have different properties, and
tests that are sensitive to differences from normality in the tails of the
distribution are typically more important for us (since deviations in the tails
are more influential than deviations in the center).
```{R}
shapiro.test(fit_c$residuals)
#library(nortest)
nortest::ad.test(fit_c$residuals)
nortest::cvm.test(fit_c$residuals)
```
### (1 p) Interpret the conclusion of the Anderson-Darling test.
[ANSWER HERE]
### Check whether populations have equal variances.
### (1 p) Interpret the standard deviations.
You may also discuss the plots of the data.
- Look at the numerical summaries below.
```{R}
# calculate summaries
dat_chds_summary <-
dat_chds %>%
group_by(smoke) %>%
summarize(
m = mean(c_bwt)
, s = sd(c_bwt)
, n = n()
, .groups = "drop_last"
) %>%
ungroup()
dat_chds_summary
```
[ANSWER HERE]
### (1 p) Interpret the result of the appropriate equal variance test.
- Formal tests for equal variances.
We can test whether the variances are equal between our three groups.
This is similar to the ANOVA hypothesis, but instead of testing means we're tesing variances.
$H_0: \sigma^2_1=\sigma^2_2=\sigma^2_3$
versus $H_A: \textrm{not } H_0$ (at least one pair of variances differ).
```{R}
## Test equal variance
# assumes populations are normal
bartlett.test(c_bwt ~ smoke, data = dat_chds)
# does not assume normality, requires car package
#library(car)
car::leveneTest(c_bwt ~ smoke, data = dat_chds)
# nonparametric test
fligner.test(c_bwt ~ smoke, data = dat_chds)
```
If normality was reasonable then use Bartlett, otherwise use Levene.
[ANSWER HERE]
## 7. Pairwise comparisons
If the ANOVA null hypothesis was rejected, then perform follow-up Post Hoc
pairwise comparison tests to determine which pairs of means are different.
There are several multiple comparison methods described in the notes.
Let's use Tukey's Honest Significant Difference (HSD) here to test which pairs of
populations differ.
__EMM plot interpretation__
This __EMM plot (Estimated Marginal Means, aka Least-Squares Means)__
is only available when conditioning on one variable.
The blue bars are confidence intervals for the EMMs;
don't ever use confidence intervals for
EMMs to perform comparisons -- they can be very misleading.
The red arrows are for the comparisons among means;
the degree to which the "comparison arrows" overlap reflects as much as
possible the significance of the comparison of the two estimates.
If an arrow from one mean overlaps an arrow from
another group, the difference is not significant, based on the adjust setting
(which defaults to "tukey").
```{R}
#| fig-width: 3
#| fig-height: 5
## CHDS
# Tukey 95% Individual p-values
#TukeyHSD(fit_c)
## Contrasts
adjust_method <- c("none", "tukey", "bonferroni")[2]
library(emmeans)
emm_cont <-
emmeans::emmeans(
fit_c
, specs = "smoke"
)
# means and CIs
emm_cont %>% print()
# pairwise differences
emm_cont %>% pairs(adjust = adjust_method) %>% print()
# plot of means, CIs, and comparison arrows
plot(
emm_cont
, comparisons = TRUE
, adjust = adjust_method
, horizontal = FALSE
, ylab = "smoke"
)
```
### (2 p) Interpret the comparisons (which pairs differ).
[ANSWER HERE]
* `0 cigs - (1-19 cigs)`: means are (not) significant different
* `0 cigs - (20+ cigs)`: means are (not) significant different
* `(1-19 cigs) - (20+ cigs)`: means are (not) significant different
### (1 p) Summarize results by ordering the means and grouping pairs that do not differ (see notes for examples).
[ANSWER HERE]
```
Replace this example with your results.
Example: Groups A and C differ, but B is not different from either.
(These groups are ordered by mean, so A has the lowest mean and C has the highest.)
Group A Group B Group C
-----------------
-----------------
```