################################################################
#
#				Mixed effects linear model: the Loire River example
#
################################################################

require(ade4)
require(entropart)
require(nlme)
#################### recall the dataset
data(macroloire)
######### isolate the faunistic table and the environmental variables
dat.1=macroloire$fau
env=macroloire$envir
############
####################### build the metacommunity and draw the entropy profile, we choose q=1 to start the example
Meta.dat.st=MetaCommunity(dat.1,Weights=apply(dat.1,2,sum))
qq=seq(0,2,by=0.1)
DivP.dat=DivProfile(qq,Meta.dat.st)
plot(DivP.dat)
#
################# Compute Alpha diversity with communities given by the monitoring stations
#
Alpha.dat.1=AlphaEntropy(q=1,MC=Meta.dat.st,Correction="ChaoShen")$Communities

#### attach the environmental table so to simplify the commands writing
attach(env)
############################
#
#		Some graphical exploration of the data to verify wether or not random effects seems suitable
#
#####################
boxplot(Alpha.dat.1~Morphoregion,col=rainbow(3))
plot(Distance,Alpha.dat.1,pch=20,xlab="distance from the source", ylab=expression(paste(alpha,"-diversity")))
abline(lm(Alpha.dat.1~Distance))
points(Distance,Alpha.dat.1,pch=20,col=c(2,3,4)[Morphoregion],cex=0.5)
legend("topright",c("regression line",levels(Morphoregion)),lty=c(1),lwd=c(1,1,1,1),col=c(1,2,3,4),cex=0.6)

plot(Altitude,Alpha.dat.1,pch=20,xlab="Altitude", ylab=expression(paste(alpha,"-diversity")))
abline(lm(Alpha.dat.1~Altitude))
points(Altitude,Alpha.dat.1,pch=20,col=c(2,3,4)[Morphoregion],cex=0.5)
legend("bottomright",c("regression line",levels(Morphoregion)),lty=c(1),lwd=c(1,1,1,1),col=c(1,2,3,4),cex=0.6)
###################################
#
#		Linear models no random effects
#
#######################################
summary(lm(Alpha.dat.1~Distance))
summary(lm(Alpha.dat.1~Altitude))
#######################################
#
#		Compare model with Distance and with Altitude
#
#######################################
mod.1d=lme(Alpha.dat.1~Distance,random=~1|Morphoregion,method="ML")
mod.1a=lme(Alpha.dat.1~Altitude,random=~1|Morphoregion,method="ML")

summary(mod.1d)
summary(mod.1a)

barplot(t(mod.1a$coefficients$random$Morphoregion),beside=T,col=rainbow(3),ylim=c(-0.25,0.3))
abline(h=0)

######################################################
#
#		The AIC of the model with Altitude is smaller than the one for the model with Distance from the source
#######################################################
######################################################
#
#		we consider the model with random intercept as well
#
#######################################################
mod.1as=lme(Alpha.dat.1~Altitude,random=~1+Altitude|Morphoregion,method="ML")
summary(mod.1as)
#This model has a larger AIC value than the one with the random intercept alone. 
#		Graphics 
#
F0=fitted(mod.1a,level=0)
F1=fitted(mod.1a,level=1)
I=order(Altitude)
DD=sort(Altitude)
plot(DD,F0[I],lwd=3,type="l",ylim=c(0,2),xlab="Altitude",ylab=expression(paste(alpha,"-diversity")))
aa=levels(Morphoregion)
colo=c(2,3,4)
for(i in 1:length(aa)){
	x1=Altitude[Morphoregion==aa[i]]
	y1=F1[Morphoregion==aa[i]]
	K=order(x1)
lines(sort(x1),y1[K],col=colo[i],lwd=4,lty=5)
abline(lm(y1[K]~sort(x1)),col=1,lwd=1)
lines(sort(x1),y1[K],col=colo[i],lwd=4)

}

points(Altitude,Alpha.dat.1,pch=25,col=c(2,3,4)[Morphoregion],cex=0.5)
legend("bottomright",c("no random effect",aa),lwd=c(3,1,1,1),col=c(1,2,3,4))


########## model with both random intercept and random slope

F0=fitted(mod.1as,level=0)
F1=fitted(mod.1as,level=1)
I=order(Altitude)
DD=sort(Altitude)
plot(DD,F0[I],type="n",ylim=c(0,2),xlab="Altitude",ylab=expression(paste(alpha,"-diversity of order 1")))
abline(lm(F0[I]~DD),lwd=4)
aa=levels(Morphoregion)
colo=c(2,3,4)
for(i in 1:length(aa)){
	x1=Altitude[Morphoregion==aa[i]]
	y1=F1[Morphoregion==aa[i]]
	K=order(x1)
#lines(sort(x1),y1[K],col=colo[i],lwd=4,lty=5)
#lines(sort(x1),y1[K],col=colo[i],lwd=4)
abline(lm(y1[K]~sort(x1)),col=1,lwd=1)
lines(sort(x1),y1[K],col=colo[i],lwd=4)

}

points(Altitude,Alpha.dat.1,pch=25,col=c(2,3,4)[Morphoregion],cex=0.5)
legend("bottomright",c("no random effect",aa),lwd=c(3,1,1,1),col=c(1,2,3,4))