setwd("C:/Users/user/R Scripts/Linear Models/Data Files")

par(pty = "s", oma = c(0., 0., 4., 0.))
library(nlme)
	react.dat <- read.table("React_time.dat", head = T)
	attach(react.dat)
	print(react.dat)

#-----------------------------------------------------------------------------
#            FIXED EFFECTS ANOVA
#-----------------------------------------------------------------------------
	par(mfrow = c(2., 2.))
	react1.ran <- aov(log(Time) ~ factor(Drug) + factor(Subject))
	plot(Subject, log(Time))
	points(Subject, react1.ran$fitted, pch = "+")
	plot(Drug, log(Time))
	points(Drug, react1.ran$fitted, pch = "+")
	print(summary(react1.ran))
#	plot(react1.ran$fitted, react1.ran$res)
	plot(react1.ran,which=1:2)

#-----------------------------------------------------------------------------
#           MIXED EFFECTS ANOVA
#-----------------------------------------------------------------------------

#	react2.ran <- aov(log(Time) ~ factor(Drug) + Error(factor(Subject)))
	react2.ran <- lme(log(Time) ~ factor(Drug), random=~1|Subject)
	print(summary(react2.ran,correlation=F))
	print(anova(react2.ran))
	par(mfrow = c(2., 2.))
	plot(Drug,log(Time))
	points(Drug,react2.ran$fitted[,1],pch="+")
      plot(fitted(react2.ran), resid(react2.ran))
      qqnorm(react2.ran$res)
      qqnorm(ranef(react2.ran)[,1])