Using contrast coding to understand ANOVA in R

In this dataset, participants are assigned to 1 of 3 conditions:

1. Control condition (recalled a neutral event)
2. Sad actor (recalled a sad event)
3. Sad partner (recalled a neutral event but was paired with the sad actor)

These data are dyadic, but for the sake of this tutorial, we’re going to ignore that and treat participants as independent.

Load the data & libraries

pkgs <- c("tidyverse", "dplyr", "haven", "foreign", "lme4", "plyr", "nlme", "lsr", "emmeans", "afex", "knitr", "kableExtra", "car")

packages <- rownames(installed.packages())
p_to_install <- pkgs[!(pkgs %in% packages)]

if(length(p_to_install) > 0){
  install.packages(p_to_install)
}

lapply(pkgs, library, character.only = TRUE)

getwd()
# Load data
df<- read_sav("lab5_anova_contrasts.sav")

---

Study Design: - Participants completed an emotion induction (recall either a sad or neutral event) and then interacted with a stranger where they asked each other questions to get acquainted.

Testable research question: - I want to know: Do people feel more (or less) sad after interacting with a sad person? Are there changes in how sad that sad individual feels after the interaction? Does this vary by gender?

Predictors: - Condition (control x sad actor x sad partner) - Gender (female x male)

Outcome: - Post-interaction sadness (1- not at all to 7- a great deal)

---

ANOVA with no data prep or pre-analyses

First off, I don’t like aov(). The summary function with aov() gives us the omnibus results with type = 1, so we’d have to use Anova from the car package anyway. That said, we’re going to use lm(). I promise, they are doing the same thing, they just present different output when you use summary().

fit <- lm(avg_sad_int ~ genderR*condition_3Level, data = df)
car::Anova(fit, type = "III")

## Anova Table (Type III tests)
## 
## Response: avg_sad_int
##                          Sum Sq  Df  F value Pr(>F)    
## (Intercept)              62.792   1 345.8616 <2e-16 ***
## genderR                   0.140   1   0.7686 0.3816    
## condition_3Level          0.086   1   0.4735 0.4921    
## genderR:condition_3Level  0.001   1   0.0032 0.9550    
## Residuals                40.668 224                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Why don’t these match?

Why are the degrees of freedom different? What does this tell us?

---

Let’s first check that the variables are in the correct class.

We’re looking at condition_3Level, genderR, and avg_sad_int.

glimpse(df)

## Rows: 230
## Columns: 6
## $ Dyad             <dbl> 100, 100, 101, 101, 102, 102, 103, 103, 105, 105, 106…
## $ ID               <dbl> 1001, 1002, 1011, 1012, 1021, 1022, 1031, 1032, 1051,…
## $ partner          <dbl> 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2,…
## $ genderR          <dbl+lbl> -1, -1, -1, -1,  1,  1, -1, -1,  1,  1,  1,  1, -…
## $ condition_3Level <dbl+lbl> 2, 3, 2, 3, 1, 1, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 1…
## $ avg_sad_int      <dbl> 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.0…

Predictors:

• condition_3Level: 1, 2, 3
• genderR: -1, 1

Change variable class

df.f <- df %>% 
  mutate(condition_3Level = as.factor(condition_3Level),
         genderR = as.factor(genderR))

# check class
class(df.f$condition_3Level)

## [1] "factor"

class(df.f$genderR)

## [1] "factor"

Now, let’s rerun the model with categorical predictors

fit2 <- lm(avg_sad_int ~ genderR*condition_3Level, data = df.f)
car::Anova(fit2, type = "III")

## Anova Table (Type III tests)
## 
## Response: avg_sad_int
##                          Sum Sq  Df  F value Pr(>F)    
## (Intercept)              73.699   1 402.5289 <2e-16 ***
## genderR                   0.375   1   2.0475 0.1539    
## condition_3Level          0.060   2   0.1629 0.8498    
## genderR:condition_3Level  0.001   2   0.0032 0.9968    
## Residuals                40.646 222                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

It still doesn’t match. What else could be going on?

What does the omnibus test tell us?

---

What is ANOVA?

The F-ratio represents how good the model is relative to how bad it is. In other words, it is the ratio of explained variance to unexplained variance.

In ANOVA, we’re asking: is predicting scores from the group means better than predicting scores from the grand mean? Are the groups’ outcomes significantly different from each other?

The equations will be variations of one basic equation:

\[deviation = \sum(observed-model)^2\]

Total Sum of Squares (SST)

This is the total amount of variation in our data. It’s the difference between the observed data point and the grand mean. Then we square these differences and add them together.

It’s the “grand variance” or the variance in responses regardless of group.

\[SST = \sum^N_{i=1}(Y_i-\overline{Y})^2\]

\[Y_i\ represents\ each \ individual\ observation. \] \[\overline{Y} \ is \ the\ grand\ mean. \] \[N\ is\ the\ total \ number\ of\ observations. \]

Between-Group Sum of Squares (SSB)

(aka model sum of squares)

How much of the variation can the model explain? In simple terms, it’s the differences between the predicted values and the grand mean.

\[SSB = \sum^k_{j=1}n_j(\overline{Y}_j-\overline{Y})^2\]

\[\overline{Y}_j\ is\ the \ mean\ of\ group\ j. \]

\[n_j \ is \ the\ number\ mean\ of\ observations\ in\ group\ j. \]

\[k\ is\ the\ total \ number\ of\ groups. \]

Within-Group Sum of Squares (SSW)

(aka residual sum of squares)

This is the variance within each group. So how much of the variance is not due to group differences, but instead due to extraneous factors like individual differences? In other words, this is the variance that cannot be explained by the model.

\[SSW =\sum^{n}_{i=1}({Y}_{ij}-\overline{Y}_j)^2\]

\[{Y}_{ij}\ represents\ each \ observation\ in\ group\ j. \]

In short, SSB tells us how much variation is due to the model and SSW tells us how much variation cannot be explained by the model (error)

Mean Squares

One issue with that formula is because the values are summed values, they will be influenced by the total number of scores summed. To get rid of this bias, we calculate the mean squares, which is the average sum of squares (SS divided by df).

Mean Squares Between

This is the average amount of variance explained by the model.

\[MSB = \frac{SSB}{df_B} \]

Mean Squares Within

This is the average amount of variance not explained by the model.

\[MSW = \frac{SSW}{df_W} \]

F-statistic

This is the ratio of variance explained by the model (systematic variance) to the variance not explained by the model (unsystematic variance)

If the F is above 1, that tells us that there was some effect of the predictor above and beyond the individual differences that could explain the outcome (but, this doesn’t tell us if it’s significant. Just the direction).

If the F is below 1, you know that there’s more error in the model than systematic variance.

\[F = \frac{MSB}{MSW} \]

tldr; the F-statistic tells us whether our model fitted to the data accounts for more variation (good) than extraneous factors (bad), but it doesn’t tell us where these differences between groups are.

For example, if we’re looking at differences in sleep quality by race and we get a significant effect, that tells us that different races do report differences in sleep quality (relative to the grand mean). However, it does not tell us which races sleep better or worse than each other.

SS Types

Imagine we’re testing the effect of various health behaviors (diet, sleep) on mood.

Type I Sums of Squares

This approach is useful if we think the sequence matters, like if we believe sleep needs to be considered before diet. Tells us the story in sequence (first sleep, then diet) and how each step adds to mood.

Type II Sum of Squares

Now, let’s say we want to understand the unique impact of diet on mood, ignoring whether or not they’re sleeping well, and vice versa. Type II sums of squares focus on the individual contribution of each factor without mixing them up. It’s like evaluating the influence of diet and sleep separately, assuming they don’t interact.

Type III Sum of Squares

Lastly, we consider the scenario where we want to understand the influence of diet and sleep on mood, including how they might interact (e.g., does the impact of diet depend on whether they slept well?). Type III sums of squares let us see the full picture, considering all possible interactions. This is best when our experiment is complex (e.g., group sizes vary, or we expect that the factors might affect each other).

I only ever use Type III.

---

Detour into contrast coding

Our main takeaway from ANOVA is that the omnibus looks at the effect of the predictor(s) relative to the grand mean (or the DV when the predictors are at their average level). With that in mind, we need to make sure our variables are actually coded in that way.

Contrast coding is a way to compare different levels/groups of a categorical variable. The most common contrasts are dummy coding and effects coding.

Dummy coding

This is the most common form of contrast coding. In dummy coding, one level of the categorical variable is chosen as a reference group. For example, if you have a variable “Race” with three levels (White, Black, Asian), you can create two dummy variables: the effect of being Black and the effect of being Asian. The reference group (White) gets coded as 0 in both new variables. The intercept represents the mean outcome for the reference group.

White = 1 Black = 2 Asian = 3

Race	Dummy1	Dummy2
White	0	0
Black	1	0
Asian	0	1

Dummy 1: outcome for Black participants vs reference (White)
Dummy 2: outcome for Asian participants vs reference (White)
Intercept = average sleep quality when predictor is zero (sleep for White participants)

Effects coding

Unlike dummy coding, where the reference category is represented by all 0’s, effects coding doesn’t leave out one category as the reference. Instead, categories are coded with numbers that balance out, such that the sum of codes for each categorical level across all coded variables equals zero. This balance allows the model’s intercept to represent the grand mean. You would use this coding scheme if you wanted to compare the mean of each group to the grand mean. The reference group gets -1 for both effects.

The sum of your contrasts should = 0.

Race	Effect1	Effect2
White	1	0
Black	0	1
Asian	-1	-1

Effect 1: compares the effect of White to the overall mean
Effect 2: compares the effect of Black to the overall mean
Intercept = grand mean across all categories of race

tldr; Dummy coding gives us the “simple effects,” comparing different levels to a reference group while effects coding gives us the “main effects,” comparing different levels to the group mean

---

Challenge!

Name the comparisons below.

1. Testing the effect of a new medication using treatment and control groups. How would you interpret the intercept? How would you interpret the effect of group?

Group	Column1
Control	0
Treatment	1

2. Now we’re testing the effect of different dosages of the same medication to the control group. How would you interpret the intercept? What does level 1 (column 1) look at? How about level 2 (column 2)?

Group	Column1	Column 2
Control	0	0
Low Dose	1	0
High Dose	0	1

3. We recoded the groups. How would you interpret the intercept? What does level 1 (column 1) look at? How about level 2 (column 2)?

Group	Column1	Column 2
Control	1	0
Low Dose	0	1
High Dose	-1	-1

4. One more…we’ve recoded the groups again. How would you interpret the intercept? What does level 1 (column 1) look at? How about level 2 (column 2)?

Group	Column1	Column 2
Control	0	2
Low Dose	1	-1
High Dose	-1	-1

---

How are our variables coded?

contrasts(df.f$genderR)

##    1
## -1 0
## 1  1

Identifiers:

• -1 = Female

• 1 = Male

Contrasts:

• 0 = Female

• 1 = Male

How does this impact our interpretation of the intercept?

• Average sadness after the interaction for women (when genderR = 0)

contrasts(df.f$condition_3Level)

##   2 3
## 1 0 0
## 2 1 0
## 3 0 1

Identifiers:

• 1 = Control dyad

• 2 = Sad actor

• 3 = Sad partner

Contrasts:

FIRST EFFECT: what is the effect of sad actor vs control condition?

SECOND EFFECT: what is the effect of sad partner vs control condition?

---

How does SPSS set up the contrasts?

Intercept

Intercept = average all levels

Gender

• female -1, but coded as 1
• male 1, but coded as -1
• when condition = 0 (average collapsing across condition)

Condition

• when gender = 0 (average collapsing across gender)
• CONTRAST 1 = control (1) vs sad partner (3)
• CONTRAST 2 = sad actor (2) vs sad partner(3)

Change reference groups

Gender contrasts

genderContrasts <- matrix(c(1,-1),
                          nrow = 2, ncol = 1, byrow = FALSE)

contrasts(df.f$genderR) = genderContrasts

print(contrasts(df.f$genderR))

##    [,1]
## -1    1
## 1    -1

Condition contrasts

effectContrasts <- matrix(c(1, 0,
                            -1, 1, 
                            -1, 0 ),
                          nrow = 3, ncol = 2, byrow = FALSE)

contrasts(df.f$condition_3Level) = effectContrasts

print(contrasts(df.f$condition_3Level))

##   [,1] [,2]
## 1    1    1
## 2    0   -1
## 3   -1    0

Rerun the omnibus test

fit3 <- lm(avg_sad_int ~ genderR*condition_3Level, data = df.f)
car::Anova(fit3, type = "III")

## Anova Table (Type III tests)
## 
## Response: avg_sad_int
##                           Sum Sq  Df   F value  Pr(>F)    
## (Intercept)              274.506   1 1499.2937 < 2e-16 ***
## genderR                    0.581   1    3.1746 0.07616 .  
## condition_3Level           0.103   2    0.2817 0.75478    
## genderR:condition_3Level   0.001   2    0.0032 0.99678    
## Residuals                 40.646 222                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Yay! It matches.

Note: technically, contrast coding shouldn't influence the omnibus test. By definition, the omnibus test is supposed to compare the group means to the grand mean. SPSS and SAS take care of this by contrast coding in the background. As you can see from the omnibus and summary statistics in R, R expects you to understand what's going on 'under the hood'

Now let’s break down these effects

summary(fit3)

## 
## Call:
## lm(formula = avg_sad_int ~ genderR * condition_3Level, data = df.f)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2738 -0.2167 -0.1576 -0.1111  2.0595 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.1972874  0.0309211  38.721   <2e-16 ***
## genderR1                   -0.0550935  0.0309211  -1.782   0.0762 .  
## condition_3Level1           0.0333985  0.0459370   0.727   0.4680    
## condition_3Level2          -0.0149933  0.0460489  -0.326   0.7450    
## genderR1:condition_3Level1 -0.0023158  0.0459370  -0.050   0.9598    
## genderR1:condition_3Level2 -0.0007076  0.0460489  -0.015   0.9878    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4279 on 222 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.02079,    Adjusted R-squared:  -0.001268 
## F-statistic: 0.9425 on 5 and 222 DF,  p-value: 0.4543

Does it match SPSS output?

No it does not. Why might that be?

How is gender coded?

How is condition coded?

Recode variables to match SPSS

Gender

genderContrasts2 <- matrix(c(1, 0), 
                           nrow = 2, ncol = 1, byrow = FALSE)
contrasts(df.f$genderR) = genderContrasts2

print(contrasts(df.f$genderR))

##    [,1]
## -1    1
## 1     0

Condition

conditionContrasts2 <- matrix(c(1, 0, 0,
                             0, 1, 0), 
                           nrow = 3, ncol = 2, byrow = FALSE)

contrasts(df.f$condition_3Level) = conditionContrasts2

print(contrasts(df.f$condition_3Level))

##   [,1] [,2]
## 1    1    0
## 2    0    1
## 3    0    0

fit4 <- lm(avg_sad_int ~ genderR*condition_3Level, data = df.f)
summary(fit4)

## 
## Call:
## lm(formula = avg_sad_int ~ genderR * condition_3Level, data = df.f)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2738 -0.2167 -0.1576 -0.1111  2.0595 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.216667   0.095679  12.716   <2e-16 ***
## genderR1                   -0.105556   0.117682  -0.897    0.371    
## condition_3Level1           0.057143   0.111463   0.513    0.609    
## condition_3Level2           0.050000   0.135311   0.370    0.712    
## genderR1:condition_3Level1 -0.010678   0.142995  -0.075    0.941    
## genderR1:condition_3Level2 -0.003216   0.166799  -0.019    0.985    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4279 on 222 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.02079,    Adjusted R-squared:  -0.001268 
## F-statistic: 0.9425 on 5 and 222 DF,  p-value: 0.4543

Woohoo! It matches.

---

Do I have to go through this process any time I want to run an ANOVA in R?

It’s important that you understand what’s going on “under the hood,” but we can shorten this process.

Let’s demonstrate the steps using our original dataframe

Step 1: Convert numeric predictors to factor variables

df.short <- df %>% 
  mutate(condition_3Level = as.factor(condition_3Level),
         genderR = as.factor(genderR))

Step 2: Contrast code for the omnibus test

We actually don’t have to create our contrasts from scratch. We can use these R functions to automatically contrast code our variables.

# effect-code
contr.sum(3)

##   [,1] [,2]
## 1    1    0
## 2    0    1
## 3   -1   -1

# dummy code
contr.treatment(3)

##   2 3
## 1 0 0
## 2 1 0
## 3 0 1

Apply contrasts (effect code)

contrasts(df.short$genderR) <- contr.sum(2)
contrasts(df.short$condition_3Level) <- contr.sum(3)

Step 3a: Run omnibus test

model_shortcut <- lm(avg_sad_int ~ genderR*condition_3Level, data = df.short)

car::Anova(model_shortcut, type = "III")

## Anova Table (Type III tests)
## 
## Response: avg_sad_int
##                           Sum Sq  Df   F value  Pr(>F)    
## (Intercept)              274.506   1 1499.2937 < 2e-16 ***
## genderR                    0.581   1    3.1746 0.07616 .  
## condition_3Level           0.103   2    0.2817 0.75478    
## genderR:condition_3Level   0.001   2    0.0032 0.99678    
## Residuals                 40.646 222                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Apply contrasts (dummy code)

contrasts(df.short$genderR) <- contr.treatment(2)
contrasts(df.short$condition_3Level) <- contr.treatment(3)

Step 4a: Run regression to get simple effects

Be mindful that this approach defaults to a different reference group than SPSS (notice Condition =1 is missing from the R output and Condition =3 is missing from the SPSS output). It’s not wrong, but good to know that their default settings are different.

model_shortcut2 <- lm(avg_sad_int ~ genderR*condition_3Level, data = df.short)

summary(model_shortcut2)

## 
## Call:
## lm(formula = avg_sad_int ~ genderR * condition_3Level, data = df.short)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2738 -0.2167 -0.1576 -0.1111  2.0595 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.157576   0.057697  20.063   <2e-16 ***
## genderR2                    0.116234   0.081230   1.431    0.154    
## condition_3Level2           0.000319   0.090261   0.004    0.997    
## condition_3Level3          -0.046465   0.089574  -0.519    0.604    
## genderR2:condition_3Level2 -0.007462   0.143426  -0.052    0.959    
## genderR2:condition_3Level3 -0.010678   0.142995  -0.075    0.941    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4279 on 222 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.02079,    Adjusted R-squared:  -0.001268 
## F-statistic: 0.9425 on 5 and 222 DF,  p-value: 0.4543

Step 3b: Change reference group to match SPSS

# base is reference (it doesn't take negative numbers, but we're telling it we want the reference to be the higher number of the two)
contrasts(df.short$genderR) <- contr.treatment(2, base = 2)

contrasts(df.short$condition_3Level) <- contr.treatment(3, base = 3)

Step 4b: Run regression to get simple effects to match SPSS

model_shortcut3 <- lm(avg_sad_int ~ genderR*condition_3Level, data = df.short)

summary(model_shortcut3)

## 
## Call:
## lm(formula = avg_sad_int ~ genderR * condition_3Level, data = df.short)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2738 -0.2167 -0.1576 -0.1111  2.0595 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.216667   0.095679  12.716   <2e-16 ***
## genderR1                   -0.105556   0.117682  -0.897    0.371    
## condition_3Level1           0.057143   0.111463   0.513    0.609    
## condition_3Level2           0.050000   0.135311   0.370    0.712    
## genderR1:condition_3Level1 -0.010678   0.142995  -0.075    0.941    
## genderR1:condition_3Level2 -0.003216   0.166799  -0.019    0.985    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4279 on 222 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.02079,    Adjusted R-squared:  -0.001268 
## F-statistic: 0.9425 on 5 and 222 DF,  p-value: 0.4543

Do I really have to do all that just for ANOVA?

No, but it’s important that you understand how ANOVA works. Here’s an actual shortcut that defaults to SS type 3 and does not require you to contrast code. You can run the pairwise comparisons as you normally would.

Omnibus shortcut

# we don't even need to convert numeric predictors to factor. It does that for us.

# omnibus test
shortcut_model <- afex::aov_car(avg_sad_int ~ genderR*condition_3Level + 
                Error(ID), # must add the error term (i.e. ID)
              data = df) # original dataframe

## Converting to factor: genderR, condition_3Level

## Contrasts set to contr.sum for the following variables: genderR, condition_3Level

print(shortcut_model)

## Anova Table (Type 3 tests)
## 
## Response: avg_sad_int
##                     Effect     df  MSE      F   ges p.value
## 1                  genderR 1, 222 0.18 3.17 +  .014    .076
## 2         condition_3Level 2, 222 0.18   0.28  .003    .755
## 3 genderR:condition_3Level 2, 222 0.18   0.00 <.001    .997
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

---

Pairwise comparisons to breakdown interactions

emm <- emmeans(model_shortcut3, specs = c("genderR", "condition_3Level"))

pairs(emm, simple = "each")

## $`simple contrasts for genderR`
## condition_3Level = 1:
##  contrast               estimate     SE  df t.ratio p.value
##  (genderR-1) - genderR1   -0.116 0.0812 222  -1.431  0.1539
## 
## condition_3Level = 2:
##  contrast               estimate     SE  df t.ratio p.value
##  (genderR-1) - genderR1   -0.109 0.1182 222  -0.920  0.3585
## 
## condition_3Level = 3:
##  contrast               estimate     SE  df t.ratio p.value
##  (genderR-1) - genderR1   -0.106 0.1177 222  -0.897  0.3707
## 
## 
## $`simple contrasts for condition_3Level`
## genderR = -1:
##  contrast                               estimate     SE  df t.ratio p.value
##  condition_3Level1 - condition_3Level2 -0.000319 0.0903 222  -0.004  1.0000
##  condition_3Level1 - condition_3Level3  0.046465 0.0896 222   0.519  0.8623
##  condition_3Level2 - condition_3Level3  0.046784 0.0975 222   0.480  0.8810
## 
## genderR = 1:
##  contrast                               estimate     SE  df t.ratio p.value
##  condition_3Level1 - condition_3Level2  0.007143 0.1115 222   0.064  0.9977
##  condition_3Level1 - condition_3Level3  0.057143 0.1115 222   0.513  0.8653
##  condition_3Level2 - condition_3Level3  0.050000 0.1353 222   0.370  0.9275
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

pairs(emm, adjust = "tukey")

##  contrast                                                       estimate     SE
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level1    -0.116234 0.0812
##  (genderR-1 condition_3Level1) - (genderR-1 condition_3Level2) -0.000319 0.0903
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level2    -0.109091 0.1117
##  (genderR-1 condition_3Level1) - (genderR-1 condition_3Level3)  0.046465 0.0896
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level3    -0.059091 0.1117
##  genderR1 condition_3Level1 - (genderR-1 condition_3Level2)     0.115915 0.0899
##  genderR1 condition_3Level1 - genderR1 condition_3Level2        0.007143 0.1115
##  genderR1 condition_3Level1 - (genderR-1 condition_3Level3)     0.162698 0.0892
##  genderR1 condition_3Level1 - genderR1 condition_3Level3        0.057143 0.1115
##  (genderR-1 condition_3Level2) - genderR1 condition_3Level2    -0.108772 0.1182
##  (genderR-1 condition_3Level2) - (genderR-1 condition_3Level3)  0.046784 0.0975
##  (genderR-1 condition_3Level2) - genderR1 condition_3Level3    -0.058772 0.1182
##  genderR1 condition_3Level2 - (genderR-1 condition_3Level3)     0.155556 0.1177
##  genderR1 condition_3Level2 - genderR1 condition_3Level3        0.050000 0.1353
##  (genderR-1 condition_3Level3) - genderR1 condition_3Level3    -0.105556 0.1177
##   df t.ratio p.value
##  222  -1.431  0.7082
##  222  -0.004  1.0000
##  222  -0.976  0.9250
##  222   0.519  0.9954
##  222  -0.529  0.9950
##  222   1.289  0.7909
##  222   0.064  1.0000
##  222   1.823  0.4530
##  222   0.513  0.9957
##  222  -0.920  0.9410
##  222   0.480  0.9968
##  222  -0.497  0.9962
##  222   1.322  0.7727
##  222   0.370  0.9991
##  222  -0.897  0.9469
## 
## P value adjustment: tukey method for comparing a family of 6 estimates

Alternative formatting

kable(pairs(emm, adjust = "tukey"), "simple", digits = 3) 

contrast	estimate	SE	df	t.ratio	p.value
(genderR-1 condition_3Level1) - genderR1 condition_3Level1	-0.116	0.081	222	-1.431	0.708
(genderR-1 condition_3Level1) - (genderR-1 condition_3Level2)	0.000	0.090	222	-0.004	1.000
(genderR-1 condition_3Level1) - genderR1 condition_3Level2	-0.109	0.112	222	-0.976	0.925
(genderR-1 condition_3Level1) - (genderR-1 condition_3Level3)	0.046	0.090	222	0.519	0.995
(genderR-1 condition_3Level1) - genderR1 condition_3Level3	-0.059	0.112	222	-0.529	0.995
genderR1 condition_3Level1 - (genderR-1 condition_3Level2)	0.116	0.090	222	1.289	0.791
genderR1 condition_3Level1 - genderR1 condition_3Level2	0.007	0.111	222	0.064	1.000
genderR1 condition_3Level1 - (genderR-1 condition_3Level3)	0.163	0.089	222	1.823	0.453
genderR1 condition_3Level1 - genderR1 condition_3Level3	0.057	0.111	222	0.513	0.996
(genderR-1 condition_3Level2) - genderR1 condition_3Level2	-0.109	0.118	222	-0.920	0.941
(genderR-1 condition_3Level2) - (genderR-1 condition_3Level3)	0.047	0.098	222	0.480	0.997
(genderR-1 condition_3Level2) - genderR1 condition_3Level3	-0.059	0.118	222	-0.497	0.996
genderR1 condition_3Level2 - (genderR-1 condition_3Level3)	0.156	0.118	222	1.322	0.773
genderR1 condition_3Level2 - genderR1 condition_3Level3	0.050	0.135	222	0.370	0.999
(genderR-1 condition_3Level3) - genderR1 condition_3Level3	-0.106	0.118	222	-0.897	0.947

kable(pairs(emm, adjust = "bonferroni"), "simple", digits = 3) 

contrast	estimate	SE	df	t.ratio	p.value
(genderR-1 condition_3Level1) - genderR1 condition_3Level1	-0.116	0.081	222	-1.431	1
(genderR-1 condition_3Level1) - (genderR-1 condition_3Level2)	0.000	0.090	222	-0.004	1
(genderR-1 condition_3Level1) - genderR1 condition_3Level2	-0.109	0.112	222	-0.976	1
(genderR-1 condition_3Level1) - (genderR-1 condition_3Level3)	0.046	0.090	222	0.519	1
(genderR-1 condition_3Level1) - genderR1 condition_3Level3	-0.059	0.112	222	-0.529	1
genderR1 condition_3Level1 - (genderR-1 condition_3Level2)	0.116	0.090	222	1.289	1
genderR1 condition_3Level1 - genderR1 condition_3Level2	0.007	0.111	222	0.064	1
genderR1 condition_3Level1 - (genderR-1 condition_3Level3)	0.163	0.089	222	1.823	1
genderR1 condition_3Level1 - genderR1 condition_3Level3	0.057	0.111	222	0.513	1
(genderR-1 condition_3Level2) - genderR1 condition_3Level2	-0.109	0.118	222	-0.920	1
(genderR-1 condition_3Level2) - (genderR-1 condition_3Level3)	0.047	0.098	222	0.480	1
(genderR-1 condition_3Level2) - genderR1 condition_3Level3	-0.059	0.118	222	-0.497	1
genderR1 condition_3Level2 - (genderR-1 condition_3Level3)	0.156	0.118	222	1.322	1
genderR1 condition_3Level2 - genderR1 condition_3Level3	0.050	0.135	222	0.370	1
(genderR-1 condition_3Level3) - genderR1 condition_3Level3	-0.106	0.118	222	-0.897	1

What else can we do with emmeans()?

Estimated marginal means

emm2 <- emmeans(model_shortcut3, specs = pairwise ~ genderR:condition_3Level)

emm2$emmeans

##  genderR condition_3Level emmean     SE  df lower.CL upper.CL
##  -1      1                  1.16 0.0577 222    1.044     1.27
##  1       1                  1.27 0.0572 222    1.161     1.39
##  -1      2                  1.16 0.0694 222    1.021     1.29
##  1       2                  1.27 0.0957 222    1.078     1.46
##  -1      3                  1.11 0.0685 222    0.976     1.25
##  1       3                  1.22 0.0957 222    1.028     1.41
## 
## Confidence level used: 0.95

Contrasts

emm2$contrasts

##  contrast                                                       estimate     SE
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level1    -0.116234 0.0812
##  (genderR-1 condition_3Level1) - (genderR-1 condition_3Level2) -0.000319 0.0903
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level2    -0.109091 0.1117
##  (genderR-1 condition_3Level1) - (genderR-1 condition_3Level3)  0.046465 0.0896
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level3    -0.059091 0.1117
##  genderR1 condition_3Level1 - (genderR-1 condition_3Level2)     0.115915 0.0899
##  genderR1 condition_3Level1 - genderR1 condition_3Level2        0.007143 0.1115
##  genderR1 condition_3Level1 - (genderR-1 condition_3Level3)     0.162698 0.0892
##  genderR1 condition_3Level1 - genderR1 condition_3Level3        0.057143 0.1115
##  (genderR-1 condition_3Level2) - genderR1 condition_3Level2    -0.108772 0.1182
##  (genderR-1 condition_3Level2) - (genderR-1 condition_3Level3)  0.046784 0.0975
##  (genderR-1 condition_3Level2) - genderR1 condition_3Level3    -0.058772 0.1182
##  genderR1 condition_3Level2 - (genderR-1 condition_3Level3)     0.155556 0.1177
##  genderR1 condition_3Level2 - genderR1 condition_3Level3        0.050000 0.1353
##  (genderR-1 condition_3Level3) - genderR1 condition_3Level3    -0.105556 0.1177
##   df t.ratio p.value
##  222  -1.431  0.7082
##  222  -0.004  1.0000
##  222  -0.976  0.9250
##  222   0.519  0.9954
##  222  -0.529  0.9950
##  222   1.289  0.7909
##  222   0.064  1.0000
##  222   1.823  0.4530
##  222   0.513  0.9957
##  222  -0.920  0.9410
##  222   0.480  0.9968
##  222  -0.497  0.9962
##  222   1.322  0.7727
##  222   0.370  0.9991
##  222  -0.897  0.9469
## 
## P value adjustment: tukey method for comparing a family of 6 estimates

Confidence intervals

emm2$contrasts %>% confint()

##  contrast                                                       estimate     SE
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level1    -0.116234 0.0812
##  (genderR-1 condition_3Level1) - (genderR-1 condition_3Level2) -0.000319 0.0903
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level2    -0.109091 0.1117
##  (genderR-1 condition_3Level1) - (genderR-1 condition_3Level3)  0.046465 0.0896
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level3    -0.059091 0.1117
##  genderR1 condition_3Level1 - (genderR-1 condition_3Level2)     0.115915 0.0899
##  genderR1 condition_3Level1 - genderR1 condition_3Level2        0.007143 0.1115
##  genderR1 condition_3Level1 - (genderR-1 condition_3Level3)     0.162698 0.0892
##  genderR1 condition_3Level1 - genderR1 condition_3Level3        0.057143 0.1115
##  (genderR-1 condition_3Level2) - genderR1 condition_3Level2    -0.108772 0.1182
##  (genderR-1 condition_3Level2) - (genderR-1 condition_3Level3)  0.046784 0.0975
##  (genderR-1 condition_3Level2) - genderR1 condition_3Level3    -0.058772 0.1182
##  genderR1 condition_3Level2 - (genderR-1 condition_3Level3)     0.155556 0.1177
##  genderR1 condition_3Level2 - genderR1 condition_3Level3        0.050000 0.1353
##  (genderR-1 condition_3Level3) - genderR1 condition_3Level3    -0.105556 0.1177
##   df lower.CL upper.CL
##  222  -0.3498    0.117
##  222  -0.2598    0.259
##  222  -0.4303    0.212
##  222  -0.2110    0.304
##  222  -0.3803    0.262
##  222  -0.1426    0.374
##  222  -0.3133    0.328
##  222  -0.0939    0.419
##  222  -0.2633    0.378
##  222  -0.4486    0.231
##  222  -0.2336    0.327
##  222  -0.3986    0.281
##  222  -0.1828    0.494
##  222  -0.3390    0.439
##  222  -0.4439    0.233
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 6 estimates

Within-group comparisons with |

emm3 <- emmeans(model_shortcut3, specs = pairwise ~ genderR|condition_3Level)

emm3 

## $emmeans
## condition_3Level = 1:
##  genderR emmean     SE  df lower.CL upper.CL
##  -1        1.16 0.0577 222    1.044     1.27
##  1         1.27 0.0572 222    1.161     1.39
## 
## condition_3Level = 2:
##  genderR emmean     SE  df lower.CL upper.CL
##  -1        1.16 0.0694 222    1.021     1.29
##  1         1.27 0.0957 222    1.078     1.46
## 
## condition_3Level = 3:
##  genderR emmean     SE  df lower.CL upper.CL
##  -1        1.11 0.0685 222    0.976     1.25
##  1         1.22 0.0957 222    1.028     1.41
## 
## Confidence level used: 0.95 
## 
## $contrasts
## condition_3Level = 1:
##  contrast               estimate     SE  df t.ratio p.value
##  (genderR-1) - genderR1   -0.116 0.0812 222  -1.431  0.1539
## 
## condition_3Level = 2:
##  contrast               estimate     SE  df t.ratio p.value
##  (genderR-1) - genderR1   -0.109 0.1182 222  -0.920  0.3585
## 
## condition_3Level = 3:
##  contrast               estimate     SE  df t.ratio p.value
##  (genderR-1) - genderR1   -0.106 0.1177 222  -0.897  0.3707

Pass these contrasts through rbind() to correct for multiple comparisons (defaults to Bonferroni, which can be a bit too conservative and too statistically underpowered for so many comparisons; see Lee & Lee, 2018; Lieberman & Cunningham, 2009)

emm2$contrasts %>% 
  rbind()

##  contrast                                                       estimate     SE
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level1    -0.116234 0.0812
##  (genderR-1 condition_3Level1) - (genderR-1 condition_3Level2) -0.000319 0.0903
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level2    -0.109091 0.1117
##  (genderR-1 condition_3Level1) - (genderR-1 condition_3Level3)  0.046465 0.0896
##  (genderR-1 condition_3Level1) - genderR1 condition_3Level3    -0.059091 0.1117
##  genderR1 condition_3Level1 - (genderR-1 condition_3Level2)     0.115915 0.0899
##  genderR1 condition_3Level1 - genderR1 condition_3Level2        0.007143 0.1115
##  genderR1 condition_3Level1 - (genderR-1 condition_3Level3)     0.162698 0.0892
##  genderR1 condition_3Level1 - genderR1 condition_3Level3        0.057143 0.1115
##  (genderR-1 condition_3Level2) - genderR1 condition_3Level2    -0.108772 0.1182
##  (genderR-1 condition_3Level2) - (genderR-1 condition_3Level3)  0.046784 0.0975
##  (genderR-1 condition_3Level2) - genderR1 condition_3Level3    -0.058772 0.1182
##  genderR1 condition_3Level2 - (genderR-1 condition_3Level3)     0.155556 0.1177
##  genderR1 condition_3Level2 - genderR1 condition_3Level3        0.050000 0.1353
##  (genderR-1 condition_3Level3) - genderR1 condition_3Level3    -0.105556 0.1177
##   df t.ratio p.value
##  222  -1.431  1.0000
##  222  -0.004  1.0000
##  222  -0.976  1.0000
##  222   0.519  1.0000
##  222  -0.529  1.0000
##  222   1.289  1.0000
##  222   0.064  1.0000
##  222   1.823  1.0000
##  222   0.513  1.0000
##  222  -0.920  1.0000
##  222   0.480  1.0000
##  222  -0.497  1.0000
##  222   1.322  1.0000
##  222   0.370  1.0000
##  222  -0.897  1.0000
## 
## P value adjustment: bonferroni method for 15 tests

---

Mini Check-in

1. What is a one-way ANOVA? How is it different from a factorial or two-way ANOVA?
2. What is an interaction? How would you interpret an interaction?
3. What kind of graph could I use to plot a one-way ANOVA?

---

Replicate Madalina’s Python code

m_df_og <- readxl::read_xlsx("data_ANOVA.xlsx")

m_df <- m_df_og %>% 
  group_by(partnum) %>% 
  dplyr::summarise(across(everything(), ~ mean(., na.rm = TRUE)) )%>% 
  ungroup()

head(m_df)

## # A tibble: 6 × 11
##   partnum RdeltaB pop_unpop sci_anec convincing  cond itemnum   age party
##     <dbl>   <dbl>     <dbl>    <dbl>      <dbl> <dbl>   <dbl> <dbl> <dbl>
## 1       0    1.91     0.5      0.531       46.1  1.56    15.5    18     1
## 2       1    2.09     0.469    0.375       53.4  1.22    15.5    19     2
## 3       2    5.34     0.656    0.469       48.0  1.59    15.5    18     1
## 4       3   -7.78     0.469    0.344        6    1.16    15.5    22     1
## 5       4    6.09     0.594    0.375       48.4  1.34    15.5    19     3
## 6       5   -1.34     0.469    0.469       46.1  1.41    15.5    20     3
## # ℹ 2 more variables: twitter <dbl>, trust <dbl>

One-way ANOVA

oneway_model <- afex::aov_car(RdeltaB ~ party + Error(partnum), 
              data = m_df) # original dataframe

## Converting to factor: party

## Contrasts set to contr.sum for the following variables: party

print(oneway_model)

## Anova Table (Type 3 tests)
## 
## Response: RdeltaB
##   Effect     df   MSE    F  ges p.value
## 1  party 2, 197 18.07 1.60 .016    .204
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Two-way ANOVA

# change to type = 2 to match python
twoway_model <- afex::aov_car(RdeltaB ~ party*twitter + Error(partnum),
                              type = "II",
                              data = m_df) 

## Converting to factor: party, twitter

## Contrasts set to contr.sum for the following variables: party, twitter

print(twoway_model)

## Anova Table (Type II tests)
## 
## Response: RdeltaB
##          Effect     df   MSE    F  ges p.value
## 1         party 2, 194 18.16 1.62 .016    .200
## 2       twitter 1, 194 18.16 0.49 .003    .486
## 3 party:twitter 2, 194 18.16 0.77 .008    .465
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Repeated Measures ANOVA

Here, I’m showing you two different ways to run a repeated measures ANOVA with the afex package. The first is the aov_car function, which takes a lm-style formula.

# within the Error statement, we're saying that within each partnum (ID), sci_anec and pop_unpop are repeated.

# remember this 'nesting' format. It will come up again when we run mixed models.

repeated_model1 <- afex::aov_car(RdeltaB ~ 1 + Error(partnum/(sci_anec*pop_unpop)),
                             data = m_df_og)

knitr::kable(nice(repeated_model1))

Effect	df	MSE	F	ges	p.value
sci_anec	1, 199	48.66	6.82 **	.008	.010
pop_unpop	1, 199	43.98	10.27 **	.011	.002
sci_anec:pop_unpop	1, 199	31.68	0.93	<.001	.336

Now, we’re going to use a different function from the afex package called aov_ez. The only difference is that aov_ez takes character arguments (like you literally enter what your between and within variables are) instead of a formula.

repeated_model2 <- afex::aov_ez(id = "partnum", 
                             dv = "RdeltaB", 
                             within = c('sci_anec', 'pop_unpop'),
                             data = m_df_og) 

knitr::kable(nice(repeated_model2))

Effect	df	MSE	F	ges	p.value
sci_anec	1, 199	48.66	6.82 **	.008	.010
pop_unpop	1, 199	43.98	10.27 **	.011	.002
sci_anec:pop_unpop	1, 199	31.68	0.93	<.001	.336

summary(repeated_model2)

## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                    Sum Sq num Df Error SS den Df F value    Pr(>F)    
## (Intercept)        4580.5      1  15779.3    199 57.7665 1.136e-12 ***
## sci_anec            331.7      1   9683.2    199  6.8166  0.009719 ** 
## pop_unpop           451.8      1   8751.2    199 10.2742  0.001572 ** 
## sci_anec:pop_unpop   29.5      1   6304.4    199  0.9309  0.335806    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

What about Mauchly’s test of sphericity?

Mauchly’s test of sphericity is only applicable when you have more than two levels in at least one of your within-participant (repeated measures) factors. This is because sphericity is a concept that applies to the variances of the differences between all combinations of levels within a factor. If a factor only has two levels, the assumption of sphericity is inherently met because there’s only one possible difference to consider, thus no variability in these differences across levels to test against.

library(ez)

mod<- ezANOVA(data = m_df_og, dv = .(RdeltaB), wid = .(partnum),
              within = .(sci_anec,pop_unpop),
              detailed = TRUE,
              type = 3)

mod

## $ANOVA
##               Effect DFn DFd       SSn      SSd          F           p p<.05
## 1           sci_anec   1 199 331.69405 9683.242  6.8166339 0.009718873     *
## 2          pop_unpop   1 199 451.81686 8751.161 10.2742429 0.001571558     *
## 3 sci_anec:pop_unpop   1 199  29.49055 6304.404  0.9308761 0.335806365      
##           ges
## 1 0.013230452
## 2 0.017935914
## 3 0.001190657

---

ANCOVA

ANOVA can be extended to include one or more continuous variables that predict the outcome. Continuous variables such as these, that are not part of the main experimental manipulation but have an influence on the dependent variable, are known as covariates and they can be included in an ANOVA analysis.

Why use ANCOVA?

To reduce within-group error variance: In the discussion of ANOVA and t-tests we got used to the idea that we assess the effect of an experiment by comparing the amount of variability in the data that the experiment can explain against the variability that it cannot explain. If we can explain some of this ‘unexplained’ variance (SSW) in terms of other variables (covariates), then we reduce the error variance, allowing us to more accurately assess the effect of the independent variable (SSB).

Elimination of confounds: In any experiment, there may be unmeasured variables that confound the results (i.e., variables other than the experimental manipulation that affect the outcome variable). If any variables are known to influence the dependent variable being measured, then ANCOVA is ideally suited to remove the bias of these variables. Once a possible confounding variable has been identified, it can be measured and entered into the analysis as a covariate.

Assumptions of the ANCOVA: (1) independence of the covariate and treatment effect and (2) homogeneity of regression slopes.

1. Basically, we want to make sure that the covariate and the predictor are not related. In other words, you should have three variances: the variance that is explained by your predictor, the variance explained by the covariate, and unexplained variance. The variance explained by the covariate should not overlap with the variance explained by the predictor. When your groups differ on the covariate, putting it into the analysis will not ‘balance out the differences.’ (e.g., anxious vs non-anxious will likely differ in depression, so controlling for depression isn’t going to give you a ‘pure’ effect of anxiety).

2. We assume that the relationship between the covariate and the outcome is the same across all groups because the slope is supposed to represent the average across groups. It assumes a linear slope. If the relationship is non-linear, you would need a different model.

Example with Viagra data

credit: Andy Field’s Discovering Statistics

Here, we want to know the effect of Viagra (none, low-dose, high-dose) on libido. However, we know that partner’s libido might impact our own libido. Therefore, we’ll need to adjust for partner’s libido.

v_df <- read.delim("ViagraCovariate.dat")

v_df$dose <- as.factor(v_df$dose)

Inspect the data

Test the independence of the covariate and predictor

contrasts(v_df$dose) <- contr.sum(3)

# is there a relationship between dose and partner libido?
test1 <- lm(partnerLibido ~ dose, data = v_df)

car::Anova(test1, type = "III")

## Anova Table (Type III tests)
## 
## Response: partnerLibido
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 234.592  1 72.7232 3.846e-09 ***
## dose         12.769  2  1.9793    0.1577    
## Residuals    87.097 27                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check the slopes

v_df %>% 
  ggplot(aes(x=partnerLibido, y=libido, color=dose)) + # Map partnerLibido to x, libido to y, and dose to color
  geom_point() +  # Use geom_point for a scatter plot
  geom_smooth(method = "lm")

Homogeneity of variances

car::leveneTest(v_df$libido, v_df$dose)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.3256 0.7249
##       27

Run the ANCOVA

v_model <- lm(libido ~ dose + partnerLibido, data = v_df)

car::Anova(v_model, type = "III")

## Anova Table (Type III tests)
## 
## Response: libido
##               Sum Sq Df F value    Pr(>F)    
## (Intercept)   76.069  1 25.0205 3.342e-05 ***
## dose          25.185  2  4.1419   0.02745 *  
## partnerLibido 15.076  1  4.9587   0.03483 *  
## Residuals     79.047 26                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

v_emm <- emmeans(v_model, specs = ~ dose)

print(v_emm)

##  dose emmean    SE df lower.CL upper.CL
##  1      2.93 0.596 26     1.70     4.15
##  2      4.71 0.621 26     3.44     5.99
##  3      5.15 0.503 26     4.12     6.18
## 
## Confidence level used: 0.95

contrasts(v_df$dose) <- contr.treatment(3)

# rerunning the model as dummy coded
v_model <- lm(libido ~ dose + partnerLibido, data = v_df)
summary(v_model)

## 
## Call:
## lm(formula = libido ~ dose + partnerLibido, data = v_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2622 -0.7899 -0.3230  0.8811  4.5699 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     1.7892     0.8671   2.063   0.0492 *
## dose2           1.7857     0.8494   2.102   0.0454 *
## dose3           2.2249     0.8028   2.771   0.0102 *
## partnerLibido   0.4160     0.1868   2.227   0.0348 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.744 on 26 degrees of freedom
## Multiple R-squared:  0.2876, Adjusted R-squared:  0.2055 
## F-statistic:   3.5 on 3 and 26 DF,  p-value: 0.02954

What if I want to compare placebo vs low + high dose and low vs high?

contrasts(v_df$dose) <- cbind(c(-2,1,1), c(0,-1,1))

print(contrasts(v_df$dose)) #placebo, low dose, then high dose

##   [,1] [,2]
## 1   -2    0
## 2    1   -1
## 3    1    1

dose1: placebo vs treatment (low + high dose)
dose2: low vs high
How would you interpret partnerLibido?

v_model2 <- lm(libido ~ dose + partnerLibido, data = v_df)

summary(v_model2)

## 
## Call:
## lm(formula = libido ~ dose + partnerLibido, data = v_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2622 -0.7899 -0.3230  0.8811  4.5699 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     3.1260     0.6250   5.002 3.34e-05 ***
## dose1           0.6684     0.2400   2.785  0.00985 ** 
## dose2           0.2196     0.4056   0.541  0.59284    
## partnerLibido   0.4160     0.1868   2.227  0.03483 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.744 on 26 degrees of freedom
## Multiple R-squared:  0.2876, Adjusted R-squared:  0.2055 
## F-statistic:   3.5 on 3 and 26 DF,  p-value: 0.02954

Post hoc tests

emm_v <- emmeans(v_model2, specs = pairwise ~ dose + partnerLibido, adjust = "Tukey")

emm_v$contrasts

##  contrast                                                                 
##  dose1 partnerLibido2.73333333333333 - dose2 partnerLibido2.73333333333333
##  dose1 partnerLibido2.73333333333333 - dose3 partnerLibido2.73333333333333
##  dose2 partnerLibido2.73333333333333 - dose3 partnerLibido2.73333333333333
##  estimate    SE df t.ratio p.value
##    -1.786 0.849 26  -2.102  0.1089
##    -2.225 0.803 26  -2.771  0.0266
##    -0.439 0.811 26  -0.541  0.8517
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Effect size

What’s the difference between eta squared and partial eta squared?

\[\eta^2 = \frac{SS_{effect}}{SS_{total}}\] Eta-squared represents the total variance explained by the effect, considering both that effect and any of its interactions.

\[ partial \ \eta^2 = \frac{SS_{effect}}{SS_{effect}+SS_{residual}}\] Partial eta-squared is the total variance explained by the effect that is not explained by other variables in the model

lsr::etaSquared(v_model2)

##                  eta.sq eta.sq.part
## dose          0.2269618   0.2416256
## partnerLibido 0.1358583   0.1601709

Contrasts effect sizes

\[ r_{contrast} = \sqrt \frac{t^2}{t^2 +df}\]

rcontrast <- function(t, df) {
  r <- sqrt(t^2/(t^2 + df))
  print(paste("r = ", r))
}

t <- c(2.785, 0.541, 2.227)
df <- 26

rcontrast(t, df)

## [1] "r =  0.47934506773475"  "r =  0.105506648839546" "r =  0.400242354364157"

The output shows that the effect of the covariate (.400) and the difference between the combined dose groups and the placebo (.479) both represent medium to large effect sizes (they’re both between .4 and .5). The difference between the high- and low-dose groups (.106) was a fairly small effect.

---

Testing the homogeneity of regression slopes

v_model_int <- lm(libido ~ dose + partnerLibido + dose:partnerLibido, #interaction term
                  data = v_df)

car::Anova(v_model_int, type = "3")

## Anova Table (Type III tests)
## 
## Response: libido
##                    Sum Sq Df F value    Pr(>F)    
## (Intercept)        53.542  1 21.9207 9.323e-05 ***
## dose               36.558  2  7.4836   0.00298 ** 
## partnerLibido      17.182  1  7.0346   0.01395 *  
## dose:partnerLibido 20.427  2  4.1815   0.02767 *  
## Residuals          58.621 24                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Uh oh… we violated that assumption. What now?

---

Your turn!

Make sure to inspect the data first for any assumption violations!

Q1: injury by superhero costume

1. Superhero.dat : Children reporting to the emergency room had the severity of their injury (injury) assessed (on a scale from 0, no injury, to 100, death). In addition, a note was taken of which superhero costume they were wearing (hero): Spiderman, Superman, the Hulk or a Teenage Mutant Ninja Turtle. Test the hypotheses that different costumes are associated with more severe injuries. Make sure to run the levene’s test and compare the different costumes!

load data

sh_df <- read.delim("/Users/kareenadelrosario/Downloads/Superhero.dat")

glimpse(sh_df)

## Rows: 30
## Columns: 2
## $ hero   <int> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, …
## $ injury <int> 51, 31, 58, 20, 47, 37, 49, 40, 69, 32, 85, 66, 58, 52, 26, 43,…

test for homogeneity of variance

# Run a levene's test
leveneTest(sh_df$injury, sh_df$hero)

## Warning in leveneTest.default(sh_df$injury, sh_df$hero): sh_df$hero coerced to
## factor.

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   0.827  0.491
##       26

run one-way ANOVA

Note that here, I’m taking the manual approach, but you can also run the afex code.

# make sure the predictor is coded correctly
sh_df$hero <- as.factor(sh_df$hero)
contrasts(sh_df$hero) <- contr.sum(4)

superModel <- lm(injury~hero, data = sh_df)
Anova(superModel, type = "3")

## Anova Table (Type III tests)
## 
## Response: injury
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept)  49402  1 294.8311 1.048e-15 ***
## hero          4181  3   8.3166 0.0004828 ***
## Residuals     4357 26                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

afex option

if you decide to use aov_car, note that you’ll need an ID variable

# create id variable for error term
sh_df_afex <- sh_df %>% 
  dplyr::mutate(id = row_number())

# make sure the predictor is coded correctly
afex::aov_car(injury~hero + Error(id), data = sh_df_afex)

## Contrasts set to contr.sum for the following variables: hero

## Anova Table (Type 3 tests)
## 
## Response: injury
##   Effect    df    MSE        F  ges p.value
## 1   hero 3, 26 167.56 8.32 *** .490   <.001
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Compare injury by costume

sh_emm <- emmeans(superModel, specs = ~ hero)

print(sh_emm)

##  hero emmean   SE df lower.CL upper.CL
##  1      41.6 4.58 26     32.2     51.0
##  2      60.3 5.28 26     49.5     71.2
##  3      35.4 4.58 26     26.0     44.8
##  4      26.2 4.58 26     16.8     35.7
## 
## Confidence level used: 0.95

Bonferroni pairwise comparisons

pairs(sh_emm) %>% rbind()

##  contrast      estimate   SE df t.ratio p.value
##  hero1 - hero2   -18.71 6.99 26  -2.676  0.0763
##  hero1 - hero3     6.25 6.47 26   0.966  1.0000
##  hero1 - hero4    15.38 6.47 26   2.376  0.1511
##  hero2 - hero3    24.96 6.99 26   3.570  0.0085
##  hero2 - hero4    34.08 6.99 26   4.875  0.0003
##  hero3 - hero4     9.12 6.47 26   1.410  1.0000
## 
## P value adjustment: bonferroni method for 6 tests

---

Q2: Hangover cure?

2. HangoverCure.dat : A marketing manager for a certain well-known drinks manufacturer was interested in the therapeutic benefit of certain soft drinks for curing hangovers. He took 15 people out on the town one night and got them drunk. The next morning as they awoke, dehydrated and feeling as though they’d licked a camel’s sandy feet clean with their tongue, he gave five of them water to drink, five of them Lucozade (in case this isn’t sold outside the UK, it’s a very nice glucose-based drink) and the remaining five a leading brand of cola (this variable is called drink). He then measured how well they felt (on a scale from 0 = I feel like death to 10 = I feel really full of beans and healthy) two hours later (this variable is called well). He wanted to know which drink produced the greatest level of wellness. However, he realized it was important to control for how drunk the person got the night before, and so he measured this on a scale from 0 = as sober as a nun to 10 = flapping about like a haddock out of water on the floor in a puddle of their own vomit. The data are in the file HangoverCure. dat. Test whether people felt better after different drinks when controlling for how drunk they were the night before.

read data

hc_df <- read.delim("/Users/kareenadelrosario/Downloads/HangoverCure.dat", header = T)

glimpse(hc_df)

## Rows: 15
## Columns: 3
## $ drink <int> 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3
## $ well  <int> 5, 5, 6, 6, 3, 5, 4, 6, 8, 6, 5, 6, 6, 6, 6
## $ drunk <int> 5, 3, 2, 1, 7, 6, 6, 4, 2, 3, 2, 3, 2, 3, 2

hc_df$drink <- factor(hc_df$drink, levels = c(1:3), labels = c("Water",
"Lucozade", "Cola"))

Does the type of drink you have the next day predict how well you feel?

oneway_hc <- lm(well ~ drink, data = hc_df)
Anova(oneway_hc, type = "3")

## Anova Table (Type III tests)
## 
## Response: well
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 125.000  1 96.1538 4.425e-07 ***
## drink         2.133  2  0.8205    0.4635    
## Residuals    15.600 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

homogeneity of variance

we want to do a Levene’s test to see whether the variance in well (how well the person feels) varies across the interaction of different types of drinks (drink) and how drunk the person was the night before (drunk). To do this we can execute:

leveneTest(hc_df$well, interaction(hc_df$drink, hc_df$drunk))

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group 10    0.34 0.9242
##        4

assumption: independence of covariate vs predictor

checkIND <- aov(drunk ~ drink, data = hc_df)

Anova(checkIND, type = "3")

## Anova Table (Type III tests)
## 
## Response: drunk
##             Sum Sq Df F value    Pr(>F)    
## (Intercept)   64.8  1 20.9032 0.0006412 ***
## drink          8.4  2  1.3548 0.2947600    
## Residuals     37.2 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

If we control for how drunk you got, does the type of drink you have the next day predict how well you feel?

contrasts(hc_df$drink) <- contr.sum(3)

ancova_hc <- lm(well ~ drink + drunk, data = hc_df)
Anova(ancova_hc, type = "3")

## Anova Table (Type III tests)
## 
## Response: well
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 145.006  1 361.4557 9.197e-10 ***
## drink         3.464  2   4.3177   0.04130 *  
## drunk        11.187  1  27.8860   0.00026 ***
## Residuals     4.413 11                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression output

contrasts(hc_df$drink) <- contr.treatment(3)

ancova_hc2 <- lm(well ~ drink + drunk, data = hc_df)
summary(ancova_hc2)

## 
## Call:
## lm(formula = well ~ drink + drunk, data = hc_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.01935 -0.37742 -0.01935  0.35806  0.99355 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.9742     0.4690  14.869 1.25e-08 ***
## drink2        1.1290     0.4054   2.785  0.01775 *  
## drink3        0.1419     0.4195   0.338  0.74149    
## drunk        -0.5484     0.1038  -5.281  0.00026 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6334 on 11 degrees of freedom
## Multiple R-squared:  0.7512, Adjusted R-squared:  0.6833 
## F-statistic: 11.07 on 3 and 11 DF,  p-value: 0.001193

Testing homogeneity of regression slopes

If the interaction is significant, we need to run a two-way ANOVA (this) and not ANCOVA.

ancova_hc3 <- update(ancova_hc, .~. + drink:drunk)
Anova(ancova_hc3, type = "III")

## Anova Table (Type III tests)
## 
## Response: well
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 59.295  1 194.1961 2.134e-07 ***
## drink        3.376  2   5.5278  0.027166 *  
## drunk        5.216  1  17.0812  0.002548 ** 
## drink:drunk  1.665  2   2.7263  0.118668    
## Residuals    2.748  9                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1