######################################################
#ORANGE TREE EXAMPLE########################################
######################################################

library(nlme)

# Orange tree data are included in nlme library
# look at Orange groupedData object:
Orange[1:3,]

# The following nnlsList function fits the fixed effects logistic model 
# to the data from each tree separately via nonlinear least squares
mOran.lis <- nlsList(circumference ~ SSlogis(age,Asym,th2,th3)|Tree, data=Orange)
coef(mOran.lis)

# Now fit one fixed-effect model to all data with a common asymptote,
# rate and shape parameter for all trees. Gives solid curve in the
# plot on p.77, The function gnls uses generalized nonlinear least
# squares as its fitting method, which means nonlinear least squares
# allowing heteroscedasticity and correlation among the errors
# Here, though, the model is just being fit by ordinary least squares
m1Oran.gnls <- gnls(circumference ~ SSlogis(age,Asym,th2,th3), data=Orange,
 params= list(Asym~1,th2~1,th3~1))

# plot of residuals from m1Oran.gnls shows cone shape suggesting 
# heteroscedasticity, but that's not the problem. Within-tree variance
# seems to be stable. Problem is that mean response from model is
# off target for each individual tree.
plot(m1Oran.gnls)

# now fit model with different asymptote parameters for each tree
m2Oran.gnls <- gnls(circumference ~ SSlogis(age,Asym,th2,th3), data=Orange,
 params= list(Asym~Tree-1,th2~1,th3~1), start=c(200,0,0,0,0,coef(m1Oran.gnls)[2:3]))
plot(m2Oran.gnls) #residual plot looks much better, but trees effects should
                  # be random


# Fit the NLMM model with random tree-specific asymptotes
m1Oran.nlme <- nlme(circumference ~ SSlogis(age,Asym,th2,th3), data=Orange,
 fixed= Asym+th2+th3~1,random= Asym~1,start=coef(m1Oran.gnls)) 

fixef(m1Oran.nlme) #compare fixed effect estimates
coef(m1Oran.gnls)

summary(m1Oran.nlme)
AIC(m1Oran.nlme,m1Oran.gnls,m2Oran.gnls) #compare AICs for the mixed, fixed models

# compare residual plots 
plot(m1Oran.gnls,grid=F,main="Residuals vs Fitteds, Model m1Oran.gnls")
plot(m2Oran.gnls,grid=F,main="Residuals vs Fitteds, Model m2Oran.gnls")
plot(m1Oran.nlme,grid=F,
 main="Residuals vs Fitteds, Mixed Effects Model (m1Oran.nlme)")

#Here's a nice plot of the predicted values from m1Oran.nlme at the
#poulation level (blue) and tree-level (pink)
plot( augPred(m1Oran.nlme,level=0:1),layout=c(5,1))

# Pink curves have asymptotes given by theta1hat+ the following BLUPs of
# the random effects
ranef(m1Oran.nlme)


plot(ACF(m2Oran.gnls,maxLag=6,form=~1|Tree),alpha=.05,
 main="ACF, Model m2Oran.gnls")
plot(ACF(m1Oran.nlme,maxLag=6),alpha=.05,main="ACF, Model m1Oran.nlme")

######################################################
# TPH EXAMPLE - PMRC SITEPREP STUDY ##################
######################################################

###Build the data set again
###Details not important
siteprep.frm <- read.table(file="c:/courses/shortcourse/TempleInland05/siteprep.dat",header=T,
  na.strings=c("."))
siteprep.frm$install <- factor(siteprep.frm$install)
siteprep.frm$site <- siteprep.frm$install
siteprep.frm$install <- NULL
siteprep.frm$plot <- factor(siteprep.frm$plot)
siteprep.frm$allplot <- factor(siteprep.frm$allplot)
siteprep.frm$soil <- factor(siteprep.frm$soil)
siteprep.frm$herb <- factor(siteprep.frm$herb)
siteprep.frm$fert <- factor(siteprep.frm$fert)
siteprep.frm$bed <- factor(siteprep.frm$bed)
siteprep.frm$chop <- factor(siteprep.frm$chop)
siteprep.frm$burn <- factor(siteprep.frm$burn)

siteprep.frm$agefac <- factor(siteprep.frm$agefac)
siteprep.frm$measocc <- factor(siteprep.frm$measocc)
siteprep.frm$soil1 <- as.numeric(siteprep.frm$soil==1)
siteprep.frm$soil2 <- as.numeric(siteprep.frm$soil==2)
siteprep.frm$soil3 <- as.numeric(siteprep.frm$soil==3)
siteprep.frm$soil4 <- as.numeric(siteprep.frm$soil==4)

siteprep.tph <- groupedData(tph100~age|site/plot,
  data=siteprep.frm[siteprep.frm$yname=="hd1",],
  labels=list(x="Age",y="Trees/Hectare"),units=list(x="(yrs)",y="(100s of trees)"))
###Done reading in and formatting data

#look at data
siteprep.tph[1:3,]

m1tph.fpl <- nlme(tph100~SSfpl(age,th1,th2,th3,th4),data=siteprep.tph,
  fixed=list(th1~1,th2~1,th3~1,th4~1),random=list(site=pdDiag(th1+th2+th3~1),
  plot=pdDiag(th1+th2+th3~1)),start=c(10,8,10,1))
#fixef(m1tph.fpl)
summary(m1tph.fpl)

plot(augPred(m1tph.fpl,level=0:2),layout=c(4,5))

plot(ranef(m1tph.fpl,aug=T,level=2),form=th1~itph100)

m2tph.fpl <- update(m1tph.fpl,
  fixed=list(th1~itph100,th2~1,th3~1,th4~1),
  start=c(0,1,fixef(m1tph.fpl)[2:4]))
anova(m2tph.fpl,type="marginal")
anova(m1tph.fpl,m2tph.fpl)

plot(ranef(m2tph.fpl,aug=T,level=2),form=th1.(Intercept)~trt)
plot(ranef(m2tph.fpl,aug=T,level=2),form=th2~trt)
plot(ranef(m2tph.fpl,aug=T,level=2),form=th3~trt)

m3tph.fpl <- update(m2tph.fpl,
  fixed=list(th1~itph100,th2~trt,th3~1,th4~1),
  start=c(fixef(m1tph.fpl)[1:3],rep(0,10),fixef(m2tph.fpl)[4:5]))
anova(m3tph.fpl,type="marginal")
anova(m3tph.fpl,m2tph.fpl)

plot(ranef(m3tph.fpl,aug=T,level=2),form=th1.(Intercept)~trt)
plot(ranef(m3tph.fpl,aug=T,level=2),form=th3~trt)

m4tph.fpl <- update(m3tph.fpl,
  fixed=list(th1~itph100,th2~trt,th3~trt,th4~1),
  start=c(fixef(m3tph.fpl)[1:14],rep(0,10),fixef(m3tph.fpl)[15]))
anova(m4tph.fpl,type="marginal")
anova(m3tph.fpl,m4tph.fpl)

m5tph.fpl <- update(m4tph.fpl,
  fixed=list(th1~itph100+trt,th2~trt,th3~trt,th4~1),
  start=c(fixef(m4tph.fpl)[1:2],rep(0,10),fixef(m4tph.fpl)[3:25]))
anova(m5tph.fpl,type="marginal")
anova(m5tph.fpl,m4tph.fpl)

plot(ranef(m5tph.fpl,aug=T,level=2),form=th1.(Intercept)~si)
plot(ranef(m5tph.fpl,aug=T,level=2),form=th2.(Intercept)~si)
plot(ranef(m5tph.fpl,aug=T,level=2),form=th3.(Intercept)~si)

plot(ranef(m5tph.fpl,aug=T,level=2),form=th1.(Intercept)~soil)
plot(ranef(m5tph.fpl,aug=T,level=2),form=th2.(Intercept)~soil)
plot(ranef(m5tph.fpl,aug=T,level=2),form=th3.(Intercept)~soil)

#examine some overall residual plots versus potential covariates
plot(m5tph.fpl,soil~resid(.,type="n"))
plot(m5tph.fpl,resid(.,type="n")~si)
plot(m5tph.fpl,herb~resid(.,type="n"))
plot(m5tph.fpl,fert~resid(.,type="n"))
plot(m5tph.fpl,bed~ resid(.,type="n"))
plot(m5tph.fpl,chop~resid(.,type="n"))
plot(m5tph.fpl,burn~resid(.,type="n"))

# residuals vs fitteds
plot(m5tph.fpl,resid(.,type="n")~fitted(.))
plot(m5tph.fpl,resid(.,type="n")~fitted(.)|site)

plot(m5tph.fpl,resid(.,type="n")~si)

m6tph.fpl <- update(m5tph.fpl, weights=varPower(form=~si),
  start=fixef(m5tph.fpl))
anova(m5tph.fpl,m6tph.fpl)

summary(m6tph.fpl)

plot(ACF(m6tph.fpl,resType="n", form=~1|site/plot),alpha=.05)

#m7tph.fpl <- update(m6tph.fpl, #This did not converge
#  corr=corAR1(value=-.4,form= ~1|site/plot),
#  start=fixef(m6tph.fpl))

#add correlation in two steps
m7tph.fpl <- update(m5tph.fpl,
  corr=corAR1(value=-.4,form= ~1|site/plot),
  start=fixef(m5tph.fpl))
m8tph.fpl <- update(m7tph.fpl,
  weights=varPower(form=~si),
  start=fixef(m7tph.fpl))
anova(m8tph.fpl,m6tph.fpl,m7tph.fpl)

plot(ACF(m8tph.fpl,resType="n", form=~1|site/plot),alpha=.05)

plot(m6tph.fpl,resid(.,type="n")~age)

m9tph.fpl <- update(m6tph.fpl,
  random=list(site=pdDiag(th1+th2+th3~1),
  plot=pdDiag(th1+th3~1)),
  start=fixef(m5tph.fpl))
anova(m9tph.fpl,m6tph.fpl)

anova(m9tph.fpl,type="marginal")
summary(m9tph.fpl)