Load packages
library(mosaic)
library(dplyr)
library(ggplot2)
library(ggvis)
library(parallel)


AxS ANOVA: Optimism Example

## Load data
optimism <- read.csv("http://www.bradthiessen.com/html5/data/optimism.csv")

## Look at variable types
str(optimism)
## 'data.frame':    54 obs. of  3 variables:
##  $ subject : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ age     : int  16 16 16 16 16 16 16 16 16 16 ...
##  $ optimism: int  35 32 33 32 31 29 29 27 27 28 ...
## Change subject and age to factor variables
optimism$age <-as.factor(optimism$age)
optimism$subject <-as.factor(optimism$subject)

## Peek at data
head(optimism)
##   subject age optimism
## 1       1  16       35
## 2       2  16       32
## 3       3  16       33
## 4       4  16       32
## 5       5  16       31
## 6       6  16       29
## Summary of data
favstats(optimism ~ age, data=optimism)
##   .group min    Q1 median    Q3 max     mean       sd  n missing
## 1     16  13 24.00     27 30.50  35 25.88889 6.578833 18       0
## 2     18  13 25.25     29 31.75  39 26.94444 7.487026 18       0
## 3     20  12 19.00     24 27.75  32 23.38889 5.922429 18       0
## Some summary plots
densityplot(~optimism|age, data=optimism, lwd=2,
    main="Optimism by age",
   layout=c(1,3))

## Let's also run a test to check for homogeneity of variances
bartlett.test(optimism~age, data=optimism)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  optimism by age
## Bartlett's K-squared = 0.9167, df = 2, p-value = 0.6323
## Regular one-way ANOVA
aov.out <- lm(optimism ~ age, data=optimism)
anova(aov.out)
## Analysis of Variance Table
## 
## Response: optimism
##           Df  Sum Sq Mean Sq F value Pr(>F)
## age        2  120.04  60.019  1.3396  0.271
## Residuals 51 2285.00  44.804
## No significant differences among age means
## Check assumptions
plot(aov.out)

##############################################

# Subject means
mean(optimism ~ subject, data=optimism)
##        1        2        3        4        5        6        7        8 
## 35.33333 32.66667 31.00000 31.00000 30.00000 29.33333 29.00000 27.66667 
##        9       10       11       12       13       14       15       16 
## 27.33333 26.33333 25.66667 25.33333 24.00000 22.66667 17.00000 16.33333 
##       17       18 
## 13.66667 13.00000
sd(optimism$optimism)
## [1] 6.736324
## AxS ANOVA
# Take a look at profile plots (optimism by subject at each age)
interaction.plot(x.factor=optimism$age, trace.factor=optimism$subject,
                 response=optimism$optimism,
                 fun=mean, type="b", legend=F)