#APTS Statistical Modelling 2017
#Practical 2 - Linear Mixed Models

hip<-read.table("hip.txt",col.names=c("y","age","sex","subj","time"))
plot(hip$time,hip$y,type="n",xlab="time",ylab="haemataocrit")
for (i in 1:30) lines(hip$time[hip$subj==i],hip$y[hip$subj==i])


#######################################################################################################
#nlme
library(nlme)
#AIC for model selection (using ML)
AIC(lme(y~factor(time)+age+sex,random= ~1|subj,data=hip,method="ML"))
AIC(lme(y~factor(time)+sex,random= ~1|subj,data=hip,method="ML"))
AIC(lme(y~factor(time),random= ~1|subj,data=hip,method="ML"))
AIC(lme(y~factor(time)+sex,random= ~1+factor(time)|subj,data=hip,method="ML"))
#AIC for linear model excluding subject - comparable to model with random interceopt
AIC(lm(y~factor(time)+sex,data=hip))

#fit model using REML
hip.r<-lme(y~factor(time)+sex,random= ~1|subj,data=hip)
summary(hip.r)

#include subject as fixed effect to produce a "shrinkage" plot
hip.f<-lm(y~factor(subj)+factor(time)+sex,data=hip)
hip.sex<-tapply(hip$sex,hip$subj,mean)#calculate the sex of each subject
#plot the fixed effects subject intercepts against the random effects subject intercepts
plot(hip.f$coef[1]+c(0,hip.f$coef[2:30]),fixef(hip.r)[1]+fixef(hip.r)[4]*hip.sex+ranef(hip.r)[[1]],ylab="intercept (random effect)",xlab="intercept (fixed effect)",pch=ifelse(hip.sex==0,1,2),xlim=c(33,47),ylim=c(33,47))
abline(0,1)
#####################################################################################################

#lme4
library(lme4)
#AIC for model selection (using ML)
AIC(lmer(y~factor(time)+age+sex+(1|subj),data=hip,REML=F))
AIC(lmer(y~factor(time)+sex+(1|subj),data=hip,REML=F))
AIC(lmer(y~factor(time)+(1|subj),data=hip,REML=F))
AIC(lmer(y~factor(time)+sex+(1+factor(time)|subj),data=hip,REML=F))
#AIC for linear model excluding subject - comparable to model with random interceopt
AIC(lm(y~factor(time)+sex,data=hip))

#fit model using REML
hip.lmer<-lmer(y~factor(time)+sex+(1|subj),data=hip)
summary(hip.lmer)

#include subject as fixed effect to produce a "shrinkage" plot
hip.f<-lm(y~factor(subj)+factor(time),data=hip)
hip.sex<-tapply(hip$sex,hip$subj,mean)#calculate the sex of each subject
#plot the fixed effects subject intercepts against the random effects subject intercepts
plot(hip.f$coef[1]+c(0,hip.f$coef[2:30]),coef(hip.lmer)$subj[,1]+hip.sex*coef(hip.lmer)$subj[,4],ylab="intercept (random effect)",xlab="intercept (fixed effect)",pch=ifelse(hip.sex==0,1,2),xlim=c(33,47),ylim=c(33,47))
abline(0,1)


