##### Mixed-effects models


###import data from web


datum=read.csv(file.choose())
head(datum)


### analyze the data without including a random effect for field


results=lm(Biomass~Treatment,data=datum)
summary(results)
#Notice t-value for effect of treatment


### Now run model with random effect for field


library(nlme)
results2=lme(Biomass~Treatment,data=datum,random=(~1|Field))
summary(results2)
#Note that t-value for effect of treatment is stronger
# also note that p-value for effect of treatment is smaller, despite lower residual df
# effect of field was 'swamping' ability to detect treatment effect


### Note that this linear model is the same as a 'paired t-test'


t.test(Biomass~Treatment,paired=TRUE,data=datum,var.equal=TRUE)


### Above model only assumes that overall biomass may vary from field to field (random intercept)
### doesn't assume that effect of treatment may vary from field to field
### which technically would be a field*Treatment interaction (random slope)


###import datum2 from web for looking at random slopes


datum2=read.csv(file.choose())
head(datum)


### Model with random intercept for field


results2=lme(Biomass~Treatment,data=datum2,random=(~1|Field))
summary(results2)
#no interaction between field and treatment in model


### Model with random intercept for field and random slope - field influences effect of treatment


results3=lme(Biomass~Treatment,data=datum2,random=(~Treatment|Field))
summary(results3)
#models interaction between treatment and field
#note that results do not affect significance of treatment


### Rarely include random slopes because don't effect significance of fixed effects
### Also complicate models


###Test significance of Treatment effect using F-drop Test


results4=lme(Biomass~1,data=datum2,random=(~1|Field))
anova(results4,results3)


###Notice warning message
###lme are fit (default) using REML, which is best for comparing models with 
###different random effects or for estimating coefficients of model (summary)
###To compare models with different fixed effects, models must be fit using
###maximum likelihood


results2a=lme(Biomass~Treatment,data=datum2,random=(~1|Field),method="ML")
results4=lme(Biomass~1,data=datum2,random=(~1|Field),method="ML")
anova(results2a,results4) # no warning message
# model with Treatment is better


### be sure to use REML when fit model to get coefficient estimates
summary(results2) ###NOT results2a