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