Purpose

Today’s lab will focus on one-way ANOVA. You can download the rmarkdown file here. We will start with the traditional way like Sara did in class on Thursday. However, the majority of the lab time afterwards will focus on running an ANOVA as a linear regression with a categorical predictor. Though we haven’t talked about it yet, the content will be covered in detail during next week’s class.

Be sure to have the following packages loaded:

library(tidyverse) # for plotting and data wrangling
library(rio) # for importing data
library(psych) # for descriptives
library(emmeans) # for t-tests

Research scenario

We’ll be looking at a dataset comparing different kinds of treatment for depression. In this study, depressed patients (n = 5 per group) were randomly assigned to one of three groups:

CBT group
Psychotherapy group
Control group (received no therapy)

After 10 weeks, participants’ depression scores were measured (on a scale of 1 = no depression to 12 = severely depressed). Our dataset will have just 2 variables: group (CBT, Psychotherapy, or Control) and depress (depression scores).

NOTE: 1 = CBT; 2 = Psychotherapy; 3 = Control

Import the data

depression <- import("https://raw.githubusercontent.com/uopsych/psy612/master/labs/lab-5/data/depression_therapy.csv")

Structure of the data

The first thing we should check is how the data are structured. We should have a factor called group and a numeric variable for depression scores called depress. Let’s use the function str().

str(depression)

## 'data.frame':    15 obs. of  2 variables:
##  $ group  : int  1 1 1 1 1 2 2 2 2 2 ...
##  $ depress: int  9 6 7 6 8 6 3 4 3 1 ...

Recode `group` as a factor

R thinks that group is an integer, so we need to make it into a factor in order to run our analysis. We can use mutate() and base R’s factor(), which can be used to turn variables into factors.

NOTE: You want to provide labels in the same order as the levels of the factor, so in this case we want them in the order CBT, Psychotherapy, and Control (see above).

depression = depression %>% 
  mutate(group = factor(group,
                        labels = c("CBT", "Psychotherapy", "Control"))) # order matters here!

Look at the structure of the data again. Now it’s clear that group is a factor variable with the correct levels.

str(depression)

## 'data.frame':    15 obs. of  2 variables:
##  $ group  : Factor w/ 3 levels "CBT","Psychotherapy",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ depress: int  9 6 7 6 8 6 3 4 3 1 ...

Descriptives

Next, we’ll get descriptives for each group using two tidyverse functions, group_by() and summarize(). First, you want to group by group, and then you want to summarize with three arguments, mean, sd, and n.

Click here to see how we did it in psy611.

depression %>% 
  group_by(group) %>% 
  summarize(mean = mean(depress, na.rm = TRUE),
            sd = sd(depress, na.rm = TRUE),
            n = n()) %>% 
  knitr::kable(digits = 2)

group	mean	sd	n
CBT	7.2	1.30	5
Psychotherapy	3.4	1.82	5
Control	9.8	1.79	5

You can also get descriptives with the describeBy() function. Try adding mat = TRUE as an argument.

describeBy(depression$depress, depression$group, mat = TRUE)

Visualizing the data

There are many ways that you can visualize the data. In ggplot, I recommend adding a geom_boxplot(), or geom_bar(), or geom_violin() layer.

This is a good resource for boxplots with ggplot, this is a good resource for bar plots with ggplot, and this is a good resource for violin plots with ggplot.

depression %>% 
  ggplot(aes(x = fct_reorder(group, depress), y = depress)) +
  geom_boxplot() +
  geom_jitter() +
  theme_minimal() +
  labs(x = NULL,
       y = "Average depression score",
       fill = "Group",
       title = "Depression scores per group")

depression %>% 
  group_by(group) %>% 
  summarise(avg_depress = mean(depress),
            sd = sd(depress)) %>% 
  ggplot(aes(x = fct_reorder(group, avg_depress), 
             y = avg_depress)) +
  geom_bar(stat = "identity",
           fill = "cornflowerblue") +
  geom_point(data = depression, 
             aes(x = group, y = depress),
             position = 'jitter')+
  theme_minimal()+
  labs(y = "Average depression score",
       x = NULL,
       title = "Average depression score per group") +
  geom_errorbar(aes(ymin=avg_depress - sd, ymax=avg_depress + sd),
                width = .1)

depression %>% 
  ggplot(aes(x = fct_reorder(group, depress), y = depress)) +
  geom_violin(aes(fill = group)) +
  geom_jitter() +
  theme_minimal() +
  labs(x = NULL,
       y = "Average depression score",
       fill = "Group",
       title = "Depression scores per group")+
  theme(legend.position = "none")

Traditional ANOVA

We will briefly review how to generate a traditional ANOVA table.

To get an ANOVA table, you can use the same lm() as regression, and use anova() to get the ANOVA summary table. Importantly, you want to make sure that the categorical IV is a factor.

model_anova <- lm(depress ~ group, data = depression)
anova(model_anova)

Question: What does the F test tell us?

To calculate pairwise comparisons between group levels, use emmeans::emmeans(). This function takes an lm() output, an equation that indicates we want to perform pairwise t-tests (pairwise) on the left, and the IV (group) on the right, as well as an adjust (e.g. “bonferroni” or “holm”).

emmeans(model_anova, pairwise ~ group, adjust = "none")

## $emmeans
##  group         emmean    SE df lower.CL upper.CL
##  CBT              7.2 0.739 12     5.59     8.81
##  Psychotherapy    3.4 0.739 12     1.79     5.01
##  Control          9.8 0.739 12     8.19    11.41
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate   SE df t.ratio p.value
##  CBT - Psychotherapy          3.8 1.05 12   3.634  0.0034
##  CBT - Control               -2.6 1.05 12  -2.487  0.0286
##  Psychotherapy - Control     -6.4 1.05 12  -6.121  0.0001

Since we are running a lot of test, we need to correct for multiple comparisons so that we don’t have an inflated type I error rate.

Bonferroni correction:multiplies the p-values by the number of comparisons.

emmeans(model_anova, pairwise ~ group, adjust = "bonferroni")

## $emmeans
##  group         emmean    SE df lower.CL upper.CL
##  CBT              7.2 0.739 12     5.59     8.81
##  Psychotherapy    3.4 0.739 12     1.79     5.01
##  Control          9.8 0.739 12     8.19    11.41
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate   SE df t.ratio p.value
##  CBT - Psychotherapy          3.8 1.05 12   3.634  0.0103
##  CBT - Control               -2.6 1.05 12  -2.487  0.0858
##  Psychotherapy - Control     -6.4 1.05 12  -6.121  0.0002
## 
## P value adjustment: bonferroni method for 3 tests

Holm correction:
- Rank order the p-values from largest to smallest.
- Start with the smallest p-value. Multiply it by its rank. (0.000052 *3 = 0.000156)
- Go to the next smallest p-value. Multiply it by its rank. If the result is larger than the previous step, keep it. Otherwise replace with the previous step adjusted p-value. (0.0034 * 2 = 0.0068)
- Repeat Step 3 for the remaining p-values. (0.0286 * 1 = 0.0286)
- Judge significance of each new p-value against α=.05

emmeans(model_anova, pairwise ~ group, adjust = "holm")

## $emmeans
##  group         emmean    SE df lower.CL upper.CL
##  CBT              7.2 0.739 12     5.59     8.81
##  Psychotherapy    3.4 0.739 12     1.79     5.01
##  Control          9.8 0.739 12     8.19    11.41
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate   SE df t.ratio p.value
##  CBT - Psychotherapy          3.8 1.05 12   3.634  0.0068
##  CBT - Control               -2.6 1.05 12  -2.487  0.0286
##  Psychotherapy - Control     -6.4 1.05 12  -6.121  0.0002
## 
## P value adjustment: holm method for 3 tests

Question: What difference do you notice about the Bonferroni vs. Holm correction?

Regression with categorical predictors

We have just reviewed how to do generate a traditional ANOVA table. We will next talk about how to perform ANOVA within a linear regression framework. This will give us the same information we already got from our regression model, and something more!

Dummy coding

In regression, categorical predictors with more than two levels are broken up into more than one predictor. This actually happens under the hood in R when our independent variable is a factor. Let’s put group as the predictor in our model as before:

model_default <- lm(depress ~ group, data = depression) 
summary(model_default)

## 
## Call:
## lm(formula = depress ~ group, data = depression)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.4   -1.2   -0.2    1.0    2.6 
## 
## Coefficients:
##                    Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)          7.2000     0.7394   9.738 0.000000477 ***
## groupPsychotherapy  -3.8000     1.0456  -3.634     0.00342 ** 
## groupControl         2.6000     1.0456   2.487     0.02861 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.653 on 12 degrees of freedom
## Multiple R-squared:  0.7595, Adjusted R-squared:  0.7195 
## F-statistic: 18.95 on 2 and 12 DF,  p-value: 0.0001934

Here is our regression equation:

\[Depression_i = \beta_0 + \beta_1psychotherapy_i + \beta_2control_i + e_i\]

R is creating two dummy variables for us under the hood: groupPsychotherapy and groupControl. This is what R is doing:

depression %>%
  mutate(groupPsychotherapy = ifelse(group == "Psychotherapy", 1, 0),
         groupControl = ifelse(group == "Control", 1, 0))

By default, R treats whatever the first level of the factor variable is as the reference group. In this case, CBT was the first level of group (because it was initially coded as 1 in our raw data), so the model treated CBT as our reference group.

Here is our regression equation, derived from the output:

\[Depression_i = 7.20 -3.80psychotherapy + 2.60control\]

Question: Determine the mean depression score of each group.

Now, let’s revisit the output.

model_default <- lm(depress ~ group, data = depression) 
summary(model_default)

## 
## Call:
## lm(formula = depress ~ group, data = depression)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.4   -1.2   -0.2    1.0    2.6 
## 
## Coefficients:
##                    Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)          7.2000     0.7394   9.738 0.000000477 ***
## groupPsychotherapy  -3.8000     1.0456  -3.634     0.00342 ** 
## groupControl         2.6000     1.0456   2.487     0.02861 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.653 on 12 degrees of freedom
## Multiple R-squared:  0.7595, Adjusted R-squared:  0.7195 
## F-statistic: 18.95 on 2 and 12 DF,  p-value: 0.0001934

Question: Interpret the intercept. What does a significant p-value signify?

Question: Interpret the slopes. What do the significant p-values signify?

Let’s imagine that we have the following (more intuitive) research questions:

Is CBT effective (relative to no therapy)?
Is psychotherapy effective (relative to no therapy)?

Question: What do want our reference group to be to answer these research questions?

Now let’s make the appropriate dummy codes. Recall that we need k-1 dummy codes (where k is the number of groups). We have 3 groups, so we need 2 dummy codes.

Remember, our decision of how to set the dummy codes (which group to set as the reference group) should be guided by our research questions.

Method 1: Create dummy variables with `mutate()`

# create new dummy variables containing 1's and 0's
depression <- depression %>%
  mutate(CBT.v.Control = 
           ifelse(group == "CBT", 1, 0), 
         Psychotherapy.v.Control = 
           ifelse(group == "Psychotherapy", 1, 0)) 

depression

Now, we can run the linear model using these new dummy variables as the IV’s.

model_dummy1 <- lm(depress ~ CBT.v.Control + Psychotherapy.v.Control, data = depression) 
summary(model_dummy1)

## 
## Call:
## lm(formula = depress ~ CBT.v.Control + Psychotherapy.v.Control, 
##     data = depression)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.4   -1.2   -0.2    1.0    2.6 
## 
## Coefficients:
##                         Estimate Std. Error t value     Pr(>|t|)    
## (Intercept)               9.8000     0.7394  13.255 0.0000000159 ***
## CBT.v.Control            -2.6000     1.0456  -2.487       0.0286 *  
## Psychotherapy.v.Control  -6.4000     1.0456  -6.121 0.0000517062 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.653 on 12 degrees of freedom
## Multiple R-squared:  0.7595, Adjusted R-squared:  0.7195 
## F-statistic: 18.95 on 2 and 12 DF,  p-value: 0.0001934

Question: What does the intercept mean? What do the slopes mean?

Method 2: Contrast Matrix

Now we’re going use another method to create our dummy codes by creating a contrast matrix. When we pass a factor to lm() it uses a contrast matrix under the hood. Although R creates the contrast matrix automatically with certain default assumptions, you can set it yourself if you have a specific set of codes in mind (which we do).

Use contrasts() to view the default contrast matrix associated with our grouping variable.

contrasts(depression$group)

##               Psychotherapy Control
## CBT                       0       0
## Psychotherapy             1       0
## Control                   0       1

Again, by default our contrasts are not how we want them because we want Control to be the reference group. To manually specify our dummy codes, we’ll use the function contr.treatment(), which will return a contrast matrix with specific characteristics, depending on what we tell it. The first argument, n is the number of levels of the factor variable. The second argument, base refers to the level of the factor that should be considered the reference group.

contr.treatment(n = 3, base = 3)

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

We can now overwrite the default contrast matrix for group with the one we just manually specified.

contrasts(depression$group) <- contr.treatment(n = 3, base = 3)

Now look at the structure of the data again. What do you notice about group?

str(depression)

## 'data.frame':    15 obs. of  4 variables:
##  $ group                  : Factor w/ 3 levels "CBT","Psychotherapy",..: 1 1 1 1 1 2 2 2 2 2 ...
##   ..- attr(*, "contrasts")= num [1:3, 1:2] 1 0 0 0 1 0
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:3] "CBT" "Psychotherapy" "Control"
##   .. .. ..$ : chr [1:2] "1" "2"
##  $ depress                : int  9 6 7 6 8 6 3 4 3 1 ...
##  $ CBT.v.Control          : num  1 1 1 1 1 0 0 0 0 0 ...
##  $ Psychotherapy.v.Control: num  0 0 0 0 0 1 1 1 1 1 ...

Let’s take a closer look at group using attributes()…

attributes(depression$group)

## $levels
## [1] "CBT"           "Psychotherapy" "Control"      
## 
## $class
## [1] "factor"
## 
## $contrasts
##               1 2
## CBT           1 0
## Psychotherapy 0 1
## Control       0 0

Now, run the linear model with group (now containing our desired dummy codes) as the IV.

model_dummy2 <- lm(depress ~ group, data = depression) 
summary(model_dummy2)

## 
## Call:
## lm(formula = depress ~ group, data = depression)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.4   -1.2   -0.2    1.0    2.6 
## 
## Coefficients:
##             Estimate Std. Error t value     Pr(>|t|)    
## (Intercept)   9.8000     0.7394  13.255 0.0000000159 ***
## group1       -2.6000     1.0456  -2.487       0.0286 *  
## group2       -6.4000     1.0456  -6.121 0.0000517062 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.653 on 12 degrees of freedom
## Multiple R-squared:  0.7595, Adjusted R-squared:  0.7195 
## F-statistic: 18.95 on 2 and 12 DF,  p-value: 0.0001934

Method 3: Re-order levels of your factor

Let’s take another look at the levels of group

levels(depression$group)

## [1] "CBT"           "Psychotherapy" "Control"

We can manually change the order of the levels by using factor() and manually specifying the levels argument to contain the order we want. Again, since R treats the first level of the factor as the reference group, and we want Control to be the reference group in order to obtain our contrasts of interests when we run the linear model, we will create a new version of group called group_relevel and make Control the first level of the factor in this new variable.

depression <- depression %>% 
  mutate(group_relevel = factor(group, levels = c("Control",
                                                  "CBT",
                                                  "Psychotherapy")))

To double check that this worked, let’s look at the levels of group_relevel

levels(depression$group_relevel)

## [1] "Control"       "CBT"           "Psychotherapy"

Now when we run the linear model, the contrast matrix that will be automatically generated will correspond to the correct contrasts we want. To verify this, let’s re-run the linear model (yet again!) with the re-leveled group_relevel variable as the IV.

model_dummy3 <- lm(depress ~ group_relevel, data = depression) 
summary(model_dummy3)

## 
## Call:
## lm(formula = depress ~ group_relevel, data = depression)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.4   -1.2   -0.2    1.0    2.6 
## 
## Coefficients:
##                            Estimate Std. Error t value     Pr(>|t|)    
## (Intercept)                  9.8000     0.7394  13.255 0.0000000159 ***
## group_relevelCBT            -2.6000     1.0456  -2.487       0.0286 *  
## group_relevelPsychotherapy  -6.4000     1.0456  -6.121 0.0000517062 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.653 on 12 degrees of freedom
## Multiple R-squared:  0.7595, Adjusted R-squared:  0.7195 
## F-statistic: 18.95 on 2 and 12 DF,  p-value: 0.0001934

Other coding schemes

There are many other coding schemes besides dummy coding, such as and deviation coding (also known as “effects coding”), which can compare individual group means to the overall mean of the dependent variable (i.e. the grand mean). You can read more about that here.

Minihacks

Research Problem:

You are interested in the effectiveness of a new drug in treating the symptoms of depression. You randomly assign patients with depression into one of 3 treatment groups with different levels of the drug (low, medium, and high), and also one control group that only received a placebo. You then measure their symptoms after 4 weeks of drug treatment. You are interested in whether the drug is effective, and what level of the drug is most effective.

Load in the raw data with the following code

depression_drug <- import("https://raw.githubusercontent.com/uopsych/psy612/master/labs/lab-5/data/depression_drug.csv")

Minihack 1

Recreate the following table.

Recreate the following boxplot.

What pattern do you notice in the data?

Minihack 2

Run an omnibus ANOVA test to determine whether there is a relationship between drug and depression. What does this tell you?
Using your preferred method of dummy coding, test whether there is a significant difference between depression scores for each drug group (low, medium, high) and the none group. Interpret all model coefficients.

Minihack 3

Conduct pairwise comparisons for depression scores across all levels of drug. Use a correction method to account for multiple comparisons.

Lab 5: Categorical Predictors/One-way ANOVA