library(Matrix)
library(car)
library(nlme)

#########
# Fundtion to Print the error covariance matrix 
########
corandcov <- function(glsob,...){
  corm <- corMatrix(glsob$modelStruct$corStruct)[[1]]
  print(corm)
  varests <- glsob$modelStruct$varStruct 
  print(varests) 
  #print(glsob$sigma^2*bdiag(1,diag(varests))%*%corm%*%bdiag(1,diag(varests)))
  }
D=read.csv(file="www.math.umt.edu/harrar/Stat542/RAAP1.CSV",header=T)

GD=groupedData(FEV1~Time|Drug,data=D,inner=~Patient)
D$AllT=factor("All Treatments");
GD1=groupedData(FEV1~Time|Drug/Patient,data=D,outer=~AllT)
plot(GD,inner=T,layout=c(3,1))
plot(GD,layout=c(3,1))
plot(GD1,display=1,collapse=1,layout=c(3,1),outer=T)
GD$DP=factor(interaction(GD$Drug,GD$Patient))
lmefitUN=gls(FEV1~factor(Time)*Drug,data=GD,method="REML",corr=corSymm(form=~1|DP),weights=varIdent(form=~1|Time))
anova(lmefitUN)
lmefitUN
corandcov(lmefitUN)
c=corMatrix(lmefitUN$modelStruct$corStruct)[[1]]
v=c(1.000000, 1.066242, 1.041196, 1.042749 ,1.128106, 1.039408, 1.048692 ,1.052557)
Sigma=lmefitUN$sigma^2*(diag(v)%*%c%*%diag(v))
plot(0:7,Sigma[,1],type="b",ylim=c(0.3,0.6))
lines(0:6,Sigma[2:8,2],type="b",col="red")
lines(0:5,Sigma[3:8,3],type="b",col="green")
lines(0:4,Sigma[4:8,4],type="b",col="cyan")
lines(0:3,Sigma[5:8,5],type="b",col="orange")
lines(0:2,Sigma[6:8,6],type="b",col="purple")
lines(0:1,Sigma[7:8,7],type="b",col="blue")
lmefitCS=gls(FEV1~factor(Time)*Drug,data=GD,method="REML",corr=corCompSymm(form=~1|DP))
anova(lmefitCS)
summary(lmefitCS)
corandcov(lmefitCS)

#########
# Selecting the  Appropriate Covaraince Model
####
lmefitUN=gls(FEV1~factor(Time)*Drug,data=GD,method="REML",corr=corSymm(form=~1|DP),weights=varIdent(form=~1|Time))
summary(lmefitUN)
corandcov(lmefitUN)
lmefitCS=gls(FEV1~factor(Time)*Drug,data=GD,method="REML",corr=corCompSymm(form=~1|DP))
summary(lmefitSC)
corandcov(lmefitCS)
#Unidentifiable model
#lmefitCSpCS=lme(FEV1~factor(Time)*Drug,data=GD,method="REML",random=~1|DP,corr=corCompSymm(form=~1|DP))
lmefitAR=gls(FEV1~factor(Time)*Drug,data=GD,method="REML",corr=corARMA(form=~1|DP,p=1,q=0))
summary(lmefitAR)
corandcov(lmefitAR)
lmefitARpCS=lme(FEV1~factor(Time)*Drug,data=GD,method="REML",random=~1|DP,corr=corARMA(form=~1|DP,p=1,q=0))
summary(lmefitARpCS)
corandcov(lmefitARpCS)
lmefitToep=lme(FEV1~factor(Time)*Drug,data=GD,method="REML",random=~1|DP,corr=corARMA(form=~1|DP,p=7,q=0))
summary(lmefitToep)
corandcov(lmefitToep)
anova(lmefitUN, lmefitAR, lmefitARpCS, lmefitToep,lmefitCS)
#########################
# The Solling Roof Project
########
D=read.csv("http://www.math.umt.edu/harrar/Stat542/So4.csv",header=TRUE)
GD=groupedData(Conc~Year|Sample,data=D,outer=~Tx)
plot(GD,outer=T)
glsfitUN=gls(Conc~factor(Year)*Tx,data=GD,method="REML",corr=corSymm(form=~1|Sample),weights=varIdent(form=~1|Year))
summary(glsfitUN)
c=corMatrix(glsfitUN$modelStruct$corStruct)[[1]]
v=c(1.0000000, 0.8465522, 0.7647957, 0.6029932, 0.4626671, 0.3797792, 0.3790406 )
Sigma=glsfitUN$sigma^2*(diag(v)%*%c%*%diag(v))
plot(0:6,Sigma[,1],type="b",ylim=c(0,30))
lines(0:5,Sigma[2:7,2],type="b",col="red")
lines(0:4,Sigma[3:7,3],type="b",col="green")
lines(0:3,Sigma[4:7,4],type="b",col="cyan")
lines(0:2,Sigma[5:7,5],type="b",col="orange")
lines(0:1,Sigma[6:7,6],type="b",col="purple")

###############
# Respiratory Ability Adjusting for Covariate
#######
lmefitUN.ABL=gls(FEV1~BASEFEV1+ factor(Time)*Drug,data=GD,method="REML",corr=corSymm(form=~1|DP),weights=varIdent(form=~1|Time))
anova(lmefitUN.ABL)
lmefitUN.ABL
corandcov(lmefitUN.ABL)
c=corMatrix(lmefitUN.ABL$modelStruct$corStruct)[[1]]
v=c(1.000000, 1.069884, 1.056962, 1.146655, 1.122590, 1.066309, 1.092259, 1.148241 )
Sigma.ABL=lmefitUN.ABL$sigma^2*(diag(v)%*%c%*%diag(v))
plot(0:7,Sigma.ABL[,1],type="b",ylim=c(0.12,0.32))
lines(0:6,Sigma.ABL[2:8,2],type="b",col="red")
lines(0:5,Sigma.ABL[3:8,3],type="b",col="green")
lines(0:4,Sigma.ABL[4:8,4],type="b",col="cyan")
lines(0:3,Sigma.ABL[5:8,5],type="b",col="orange")
lines(0:2,Sigma.ABL[6:8,6],type="b",col="purple")
lines(0:1,Sigma.ABL[7:8,7],type="b",col="blue")
lmefitARpCS.ABL=lme(FEV1~BASEFEV1+factor(Time)*Drug,data=GD,method="REML",random=~1|DP,corr=corARMA(form=~1|DP,p=1,q=0))
anova(lmefitARpCS.ABL)
summary(lmefitARpCS.ABL)
interaction.plot(GD$Time,GD$Drug, lmefitARpCS.ABL$fitted[,1]-lmefitARpCS.ABL$coefficients$fixed[[2]]*GD$BASEFEV1)


########################
#Unequally Spaced Repeated Measures
############
D=read.table(file="www.math.umt.edu/harrar/Stat542/Heart.dat",header=T)
GD=groupedData(HeartRate~Time|Id,data=D,outer=~Drug)
plot(GD,outer=T)
lmefit=lme(HeartRate~ Baseline + Drug*factor(Time),data=GD, random=~1|factor(Id), correlation=corGaus(form=~Time|factor(Id)),method="ML")