A psychologist was interested whether certain TV shows result in a more positive attitude towards life. In a study, the participants were randomly allocated to four rooms to view one of the following: The Muppet Show, Game of Thrones, ONENEWS, no show (control). After the experiment a blood sample was taken and serotonin levels (happy hormone) were measured (the more serotonin the happier!).  

  1. Perform an exploratory data analysis (draw a boxplot).  

  2. Run an analysis of variance (ANOVA) model using the lm function.  

  3. Check the model diagostic plots for violoations of the variance homogeneity and normality assumptions.  

  4. Perform a post-hoc test to determine which groups differ significantly from each other (use the emmeans package and/or the multcomp package).  

  5. Create a publication-quality figure containing statistical annotation (e.g. lower case letters indicating significant differences between groups).  

show <- gl(n = 4, k = 7, labels = c("Control", "The Muppet Show", "Game of Thrones", "ONENEWS"))
serotonin <- c(7, 7, 5, 4, 3, 4, 4, 11, 7, 8, 14, 11, 10, 5, 4, 8, 6, 11, 9, 8, 7, 4, 3, 2, 2, 3, 6, 3)
tv.att <- data.frame(serotonin, show)

## Graphical data exploration
## Base graphics
boxplot(serotonin ~ show, data = tv.att)

## ggplot
library(ggplot2)

ggplot(data = tv.att, mapping = aes(x = show, y = serotonin)) + geom_boxplot()

## ANOVA model
m1 <- lm(serotonin ~ show, data = tv.att)
summary.aov(m1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## show         3  158.0   52.67   11.52 7.13e-05 ***
## Residuals   24  109.7    4.57                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m1)
## 
## Call:
## lm(formula = serotonin ~ show, data = tv.att)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4286 -1.2857 -0.2857  1.4643  4.5714 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           4.8571     0.8081   6.010 3.32e-06 ***
## showThe Muppet Show   4.5714     1.1429   4.000 0.000527 ***
## showGame of Thrones   2.7143     1.1429   2.375 0.025886 *  
## showONENEWS          -1.5714     1.1429  -1.375 0.181833    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.138 on 24 degrees of freedom
## Multiple R-squared:  0.5902, Adjusted R-squared:  0.539 
## F-statistic: 11.52 on 3 and 24 DF,  p-value: 7.129e-05
## Post-hoc comparison
library(emmeans)
emms <- emmeans(object = m1, specs = "show") # calculate least-squares estimates of the means
pairs(emms)
##  contrast                           estimate       SE df t.ratio p.value
##  Control - The Muppet Show         -4.571429 1.142857 24  -4.000  0.0028
##  Control - Game of Thrones         -2.714286 1.142857 24  -2.375  0.1093
##  Control - ONENEWS                  1.571429 1.142857 24   1.375  0.5264
##  The Muppet Show - Game of Thrones  1.857143 1.142857 24   1.625  0.3842
##  The Muppet Show - ONENEWS          6.142857 1.142857 24   5.375  0.0001
##  Game of Thrones - ONENEWS          4.285714 1.142857 24   3.750  0.0051
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
library(multcomp)
# summary(as.glht(pairs(lsm))) # single-step P-value adjustment
summary(as.glht(pairs(emms)), test = adjusted(type = "BH")) # Benjamini & Hochberg (1995) adjustment
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Linear Hypotheses:
##                                        Estimate Std. Error t value
## Control - The Muppet Show == 0           -4.571      1.143  -4.000
## Control - Game of Thrones == 0           -2.714      1.143  -2.375
## Control - ONENEWS == 0                    1.571      1.143   1.375
## The Muppet Show - Game of Thrones == 0    1.857      1.143   1.625
## The Muppet Show - ONENEWS == 0            6.143      1.143   5.375
## Game of Thrones - ONENEWS == 0            4.286      1.143   3.750
##                                        Pr(>|t|)    
## Control - The Muppet Show == 0          0.00158 ** 
## Control - Game of Thrones == 0          0.03883 *  
## Control - ONENEWS == 0                  0.18183    
## The Muppet Show - Game of Thrones == 0  0.14067    
## The Muppet Show - ONENEWS == 0         9.68e-05 ***
## Game of Thrones - ONENEWS == 0          0.00198 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- BH method)
## Publication-quality figure
bp <- boxplot(serotonin ~ show, data = tv.att, col = "grey", ylim = c(0, 16), las = 1, pars = list(boxwex = 0.6, whisklty = 1, staplewex = 0.3))
mtext("Serotonin level", side = 2, line = 2.5) # side = 2 specifies the y-axis, line specifies the distance of the axis name to the actual axis.
text(x = 1:4, y = bp$stats[5, ] + 1, labels = c("a", "b", "b", "a"), cex = 1.2) # cex = character expansion (font size)