Introduction to Biostatistics

1 Learning objectives

1.1 Learning objectives

  • Describe the logic of ANOVA
  • Conduct and interpret the overall ANOVA test
  • Conduct post hoc tests

2 Extending the \(t\)-test

2.1 Review: \(t\)-test

  • Compare exactly 2 groups
  • Focus on difference between group means
  • \(t\) distribution
  • Statistical test: difference / standard error

2.2 ANOVA compares 3+ groups

  • Medication: Control, low dose, high dose
  • Treatment: Usual care, new treatment 1, new treatment 2
  • Modality: In-person, online, hybrid

2.3 2 group approaches…

Code 
ggplot(data.frame(x = c(-3, 7)), aes(x)) +
  stat_function(fun = dnorm, args = list(mean = 0, sd = 1), 
                geom = "polygon", 
                alpha = 0.3, 
                fill = "blue") +
  stat_function(fun = dnorm, args = list(mean = 2, sd = 1), 
                geom = "polygon", 
                alpha = 0.3,
                fill = "red") +
  ylim(0,.5) + 
  scale_x_continuous(breaks = -3:7) +
  labs(x = "X", y = "") +
  annotate("text", x = 0, y = 0.2, label = "Group 1", color = "blue") +
  annotate("text", x = 2, y = 0.2, label = "Group 2", color = "red") +
  annotate("segment", x = 0,, xend = 2, y = .4, yend = .4,
           arrow = arrow(ends = "both", angle = 90, length = unit(.2,"cm"))) +
  annotate("text", x = 1, y = .47, label = "Difference")

2.4 … don’t extend well to 3+ groups

Code 
ggplot(data.frame(x = c(-3, 7)), aes(x)) +
  stat_function(fun = dnorm, args = list(mean = 0, sd = 1), 
                geom = "polygon", 
                alpha = 0.3, 
                fill = "blue") +
  stat_function(fun = dnorm, args = list(mean = 2, sd = 1), 
                geom = "polygon", 
                alpha = 0.3,
                fill = "red") +
  stat_function(fun = dnorm, args = list(mean = 4, sd = 1), 
                geom = "polygon", 
                alpha = 0.3,
                fill = "forestgreen") +
  ylim(0,.5) + 
  scale_x_continuous(breaks = -3:7) +
  labs(x = "X", y = "") +
  annotate("text", x = 0, y = 0.2, label = "Group 1", color = "blue") +
  annotate("text", x = 2, y = 0.2, label = "Group 2", color = "red") +
  annotate("text", x = 4, y = 0.2, label = "Group 3", color = "forestgreen") +
  annotate("segment", x = 0, xend = 1.95, y = .4, yend = .4,
           arrow = arrow(ends = "both", angle = 90, length = unit(.2,"cm"))) +
  annotate("text", x = 2, y = .47, label = "Differences?") +
  annotate("segment", x = 2.05, xend = 4, y = .4, yend = .4,
           arrow = arrow(ends = "both", angle = 90, length = unit(.2,"cm"))) +
  annotate("segment", x = 0, xend = 4, y = .43, yend = .43,
           arrow = arrow(ends = "both", angle = 90, length = unit(.2,"cm")))

3 Between-subjects ANOVA

3.1 Logic

  • Extension of independent samples \(t\)-test
  • Sources of variability instead of differences
    • Variability across groups
    • Variability across individuals within a group (random error)
  • Is the variability across groups (i.e., mean differences) large in the context of the within-group variability?
    • Signal to noise ratio

3.2 Between > within

Code 
ggplot(data.frame(x = c(-3, 7)), aes(x)) +
  stat_function(fun = dnorm, args = list(mean = 0, sd = 1), 
                geom = "polygon", 
                alpha = 0.3, 
                fill = "blue") +
  stat_function(fun = dnorm, args = list(mean = 4, sd = 1), 
                geom = "polygon", 
                alpha = 0.3,
                fill = "red") +
  ylim(0,.5) + 
  scale_x_continuous(breaks = -3:7) +
  labs(x = "X", y = "")

3.3 Within > between

Code 
ggplot(data.frame(x = c(-3, 7)), aes(x)) +
  stat_function(fun = dnorm, args = list(mean = 1.5, sd = 2), 
                geom = "polygon", 
                alpha = 0.3, 
                fill = "blue") +
  stat_function(fun = dnorm, args = list(mean = 2.5, sd = 2), 
                geom = "polygon", 
                alpha = 0.3,
                fill = "red") +
  ylim(0,.5) + 
  scale_x_continuous(breaks = -3:7) +
  labs(x = "X", y = "")

3.4 Total variance

Code 
library(Stat2Data)
data(Amyloid)
#head(Amyloid)
ggplot(data = Amyloid, aes(x = Abeta)) +
  geom_dotplot(method = "histodot") +
  theme(legend.position="none")

3.5 Partitioned variance

Code 
ggplot(data = Amyloid, aes(x = Abeta, fill= Group)) +
  geom_dotplot(method = "histodot", stackgroups = TRUE) +
  theme(legend.position="none")

3.6 Partitioned variance

Code 
means <- Amyloid %>% group_by(Group) %>% summarize(meanAbeta = mean(Abeta))
ggplot(data = Amyloid, aes(x = Abeta, fill= Group)) +
  geom_dotplot(method = "histodot", stackgroups = TRUE) +
  theme(legend.position="none") +
  geom_vline(aes(xintercept = means$meanAbeta[1]), 
             color = "#F8766D",
             linewidth = 1) +
  geom_vline(aes(xintercept = means$meanAbeta[2]), 
             color = "#00BA38",
             linewidth = 1) +
  geom_vline(aes(xintercept = means$meanAbeta[3]), 
             color = "#619CFF",
             linewidth = 1)

3.7 Assumptions

  • One categorical (grouping) variable with 3+ categories
  • One continuous (i.e., interval or ratio) variable
  • Data are randomly sampled from the population
  • Data are independent (within and across groups)
  • Data are approximately normally distributed OR sample size is large enough for normally distributed sampling distribution
    • Residuals from model are normally distributed
  • Population variance is unknown and same in all groups

3.8 Hypotheses

  • \(H_0\): All means are equal
    • \(\mu_1 = \mu_2 = \mu_3\)
  • \(H_1\): At least one mean is different from the others
    • (NOT \(\mu_1 \ne \mu_2 \ne \mu_3\))

3.9 ANOVA summary table

Source of variation Sums of squares (SS) df Mean squares (MS) F
Between \(SS_{between}\) \(k - 1\) \(MS_{between}\) \(F\)
Error \(SS_{error}\) \(N - k\) \(MS_{error}\)
Total \(SS_{total}\) \(N - 1\)
  • \(k\) = number of groups
  • \(N\) = total number of participants
  • \(n\) = number of participants in a group

3.10 Test statistic 1: Sums of squares

  • \(SS_{between}\) = variation of \(\color{red}{group\ means}\) around \(\color{blue}{grand\ mean}\)
    • \(SS_{between} = \Sigma_{1}^{k} (\color{red}{\bar{X}_j} - \color{blue}{\bar{X}})^2\)
  • \(SS_{error}\) = variation of \(individual\ values\) around \(\color{red}{group\ means}\)
    • \(SS_{error} = \Sigma_{1}^{k}\Sigma_{1}^{n} (X - \color{red}{\bar{X}_j})^2\)
  • \(SS_{total}\) = variation of \(individual\ values\) around \(\color{blue}{grand\ mean}\)
    • \(SS_{total} = \Sigma_{1}^N (X - \color{blue}{\bar{X}})^2\)
  • \(SS_{total} = SS_{between} + SS_{error}\)

3.11 Test statistic 2: Mean squares

Mean squares = Sum of squares divided by degrees of freedom

  • \(MS_{between}\) = between group variance
    • \(MS_{between} = \frac{SS_{between}}{k - 1}\)
  • \(MS_{error}\) = within group variance
    • \(MS_{error} = \frac{SS_{error}}{N-k}\)

3.12 Test statistic 3: \(F\)-statistic

Analysis of “variance” = compare variances

  • \(F\)-statistic = between group variance / within group variance
    • \(F = \frac{MS_{between}}{MS_{error}}\)
  • Compare to critical \(F\) from a distribution of \(k - 1\) and \(N - k\) degrees of freedom
    • If the observed \(F\) statistic is larger, reject \(H_0\)
    • Otherwise, retain \(H_0\)

3.13 \(F\) distribution

  • Two different degrees of freedom values
    • First is about the number of groups: \(k - 1\)
    • Second is about the number of subjects (and groups): \(N - k\)

Note

  • Compare to \(t\) distribution
    • Just a single df value
    • Always only comparing 2 groups
    • \([t(df)]^2 = F(1, df)\)

3.14 \(F\) distribution

Code 
ggplot(data.frame(x = c(0, 5)), aes(x)) +
  stat_function(fun = df, 
                args = list(df1 = 2, df2 = 47), 
                geom = "line",
                linewidth = 1,
                color = "red") +
  stat_function(fun = df, 
                args = list(df1 = 4, df2 = 45), 
                geom = "line",
                linewidth = 1,
                 color = "blue") +
stat_function(fun = df,
              args = list(df1 = 9, df2 = 40),
              geom = "line",
              linewidth = 1,
              color = "forestgreen") +
  scale_x_continuous(breaks = 0:5) +
  labs(x = "F", y = "") +
  annotate("text", 
           x = 3, y = 0.8, 
           label = "F(2,47)", 
           color = "red",
           size = 6) +
  annotate("text", 
           x = 3, y = 0.6, 
           label = "F(4,45)", 
           color = "blue",
           size = 6) +
  annotate("text",
           x = 3, y = 0.4,
           label = "F(9,40)",
           color = "forestgreen",
           size = 6)
  • 50 subjects in 3 vs 5 vs 10 groups

3.15 Example data

  • \(N\) = 57 observations on the following 2 variables.
    • Group
      • mAD=Alzheimer’s disease
      • MCI=mild impairment
      • NCI=no impairment
    • Abeta
      • Amount of amyloid beta from the posterior cingulate cortex (pmol/g tissue)

3.16 Example data

Code 
Amyloid %>% group_by(Group) %>%
  summarize(mean = mean(Abeta),
            variance = var(Abeta)) %>%
  kable()
Group mean variance
mAD 761.2941 182068.0
MCI 341.0476 165168.4
NCI 336.2632 189755.8

3.17 Check assumptions: Equal variance

  • Graphically
Code 
ggplot(data = Amyloid,
       aes(x = Group, y = Abeta)) +
  geom_boxplot(alpha = 0.5) +
  geom_jitter(width = 0.2)

  • leveneTest() in car package
Code 
library(car)
leveneTest(Abeta ~ Group, data = Amyloid)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  0.3505 0.7059
      54
  • \(H_0\): No differences in variances across groups
    • NS result = equal variances

3.18 Example: ANOVA

Code 
m1 <- lm(Abeta ~ Group,
         contrasts=list(Group='contr.sum'),
         data = Amyloid)
anova(m1)
Analysis of Variance Table

Response: Abeta
          Df  Sum Sq Mean Sq F value   Pr(>F)   
Group      2 2129969 1064985  5.9706 0.004544 **
Residuals 54 9632060  178371                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • We reject the null hypothesis that all means are equal, F(2, 54) = 5.9706, p = .0045. At least one group mean is different from the others.

3.19 Type I vs III sums of squares

  • Type III sums of squares are “corrected sums of squares”
    • How ANOVA is taught, default in every software package except R
  • Type I sums of squares are “uncorrected sums of squares”
    • Default in R
  • Different interpretations
    • But not different numerically if all groups are the same size
    • Type III recommended by almost everyone when groups are different sizes

3.20 Assumptions: Normality

  • Graphically
Code 
m1resid <- data.frame(residuals(m1))
ggplot(data = m1resid, 
       aes(sample = residuals.m1.)) +
  geom_qq() +
  geom_qq_line()

Code 
m1resid <- data.frame(residuals(m1))
ggplot(data = m1resid, 
       aes(x = residuals.m1.)) +
  geom_histogram()

  • shapiro.test() in stats
Code 
shapiro.test(resid(m1))

    Shapiro-Wilk normality test

data:  resid(m1)
W = 0.88246, p-value = 0.00004857
  • \(H_0\): Observations come from normally distributed population
    • NS result = normally distributed residuals

3.21 Post hoc tests

  • Significant overall (omnibus) test for ANOVA means
    • At least one group mean is different from the others
    • But which one(s)?
  • Post hoc tests
    • Figure out which means are different from one another
    • But not all possible tests, only specific ones

3.22 Post hoc tests

  • Conceptually, a series of \(t\)-tests
    • Use standard error estimate based on whole model instead of standard error estimate from just the two groups
    • Context of whole model
      • Not just some random set of 2 group comparisons

3.23 Contrasts for post hoc tests

  • Many built-in contrasts for specific patterns
    • Four common: dummy, deviation, polynomial, Helmert
    • Also: User-specified contrasts
  • Contrasts / weights: UCLA Statistical Consulting
  • ANOVA table is the same regardless of which contrasts you use
    • summary() of model will give you mean differences

3.24 Example again: ANOVA

Code 
m1 <- lm(Abeta ~ Group,
         contrasts=list(Group='contr.sum'),
         data = Amyloid)
anova(m1)
Analysis of Variance Table

Response: Abeta
          Df  Sum Sq Mean Sq F value   Pr(>F)   
Group      2 2129969 1064985  5.9706 0.004544 **
Residuals 54 9632060  178371                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.25 Contrasts: Dummy / treatment

  • Compare each group to a specific group (e.g., control)
  mAD v MCI mAD v NCI
mAD 0 0
MCI 1 0
NCI 0 1
  • Default reference: First

3.26 Contrasts: Dummy / treatment

Code 
m1_trt <- lm(Abeta ~ Group,
         contrasts=list(Group='contr.treatment'),
         data = Amyloid)
#anova(m1_trt)
summary(m1_trt)

Call:
lm(formula = Abeta ~ Group, data = Amyloid, contrasts = list(Group = "contr.treatment"))

Residuals:
   Min     1Q Median     3Q    Max 
-623.3 -322.3 -108.3  227.9 1269.0 

Coefficients:
            Estimate Std. Error t value       Pr(>|t|)    
(Intercept)    761.3      102.4   7.432 0.000000000819 ***
GroupMCI      -420.2      137.8  -3.050        0.00354 ** 
GroupNCI      -425.0      141.0  -3.014        0.00392 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 422.3 on 54 degrees of freedom
Multiple R-squared:  0.1811,    Adjusted R-squared:  0.1508 
F-statistic: 5.971 on 2 and 54 DF,  p-value: 0.004544

3.27 Contrasts: Deviation / sum

  • Compare each group (except one) to unweighted grand mean
    • \((\bar{X}_1 + \bar{X}_2 + \bar{X}_3)/3\)
  mAD v GM MCI v GM
mAD 1 0
MCI 0 1
NCI -1 -1
  • Default left out: Last

3.28 Contrasts: Deviation / sum

Code 
m1_sum <- lm(Abeta ~ Group,
         contrasts=list(Group='contr.sum'),
         data = Amyloid)
#anova(m1_sum)
summary(m1_sum)

Call:
lm(formula = Abeta ~ Group, data = Amyloid, contrasts = list(Group = "contr.sum"))

Residuals:
   Min     1Q Median     3Q    Max 
-623.3 -322.3 -108.3  227.9 1269.0 

Coefficients:
            Estimate Std. Error t value        Pr(>|t|)    
(Intercept)   479.53      56.15   8.540 0.0000000000134 ***
Group1        281.76      81.55   3.455         0.00108 ** 
Group2       -138.49      77.36  -1.790         0.07902 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 422.3 on 54 degrees of freedom
Multiple R-squared:  0.1811,    Adjusted R-squared:  0.1508 
F-statistic: 5.971 on 2 and 54 DF,  p-value: 0.004544

3.29 Contrasts: Polynomial / poly

  • Linear, quadatric, etc. effects (orthogonal)
  Linear Quadratic
mAD -1 1
MCI 0 -2
NCI 1 1
  • Only useful with ordered levels (none, low, high)

3.30 Contrasts: Polynomial / poly

Code 
m1_poly <- lm(Abeta ~ Group,
              contrasts=list(Group='contr.poly'),
              data = Amyloid)
#anova(m1_poly)
summary(m1_poly)

Call:
lm(formula = Abeta ~ Group, data = Amyloid, contrasts = list(Group = "contr.poly"))

Residuals:
   Min     1Q Median     3Q    Max 
-623.3 -322.3 -108.3  227.9 1269.0 

Coefficients:
            Estimate Std. Error t value        Pr(>|t|)    
(Intercept)   479.53      56.15   8.540 0.0000000000134 ***
Group.L      -300.54      99.70  -3.014         0.00392 ** 
Group.Q       169.61      94.74   1.790         0.07902 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 422.3 on 54 degrees of freedom
Multiple R-squared:  0.1811,    Adjusted R-squared:  0.1508 
F-statistic: 5.971 on 2 and 54 DF,  p-value: 0.004544

3.31 Contrasts: Helmert / helmert

  • Compare 1st to mean of 2 to k, 2nd to mean of 3 to k, etc.
  mAD v MCI mAD v NCI
mAD 2/3 0
MCI -1/3 1/2
NCI -1/3 -1/2

3.32 Contrasts: Helmert / helmert

Code 
m1_helm <- lm(Abeta ~ Group,
              contrasts=list(Group='contr.helmert'),
              data = Amyloid)
#anova(m1_helm)
summary(m1_helm)

Call:
lm(formula = Abeta ~ Group, data = Amyloid, contrasts = list(Group = "contr.helmert"))

Residuals:
   Min     1Q Median     3Q    Max 
-623.3 -322.3 -108.3  227.9 1269.0 

Coefficients:
            Estimate Std. Error t value        Pr(>|t|)    
(Intercept)   479.53      56.15   8.540 0.0000000000134 ***
Group1       -210.12      68.90  -3.050         0.00354 ** 
Group2        -71.64      39.63  -1.808         0.07623 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 422.3 on 54 degrees of freedom
Multiple R-squared:  0.1811,    Adjusted R-squared:  0.1508 
F-statistic: 5.971 on 2 and 54 DF,  p-value: 0.004544

3.33 Mean comparisons

Code 
ggplot(data = Amyloid,
       aes(x = Group, 
           y = Abeta, 
           color = Group)) +
  geom_jitter(width = 0.2) +
  geom_point(data = means, 
             aes(x = Group, 
                 y = meanAbeta, 
                 color = Group,
                 size = 10)) +
  geom_hline(yintercept = sum(means$meanAbeta)/3, 
             linetype = "dashed") +
  theme(legend.position="none")

  • Dummy: mAD v MCI, mAD v NCI
  • Sum: mAD v GM, not MCI v GM
  • Poly: Linear, not quadratic
  • Helmert: 1 v 2&3, not 2 v 3

3.34 Should you correct for that?

  • Arguments for YES
    • You’re doing multiple tests
    • Inflated type I error rate
  • Arguments for NO
    • Significant overall ANOVA means “there is a difference”
    • Specific, theoretically-based comparisons
  • If you do correct, use methods we’ve already talked about

3.35 Welch ANOVA

  • Accounts for unequal variances across groups: Heterogeneity
  • oneway.test() in stats package
oneway.test(Abeta ~ Group, data = Amyloid, var.equal = FALSE)

    One-way analysis of means (not assuming equal variances)

data:  Abeta and Group
F = 5.7936, num df = 2.000, denom df = 35.141, p-value = 0.006689

3.36 Kruskal Wallis

  • Non-parametric version of ANOVA
  • kruskal.test() in stats package
kruskal.test(Abeta ~ Group, data = Amyloid)

    Kruskal-Wallis rank sum test

data:  Abeta by Group
Kruskal-Wallis chi-squared = 11.967, df = 2, p-value = 0.00252

4 In-class activities

4.1 In-class activities

  • Conduct an ANOVA
  • Do some post hoc tests
  • Compare some means

4.2 Next week

  • Last week of lecture! Woo!
  • Factorial ANOVA
    • Multiple grouping variables
      • Condition and sex
      • Medication 1 dose and medication 2 dose
  • Within-subjects ANOVA
    • Extension of paired \(t\)-test