###################################one way random effects model################################################
##install package lme4
library(nlme)
ex.data <- read.table(file="C:/jenn/teaching/stat579/data/influent.txt",header=T)
attach(ex.data)
summary(ex.data)
nrow(ex.data)
head(ex.data)
nrow(ex.data)
head(ex.data)

##descriptive statistics
tapply(nitrogen,influent,mean)
tapply(nitrogen,influent,sd)
tapply(nitrogen,influent,length) #unbalanced design

#declare grouping variable as a factor
influent.f<-as.factor(influent)
boxplot(nitrogen~influent.f,main="Boxplot",xlab = "Influent",
        ylab = "Nitrogen")

# fit one-way ANOVA, assuming only random error
fit1 <- aov(nitrogen~influent.f)
summary(fit1)
plot(fit1$resid~fit1$fitted)
qqnorm(fit1$resid)
qqline(fit1$resid)

# fit linear regression model, see the difference and similarities
fit2<-lm(nitrogen ~ influent.f)
anova(fit2) #compare to summary(fit1)

# fit a random effects model, assuming two nested errors, defaut is REML
fit3 <- lme(nitrogen~1,random=~1|influent.f)
summary(fit3)
intervals(fit3)
fixef(fit3) #overall mean, intercept
ranef(fit3) #alpha_i hat
yhat<-predict(fit3)
resid<-resid(fit3)
plot(resid~yhat)
par(mfrow=c(2,2))
plot(fit3) #standardized residuals v.s fitted values
qqnorm(fit3$resid) #check normality of \epsilon
qqline(fit3$resid)


#refitting model using ML
fit4 <- lme(nitrogen~1,random=~1|influent.f,method="ML") 
summary(fit4)
#results should be similar to fit 3 

##using another package lme4, lmerTest
library(lme4)
library(lmerTest)
fit5<-lmer(nitrogen~(1|influent.f), REML=TRUE,data=ex.data)
summary(fit5)
anova(fit5,df="Kenward-Roger")