###################################################### #MFA EXAMPLE######################################## ###################################################### library(nlme) #first read mfa data into a data frame and look at it mfa.frm <- read.table(file="e:/TempleInland05/mfadata.dat",header=T) mfa.frm[1:3,] mfa <- groupedData(mfa~diskht|tree,data=mfa.frm, labels=list(x="Height on stem", y="Whole-disk cross-sectional microbril angle"), units=list(x="(m)",y="(deg)"), order.groups=F)#make mfa into a groupedData object mfa[1:3,] #look at it mfa\$disk <- as.factor(mfa\$disk) mfa\$regionno <- as.factor(mfa\$regionno) mfa\$stand <- as.factor(mfa\$stand) mfa\$tree <- as.factor(mfa\$tree) lapply(mfa,is.factor) lapply(mfa,is.ordered) #make sure factors are all treated as unordered #plot the data plot(mfa,outer= ~regname,key=F) mfa.lme1 <- lme(mfa ~ regname + diskht -1 , data=mfa, random= ~1|tree) summary(mfa.lme1) mfa.lme2 <- update(mfa.lme1, random= ~diskht|tree) summary(mfa.lme2) # define a contrast matrix A to test the hypothesis A*beta=0 # with a Wald-based F test # Here, the hypothesis is that all regions have the same intercept # or, equivalently, the same mean MFA A <- matrix( c(1,-1,0,0,0, 0,1,-1,0,0, 0,0,1,-1,0),nrow=3,byrow=T) A fixef(mfa.lme1) anova(mfa.lme1,L=A,type="marginal") # use marginal, or "Type 3" SSs # Notice above that the denominator df in this F test is based upon # the number of trees (59), not the total number of observations # Trees here are playing the role of whole plots. With respect to # region, which is a whole plot factor, repeated measures withom each # tree are pseudoreplicates # Another way to get the same result is to reparameterize # the model as y_ijk = mu + alpha_i + beta*height + b_ij + e_ijk # and then get the anova table as follows: mfa.lme1a <- lme(mfa ~ regname + diskht, data=mfa, random= ~1|tree) anova( mfa.lme1a, type="marginal") # or can use the Terms argument to just test region anova( mfa.lme1a, type="marginal", Terms=2) # We can get t-based confidence intervals for the means in each region # (or any fixed effects) # via the intervals() function. The intervals given for the variance # components here are Wald-based too. They are very approximate though # and should be treated with caution intervals(mfa.lme1) # Now consider random intercept model with quadratic effect of height # on MFA. Fit with ML mfa.lme3.ML <- lme(mfa ~ regname + diskht + I(diskht^2) -1 , data=mfa, random= ~1|tree,method="ML") #full model mfa.lme1.ML <- update(mfa.lme3.ML, fixed= ~ regname + diskht -1 )#partial model anova(mfa.lme3.ML,mfa.lme1.ML) # do LRT of diskht^2 2*(logLik(mfa.lme3.ML)[1]-logLik(mfa.lme1.ML)[1]) summary(mfa.lme3.ML)\$tTable #gives Wald-based t-test on diskht^2, alternative to LRT anova(mfa.lme3.ML,Terms=3) #Wald F test, just square of previous #refit random intercept, quadratic in height model (model mfa.lme3.ML) #with REML, and then consider whether linear height effect, quadratic #height effects should have random effects so that they vary form tree #to tree mfa.lme3 <- update(mfa.lme3.ML, method="REML" ) mfa.lme4 <- update(mfa.lme3, random=~diskht|tree ) mfa.lme5 <- update(mfa.lme3, random=~diskht+I(diskht^2)|tree ) anova(mfa.lme3,mfa.lme4,mfa.lme5) # Returning to the quadratic in height, random intercept model # let's examine the predicted random effects ranef(mfa.lme3) plot(augPred(mfa.lme3,level=0:1)) #Plots data from each tree along with #estimated/predicted mean value (blue) #and mean+predicted random effect #for each tree (pink) #Now do same for just Hilly region (trees 51-59) ranef(mfa.lme3)[51:59,] plot(augPred(mfa.lme3,level=0:1)[augPred(mfa.lme3,level=0:1)\$.groups>50,], main="Data from Hilly Region only") # residual plots plot(mfa.lme3) #standardized (level 1) resids vs fitteds plot(mfa.lme3,resid(.,type="p")~fitted(.),id=.01,adj=-.3,abline=0) #same thing plot(mfa.lme3,resid(.,type="p")~fitted(.)|regname,id=.01,adj=-.3,abline=0)#same thing plotted by region plot(mfa.lme3,regname~resid(.,type="p") ,abline=0) #same thing but boxplots plot(mfa.lme3,mfa~fitted(.)) #response vs fitteds (should be straight line) qqnorm(mfa.lme3, ~resid(.)) #normal quantile-quantile plot of resids qqnorm(mfa.lme3, ~resid(.)|regname) #normal q-q plots by region plot(mfa.lme3,resid(.,type="p")~diskht,id=.01,adj=-.3,abline=0) #resids vs diskht plot(mfa.lme3,resid(.,type="p")~diskht|regname,id=.01,adj=-.3,abline=0) #resids vs diskht # predicted random effect plots plot(ranef(mfa.lme3)) #compare random effects across trees qqnorm(mfa.lme3,~ranef(.)) #qq plot of random effects (checks normality) # Let linear diskht effect differ across regions mfa.lme6 <- update(mfa.lme3, fixed= ~ -1+regname+diskht+I(diskht^2)+regname:diskht) summary(mfa.lme6) anova(mfa.lme6,type="marginal") #this model is better than mfa.lme3 # Let linear and quadratic diskht effects differ across regions mfa.lme7 <- update(mfa.lme3, fixed= ~ -1+regname+diskht+I(diskht^2)+regname:diskht+regname:I(diskht^2)) summary(mfa.lme7) anova(mfa.lme7,type="marginal") #this model is no better than mfa.lme6 #look at diagnostics for mfa.lme6 plot(mfa.lme6,resid(.,type="p")~fitted(.),id=.01,adj=-.3,abline=0) #same thing plot(mfa.lme6,resid(.,type="p")~fitted(.)|regname,id=.01,adj=-.3,abline=0)#same thing plotted by region qqnorm(mfa.lme6, ~resid(.)) #normal quantile-quantile plot of resids plot(mfa.lme6,resid(.,type="p")~diskht|regname,id=.01,adj=-.3,abline=0) #resids vs diskht # model still has problems # Now just treat diskht as a factor and use two-way anova model with # random tree effects mfa.lme8 <- update(mfa.lme6, fixed= ~ regname+factor(diskht)+regname:factor(diskht)) summary(mfa.lme8) anova(mfa.lme8,type="marginal") #check residual plots again plot(mfa.lme8,resid(.,type="p")~fitted(.),id=.01,adj=-.3,abline=0) #same thing plot(mfa.lme8,resid(.,type="p")~fitted(.)|regname,id=.01,adj=-.3,abline=0)#same thing plotted by region qqnorm(mfa.lme8, ~resid(.)) #normal quantile-quantile plot of resids plot(mfa.lme8,resid(.,type="p")~diskht|regname,id=.01,adj=-.3,abline=0) #resids vs diskht #these finally look ok