library(nlme)
#######
#Wheat Varieties
############
D=read.table(file="Wheat.dat",header=T)
GD=groupedData(Yield~Moist|Variety,data=D)
plot(GD)
lmfit1=lmList(Yield~Moist|Variety,data=GD)
plot(intervals(lmfit1))
lmefit1=lme(Yield~Moist,data=GD,random=~Moist|Variety)
summary(lmefit1)
ranef(lmefit1) 
coef(lmefit1)
lmefit2=lme(Yield~Moist,data=GD,random=list(Variety=pdDiag(form=~Moist)))
anova(lmefit1,lmefit2)
plot(augPred(lmefit2),aspect="xy",grid=T)


###############
# Pig Weights
#####
D=read.table(file="Pigs.dat",header=T)
GD=groupedData(Weight~Day|Pig,data=D,outer=~factor(Treat))
plot(GD,outer=T)
lmefit1=lme(Weight~factor(Treat)*Day+factor(Treat)*I(Day^2),data=GD,random=list(Pig=pdDiag(form=~Day+I(Day^2))))
lmefit2=lme(Weight~factor(Treat)*Day+factor(Treat)*I(Day^2),data=GD,random=list(Pig=pdDiag(form=~Day+I(Day^2))), corr=corARMA(form=~1|Pig,p=1,q=0))
anova(lmefit1,lmefit2)
plot(augPred(lmefit1),aspect="xy",grid=T)
lmfit1=lmList(Weight~Day+I(Day^2),data=GD)
plot(comparePred(lmefit1,lmfit1))



########################
# Relating Discharge to Temprature and Precipitation for 15 randomly selected streams in MT observed for over 50 years
##############
D=read.csv(file="Discharge.csv",header=T)
D$Discharge=log(D$Discharge)
library(nlme)
GD1=groupedData(Discharge~Temp|SN,data=D)
plot(GD1)
GD2=groupedData(Discharge~Precip|SN,data=D)
plot(GD2)
lmListfit=lmList(Discharge~ Temp + Precip + WS,data=GD1)
plot(intervals(lmListfit))
lmefit1=lme(Discharge ~ Temp + Precip + WS,data=GD1, random=list(SN=pdDiag(~1+Precip + Temp)))
#plot(augPred(lmefit1),grid=T)
cf=compareFits(coef(lmefit1),coef(lmListfit))
plot(cf,mark=fixef(lmefit1))
#plot(comparePred(lmefit1,lmListfit,length.out=2))
#plot(lmefit1,resid(.)~fitted(.)|SN)
plot(lmefit1,resid(.,type="p")~fitted(.)|SN) #Standardized Residuals
plot(lmefit1,SN~resid(.))
plot(lmefit1,resid(.)~Temp|SN)
plot(lmefit1,resid(.)~Precip|SN)
qqnorm(lmefit1,~resid(.),id=.05,adj=-0.75)
qqnorm(lmefit1,~resid(.)|SN)
lmefit2=update(lmefit1,weights=varPower())
plot(lmefit2,resid(.,type="p")~fitted(.)|SN)
plot(lmefit2,SN~resid(.,type="p"))
plot(lmefit2,resid(.,type="p")~Temp|SN)
plot(lmefit2,resid(.,type="p")~Precip|SN)
qqnorm(lmefit2,~resid(.,type="p"),id=.05,adj=-0.75)
qqnorm(lmefit2,~resid(.,type="p")|SN)
plot(ACF(lmefit2,maxLag=10),alpha=0.05)
lmefit3=update(lmefit2,corr=corARMA(form=~1|SN,p=0,q=4))
plot(ACF(lmefit3,maxLag=10),alpha=0.05)
anova(lmefit2,lmefit3)
StationNames=as.character(unique(GD1$SN))
par(mfrow=c(3,5))
par(cex.axis=1.10,cex.lab=1.25)
for(i in 1:15){
	year=GD1$Year[GD1$SN==StationNames[i]]
	y0=GD1$Discharge[GD1$SN==StationNames[i]]
	y2=lmefit3$fitted[,2][lmefit3$groups==StationNames[i]]
	y1=lmefit3$fitted[,1][lmefit3$groups==StationNames[i]]
	plot(year,y0,type="l",main=StationNames[i],ylab=expression(log),ylim=c(0,10))
	lines(year,y1, type="l",col="green")
	lines(year,y2,type="l", col="blue")
}

########################
# Relating Normalized Discharge to Temprature and Precipitation for 15 randomly selected streams in MT observed for over 50 years
##############
D=read.csv(file="Discharge.csv",header=T)
library(nlme)
GD1=groupedData(Normalized~Temp|SN,data=D)
plot(GD1)
GD2=groupedData(Normalized~Precip|SN,data=D)
plot(GD2)
lmListfit=lmList(Normalized~ Temp + Precip + WS,data=GD1)
plot(intervals(lmListfit))
lmefit1=lme(Normalized ~ Temp + Precip + WS,data=GD1, random=list(SN=pdDiag(~1+Precip + Temp)))
plot(augPred(lmefit1),grid=T)
cf=compareFits(coef(lmefit1),coef(lmListfit))
plot(cf,mark=fixef(lmefit1))
plot(comparePred(lmefit1,lmListfit,length.out=2))
plot(lmefit1,resid(.)~fitted(.)|SN)
plot(lmefit1,resid(.,type="p")~fitted(.)|SN) #Standardized Residuals
plot(lmefit1,SN~resid(.))
plot(lmefit1,resid(.)~Temp|SN)
plot(lmefit1,resid(.)~Precip|SN)
qqnorm(lmefit1,~resid(.),id=.05,adj=-0.75)
qqnorm(lmefit1,~resid(.)|SN)
lmefit2=update(lmefit1,weights=varPower())
plot(lmefit2,resid(.,type="p")~fitted(.)|SN)
plot(lmefit2,SN~resid(.,type="p"))
plot(lmefit2,resid(.,type="p")~Temp|SN)
plot(lmefit2,resid(.,type="p")~Precip|SN)
qqnorm(lmefit2,~resid(.,type="p"),id=.05,adj=-0.75)
qqnorm(lmefit2,~resid(.,type="p")|SN)
plot(ACF(lmefit2,maxLag=10),alpha=0.05)
lmefit3=update(lmefit2,corr=corARMA(form=~1|SN,p=0,q=4))
plot(ACF(lmefit3,maxLag=10),alpha=0.05)
anova(lmefit2,lmefit3)
StationNames=as.character(unique(GD1$SN))
par(mfrow=c(3,5))
par(cex.axis=1.10,cex.lab=1.25)
for(i in 1:15){
	year=GD1$Year[GD1$SN==StationNames[i]]
	y0=GD1$Normalized[GD1$SN==StationNames[i]]
	y2=lmefit3$fitted[,2][lmefit3$groups==StationNames[i]]
	y1=lmefit3$fitted[,1][lmefit3$groups==StationNames[i]]
	plot(year,y0,type="l",main=StationNames[i],ylim=c(0,3.5),ylab="CFS")
	lines(year,y1, type="l",col="green")
	lines(year,y2,type="l", col="blue")
}