#LMM 05 CONTROLLING VARIANCE WITH GLS
library(nlme)
setwd("c:/DATA/Models")
#READING DATA AND DEFINING MONTH AS A FACTOR
SQ=read.table("Squid.txt",header=T)
SQ$fMONTH=factor(SQ$MONTH) #PUTS FACTOR INSIDE SQ
SQ

M <- lm(Testisweight ~ DML * fMONTH,data=SQ)
anova(M)
summary(M)

M.LM=gls(Testisweight~DML*fMONTH,data=SQ)
anova(M.LM)
summary(M.LM)

#VALIDATION PLOTS
op=par(mfrow = c(2,2), mar=c(4,4,2,2))
plot(M)
par(op)

#HOMEMADE PLOTS
op=par(mfrow = c(2,2), mar=c(4,4,2,2))
plot(fitted(M),resid(M),col='blue',xlab='Fitted Values',ylab='Residuals', pch=20)
plot(SQ$fMONTH, resid(M.LM),pch=19,col='green',xlab="Month",ylab="Residuals")
plot(SQ$DML, resid(M.LM),col='red',xlab="DML",ylab="Residuals",pch=20)
hist(resid(M.LM),col="brown",main='',xlab="Residuals")
par(op)

M.varFixed=gls(Testisweight~DML*fMONTH,
              weights=varFixed(~DML),
			  data=SQ)
			  
M.varIdent=gls(Testisweight ~ DML*fMONTH,
              weights=varIdent(form= ~ 1|fMONTH), 
			  data =SQ)

M.varPowerA=gls(Testisweight ~ DML * fMONTH,
              weights = varPower(form =~ DML),
			  data=SQ)

M.varPowerB=gls(Testisweight ~ DML * fMONTH,
              weights = varPower(form=~ DML | fMONTH),
			  data = SQ)


M.varExp=gls(Testisweight ~ DML * fMONTH,
              weights = varExp(form =~ DML), 
			  data = SQ)

M.varConstPowerA=gls(Testisweight ~ DML * fMONTH,
            weights = varConstPower(form =~ DML), 
			data = SQ)


M.varConstPowerB=gls(Testisweight ~ DML * fMONTH,
              weights = varConstPower(form =~ DML | fMONTH), 
			  data = SQ)

M.varComb=gls(Testisweight ~ DML * fMONTH,
              weights = varComb(varIdent(form= ~ 1 | fMONTH),
						varExp(form =~ DML) ), 
			  data = SQ)

anova(	M.LM,
		M.varFixed,
		M.varIdent,
		M.varPowerA,
		M.varPowerB,
		M.varExp,
		M.varConstPowerA,
		M.varConstPowerB,
		M.varComb)

AIC(	M.LM,
		M.varFixed,
		M.varIdent,
		M.varPowerA,
		M.varPowerB,
		M.varExp,
		M.varConstPowerA,
		M.varConstPowerB,
		M.varComb)
		
#TESTING IMPROVEMENT OF M.varPowerB MODEL
anova(M.LM,M.varPowerB)

#VALIDATION PLOTS
M2=M.varPowerB
op=par(mfrow = c(2,2), mar=c(4,4,2,2))
plot(fitted(M2),resid(M2),col='blue',xlab='Fitted Values',ylab='Residuals', pch=20)
plot(SQ$fMONTH, resid(M2),col='green',xlab="Month",ylab="Residuals")
plot(SQ$DML, resid(M2),col='red',xlab="DML",ylab="Residuals",pch=20)
hist(resid(M2),col="brown",main='',xlab="Residuals")
par(op)

#PLOTTING STANDARDIZED RESIDUALS
plot(M2,pch=20,col='blue')