###################################################### #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)