########################################################
### Example: One-way and two-way ANOVA with B&O data ###
########################################################

# Get the B&O data from the lmerTest-package
library(lmerTest)
data(TVbo)
head(TVbo) # First rows of the data

# Define factor identifying the 12 TV set and picture combinations 
TVbo$TVPic <- factor(TVbo$TVset:TVbo$Picture)
 
# Each of 8 assessors scored each of the 12 combinations twice.
# Average the two replicates for each assessor and combination of 
# TV set and picture
library(doBy)
TVbonoise <- summaryBy(Noise ~ Assessor + TVPic, data = TVbo, 
                       keep.names = T)

# One-way ANOVA of the noise (not the correct analysis!)
anova(lm(Noise ~ TVPic, data = TVbonoise))

# Two-way ANOVA of the noise (better analysis, week 12)
anova(lm(Noise ~ Assessor + TVPic, data = TVbonoise))

##############################################
### Simple example: Plots of data by group ###
##############################################

# Input data
y <- c(2.8, 3.6, 3.4, 2.3,
       5.5, 6.3, 6.1, 5.7,
       5.8, 8.3, 6.9, 6.1)

## Define treatment groups
treatm <- factor(c(1, 1, 1, 1,
                   2, 2, 2, 2,
                   3, 3, 3, 3))

## Plot data by treatment groups
par(mfrow = c(1,2))
plot(y ~ as.numeric(treatm), xlab = "Treatment", ylab = "y")
boxplot(y ~ treatm, xlab = "Treatment", ylab = "y")

###############################################
###  Plot F-distribution and critical value ###
###############################################

# Remember, this is "under H0" (i.e. we compute as if H0 is true)

# Number of groups
k <- 3

# Total number of observations
n <- 12

# Sequence for plot
xseq <- seq(0, 10, by = 0.1)

# Plot density of the F-distribution
plot(xseq, df(xseq, df1 = k-1, df2 = n-k), type = "l")

# Plot critical value for significance level 5%
cr <- qf(0.95, df1 = k-1, df2 = n-k)
abline(v = cr, col = "red") 

############################################
### One-way ANOVA using anova() and lm() ###
############################################

anova(lm(y ~ treatm))

###############################
### One-way ANOVA 'by hand' ###
###############################

k <- 3; n <- 12  # Number of groups k, total number of observations n

# Total variation, SST
(SST <- sum( (y - mean(y))^2 ))

# Residual variance after model fit, SSE
y1 <- y[1:4]; y2 <- y[5:8]; y3 <- y[9:12]

(SSE <- sum( (y1 - mean(y1))^2 ) + 
        sum( (y2 - mean(y2))^2 ) + 
        sum( (y3 - mean(y3))^2 ))

# Variance explained by the model, SS(Tr)
(SSTr <- SST - SSE)

# Test statistic
(Fobs <- (SSTr/(k-1)) / (SSE/(n-k)))

# P-value
(1 - pf(Fobs, df1 = k-1, df2 = n-k))

########################
### Model validation ###
########################

# Check assumption of homogeneous variance using, e.g., 
# a box plot.
plot(treatm, y)

# Check normality of residuals using a normal QQ-plot
fit1 <- lm(y ~ treatm)
qqnorm(fit1$residuals)
qqline(fit1$residuals)