# R code for Lazic SE (2008). Why we should use simpler models
# if the data allow this: relevance for ANOVA designs
# in experimental biology. BMC Physiology.


#create vectors for dose and immobility times
dose <- c(rep(0,5), rep(80, 5), rep(160,5), rep(240,5))
time.immob <- c(182, 112, 206, 170, 164, 158, 165, 168, 182,
	97, 140, 135, 110, 128, 155, 163, 183,  25, 100, 61)


mod1 <- lm(time.immob~factor(dose)) #"ANOVA" analysis
summary.aov(mod1)		#provides ANOVA table
par(mfrow=c(2,2)); plot(mod1) #diagnostic plots


mod2 <- lm(time.immob~dose) #"Regression" analysis
summary.aov(mod2)
par(mfrow=c(2,2)); plot(mod2)


#testing for lack of fit in the regression model
anova(mod2,mod1) #the ANOVA model is not a significantly better fit over the regression model (p = 0.929)


#testing for lack of fit "by hand", see equation 3 in the paper
(28302.6-28043.2)/(18-16) / (28043.2/16)  #same as above
pf(0.07400011, 2,16, lower.tail=F)        #same as above


#model comparison with Akaike's Information Criterion
AIC(mod1,mod2) # the regression model has the lower AIC and is therefore preferred


#model comparison with the Bayesian Information Criterion
-2*logLik(mod1) + 5*log(20) # ANOVA model
-2*logLik(mod2) + 3*log(20) # regression (preferred)


#multiple comparisons
pairwise.t.test(time.immob,factor(dose), p.adjust.method="none")  # 0 vs. 240 is sig (p = 0.037) with no correction
pairwise.t.test(time.immob,factor(dose), p.adjust.method="bonf")  # 0 vs. 240 is not sig (p = 0.22) with Bonferroni correction  

fac.dose <- factor(dose) #need to make dose into a new variable (factor)
mod1.multi <- aov(time.immob~fac.dose)

library(multcomp) #need to load "multcomp" package
multi.mod <- glht(mod1.multi, linfct = mcp(fac.dose="Dunnett"))
summary(multi.mod) # 0 vs. 240 is not sig (p = 0.090) with Dunnett's correction