MVS Module 4 Handout: Advanced ANOVA

Example #1: ANCOVA with an experimental design and pretest as covariate

Load in some helpful packages

library(tidyverse)
library(haven)

Load in the dataset

hap_ancova <- read_dta("hap-ancova1.dta")

Explore your data

glimpse(hap_ancova)

Rows: 82
Columns: 4
$ id       <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
$ hap_pre  <dbl> 40, 51, 54, 58, 52, 44, 58, 41, 42, 66, 47, 66, 47, 54,…
$ hap_post <dbl> 49, 56, 51, 52, 68, 53, 47, 49, 47, 58, 36, 59, 49, 54,…
$ treat    <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …

Clean your data

hap_ancova.clean <- hap_ancova %>%
mutate(.,
  treat.fac = as_factor(treat))

Check out descriptive statistics for the treatment factor and continuous variables

The summary package in base R is a quick and easy way to do this!

summary(hap_ancova.clean)

       id           hap_pre        hap_post         treat    
 Min.   : 1.00   Min.   :28.0   Min.   :30.00   Min.   :1.0  
 1st Qu.:21.25   1st Qu.:44.0   1st Qu.:47.00   1st Qu.:1.0  
 Median :41.50   Median :51.5   Median :51.00   Median :1.5  
 Mean   :41.50   Mean   :51.3   Mean   :53.46   Mean   :1.5  
 3rd Qu.:61.75   3rd Qu.:58.0   3rd Qu.:59.00   3rd Qu.:2.0  
 Max.   :82.00   Max.   :81.0   Max.   :77.00   Max.   :2.0  
                  treat.fac 
 Control group         :41  
 Optimism therapy group:41  
                            
                            
                            
                            

Correlate pre and post test results- how similar are they?

Hopefully, somewhat similar! If not, there may not be systematic change over time. The cor.test function in base R works well for simple correlations.

cor.test(hap_ancova.clean$hap_pre, hap_ancova.clean$hap_post)

    Pearson's product-moment correlation

data:  hap_ancova.clean$hap_pre and hap_ancova.clean$hap_post
t = 6.3349, df = 80, p-value = 1.302e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4127435 0.7063892
sample estimates:
      cor 
0.5779818 

Run a basic ANCOVA

Here, we have hap_post as the outcome, treat as the factor, and hap_pre as the covariate. Since we cleaned our data at the beginning, R recognizes that treat is a factor and hap_pre is continuous so it will run as an ANCOVA and not a two-way ANOVA.

ancova1<-aov(hap_post~ treat.fac + hap_pre,data=hap_ancova.clean)
summary(ancova1)

            Df Sum Sq Mean Sq F value   Pr(>F)    
treat.fac    1    226   225.6   3.653   0.0596 .  
hap_pre      1   2700  2700.4  43.730 4.07e-09 ***
Residuals   79   4878    61.8                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

If you had more than 2 groups, you could do post-hoc comparisons:

TukeyHSD(ancova1)

non-factors ignored: hap_pre'which' specified some non-factors which will be dropped

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = hap_post ~ treat.fac + hap_pre, data = hap_ancova.clean)

$treat.fac
                                         diff        lwr      upr
Optimism therapy group-Control group 3.317073 -0.1375479 6.771694
                                         p adj
Optimism therapy group-Control group 0.0596061

Get measures of effect size:

library(effectsize)
eta_squared(ancova1)

 treat.fac    hap_pre 
0.02890181 0.34601162 

omega_squared(ancova1)

Parameter | Omega2 (partial) |       90% CI
-------------------------------------------
treat.fac |             0.03 | [0.00, 0.12]
hap_pre   |             0.34 | [0.21, 0.46]

What if we didn’t have the pretest score?

Now, see the difference without the covariate- just an ANOVA for treat:

anova1<-aov(hap_post~ treat.fac, data=hap_ancova.clean)
summary(anova1)

            Df Sum Sq Mean Sq F value Pr(>F)
treat.fac    1    226  225.56   2.381  0.127
Residuals   80   7579   94.74               

And, of course- how to run this as a multiple regression!

regression1 <- lm(hap_post ~ treat.fac + hap_pre, data = hap_ancova.clean)
summary(regression1)

Call:
lm(formula = hap_post ~ treat.fac + hap_pre, data = hap_ancova.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.7017  -3.7916  -0.0503   4.2741  18.3137 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      20.5479     4.8834   4.208 6.76e-05 ***
treat.facOptimism therapy group   3.9496     1.7382   2.272   0.0258 *  
hap_pre                           0.6031     0.0912   6.613 4.07e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.858 on 79 degrees of freedom
Multiple R-squared:  0.3749,    Adjusted R-squared:  0.3591 
F-statistic: 23.69 on 2 and 79 DF,  p-value: 8.702e-09

Example #2: ANCOVA with experimental design and related covariates (not pretest)

It might be better for external validity to use related measures as covariates, instead of a pretest. ## Load in the dataset:

hap_ancova2 <- read_dta("hap-ancova2.dta")

Explore your data:

glimpse(hap_ancova2)

Rows: 96
Columns: 6
$ id        <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,…
$ hap       <dbl> 34, 44, 58, 58, 50, 51, 48, 64, 36, 54, 42, 55, 41, 40…
$ treat     <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ satfam    <dbl+lbl> 5, 6, 6, 5, 3, 4, 5, 5, 6, 4, 3, 5, 5, 2, 4, 3, 4,…
$ health    <dbl+lbl> 2, 5, 5, 4, 5, 3, 5, 4, 4, 3, 3, 3, 3, 4, 3, 3, 4,…
$ socfriend <dbl+lbl> 6, 4, 3, 2, 5, 5, 5, 7, 5, 4, 4, 4, 5, 6, 4, 5, 5,…

Clean your data:

hap_ancova2.clean <- hap_ancova2 %>%
  mutate(.,
         treat.fac = as_factor(treat))

Explore your data

summary(hap_ancova2.clean)

       id             hap            treat         satfam     
 Min.   : 1.00   Min.   :28.00   Min.   :1.0   Min.   :1.000  
 1st Qu.:24.75   1st Qu.:42.00   1st Qu.:1.0   1st Qu.:3.000  
 Median :48.50   Median :50.50   Median :1.5   Median :4.000  
 Mean   :48.50   Mean   :50.11   Mean   :1.5   Mean   :3.854  
 3rd Qu.:72.25   3rd Qu.:58.00   3rd Qu.:2.0   3rd Qu.:5.000  
 Max.   :96.00   Max.   :81.00   Max.   :2.0   Max.   :7.000  
     health        socfriend                      treat.fac 
 Min.   :2.000   Min.   :2.000   Control group         :48  
 1st Qu.:3.000   1st Qu.:4.000   Optimism therapy group:48  
 Median :4.000   Median :5.000                              
 Mean   :3.708   Mean   :4.594                              
 3rd Qu.:4.000   3rd Qu.:5.000                              
 Max.   :7.000   Max.   :7.000                              

How much variance does our set of covariates explain?

regression2 <- lm(hap ~ satfam + health + socfriend, data = hap_ancova2.clean)
summary(regression2)

Call:
lm(formula = hap ~ satfam + health + socfriend, data = hap_ancova2.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.9955  -6.1508   0.0402   6.5118  19.5972 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  22.1176     5.0588   4.372 3.23e-05 ***
satfam        2.5023     0.6515   3.841 0.000225 ***
health        3.0764     0.8946   3.439 0.000880 ***
socfriend     1.5117     0.7884   1.918 0.058272 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.888 on 92 degrees of freedom
Multiple R-squared:  0.2878,    Adjusted R-squared:  0.2646 
F-statistic: 12.39 on 3 and 92 DF,  p-value: 7.035e-07

Now, plug them into the ANCOVA model:

ancova2 <- aov(hap ~ treat.fac + satfam + health + socfriend, data = hap_ancova2.clean)
summary(ancova2)

            Df Sum Sq Mean Sq F value   Pr(>F)    
treat.fac    1    298   297.5   3.960   0.0496 *  
satfam       1   1504  1503.7  20.014 2.22e-05 ***
health       1   1280  1279.5  17.030 8.14e-05 ***
socfriend    1    286   285.6   3.801   0.0543 .  
Residuals   91   6837    75.1                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Regression model with the variable treat.fac

regression3 <- lm(hap ~ treat.fac + satfam + health + socfriend, data = hap_ancova2.clean)
summary(regression3)

Call:
lm(formula = hap ~ treat.fac + satfam + health + socfriend, data = hap_ancova2.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-18.9098  -6.6984   0.6385   5.9155  17.7141 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      19.1636     5.0862   3.768 0.000292 ***
treat.facOptimism therapy group   4.3025     1.7993   2.391 0.018849 *  
satfam                            2.3761     0.6376   3.727 0.000337 ***
health                            3.4395     0.8856   3.884 0.000195 ***
socfriend                         1.4992     0.7689   1.950 0.054288 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.668 on 91 degrees of freedom
Multiple R-squared:  0.3299,    Adjusted R-squared:  0.3005 
F-statistic:  11.2 on 4 and 91 DF,  p-value: 1.964e-07

Example #3: Repeated Measures ANOVA

Reshaping data from wide to long

Hint for this week’s content review!

sleep <- read_dta("sleep_wide.dta")

glimpse(sleep)

Rows: 100
Columns: 4
$ id     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
$ sleep1 <dbl> 303, 331, 374, 374, 349, 353, 343, 390, 307, 360, 327, 36…
$ sleep2 <dbl> 349, 380, 400, 385, 442, 365, 370, 399, 410, 346, 371, 34…
$ sleep3 <dbl> 382, 350, 345, 356, 366, 372, 354, 352, 375, 403, 336, 40…

summary(sleep)

       id             sleep1          sleep2          sleep3     
 Min.   :  1.00   Min.   :276.0   Min.   :308.0   Min.   :297.0  
 1st Qu.: 25.75   1st Qu.:326.5   1st Qu.:348.0   1st Qu.:348.8  
 Median : 50.50   Median :349.0   Median :368.0   Median :366.5  
 Mean   : 50.50   Mean   :348.2   Mean   :369.3   Mean   :366.7  
 3rd Qu.: 75.25   3rd Qu.:369.2   3rd Qu.:389.0   3rd Qu.:382.2  
 Max.   :100.00   Max.   :434.0   Max.   :466.0   Max.   :437.0  

Using pivot_longer to reshape from wide to long

sleep.long <- sleep %>%
pivot_longer(.,
              cols = starts_with("sleep"),
              names_to = "month",
              values_to = "sleep",
              names_prefix = "sleep",
              )
glimpse(sleep.long)

Rows: 300
Columns: 3
$ id    <dbl> 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7…
$ month <chr> "1", "2", "3", "1", "2", "3", "1", "2", "3", "1", "2", "3"…
$ sleep <dbl> 303, 349, 382, 331, 380, 350, 374, 400, 345, 374, 385, 356…

Reshape back to wide using pivot_wider

sleep.wide <- sleep.long %>%
  pivot_wider(.,
              id_cols = c("id"),
              names_from = month,
              values_from = sleep,
              names_prefix = "sleep",
              )

glimpse(sleep.wide)

Rows: 100
Columns: 4
$ id     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
$ sleep1 <dbl> 303, 331, 374, 374, 349, 353, 343, 390, 307, 360, 327, 36…
$ sleep2 <dbl> 349, 380, 400, 385, 442, 365, 370, 399, 410, 346, 371, 34…
$ sleep3 <dbl> 382, 350, 345, 356, 366, 372, 354, 352, 375, 403, 336, 40…

Density Plot for Month 1

kdensity1 <- ggplot(sleep, aes(x=sleep1)) + 
  geom_density() + 
  stat_function(fun = dnorm, n = 100, args = list(mean = 350, sd = 30), color = "red", linetype = "dashed") +
  labs(title="Distribution of Sleep in Month 1",x="Sleep (minutes per night)", y = "Density", caption = "N = 100 participants. The red line indicates a normal distribution. The black lines represents the observed distribution.")
kdensity1

Density Plot for Month 2

kdensity2 <- ggplot(sleep, aes(x=sleep2)) + 
  geom_density() + 
  stat_function(fun = dnorm, n = 100, args = list(mean = 370, sd = 30), color = "red", linetype = "dashed") +
  labs(title="Distribution of Sleep in Month 2",x="Sleep (minutes per night)", y = "Density", caption = "N = 100 participants. The red line indicates a normal distribution. The black lines represents the observed distribution.")
kdensity2

Density Plot for Month 3

kdensity3 <- ggplot(sleep, aes(x=sleep2)) + 
  geom_density() + 
  stat_function(fun = dnorm, n = 100, args = list(mean = 370, sd = 30), color = "red", linetype = "dashed") +
  labs(title="Distribution of Sleep in Month 3",x="Sleep (minutes per night)", y = "Density", caption = "N = 100 participants. The red line indicates a normal distribution. The black lines represents the observed distribution.")
kdensity3

Using data in long form, we can facet by month to include all three graphs in one table

all_sleep<-ggplot(sleep.long, aes(x=sleep))+
  geom_density()+facet_wrap(month ~ ., ncol = 1)
all_sleep

Summary Statistics of Sleep by Month

library(rstatix)
library(ggpubr)
sleep.long %>%
  group_by(month) %>%
  get_summary_stats(sleep, type = "mean_sd")

Code for basic repeated measures ANOVA

Note:id is participant id in this case.

sleep.aov <- anova_test(data = sleep.long, dv = sleep, wid = id, within = month)
summary(sleep.aov)

                              Length Class      Mode
ANOVA                         7      data.frame list
Mauchly's Test for Sphericity 4      data.frame list
Sphericity Corrections        9      data.frame list

Here is the basic ANOVA table:

get_anova_table(sleep.aov)

ANOVA Table (type III tests)

  Effect DFn DFd      F        p p<.05  ges
1  month   2 198 17.938 6.92e-08     * 0.09

Here are the corrected p-values with the sphericity adjustments:

sleep.aov$`Sphericity Corrections`

Post-hoc comparisons by month, using Tukey’s method:

pairwise_t_test(
    formula = sleep ~ month, paired = TRUE,
    p.adjust.method = "bonferroni",
    data = sleep.long
    )

Check for sphericity and multivariate normality

Start with univariate normality, plots and statistics

library(MVN)
sleep_univariate <- mvn(data = sleep.wide, mvnTest = NULL, univariatePlot = "qqplot")

sleep_univariate2 <- mvn(data = sleep.wide, mvnTest = NULL, univariateTest = "SW", desc = TRUE)
sleep_univariate2$univariateNormality

Tests for MVN

Doornik-Hansen Test

The last column indicates whether dataset follows a multivariate normality or not (i.e, YES or NO) at significance level 0.05.

sleep_mvn <- mvn(data = sleep.wide, mvnTest = "dh")
sleep_mvn$multivariateNormality

Mardia’s MVN Test

This function performs multivariate skewness and kurtosis tests at the same time and combines test results for multivariate normality. If both tests indicate multivariate normality, then data follows a multivariate normal distribution at the 0.05 significance level.

sleep_mvn2 <- mvn(data = sleep.wide, mvnTest = "mardia")
sleep_mvn2$multivariateNormality

Check the sphericity assumption

sleep.aov$`Mauchly's Test for Sphericity`