R routines for modeling coarse woody debris dynamics

By |2018-02-12T10:40:12+00:00February 14th, 2016|My Blog, Programming Language|0 Comments

bioAgro = read.table(“BioAgro.txt”,header=T,sep=”\t”)
bioAgro$Block=as.factor(bioAgro$Block)
str(bioAgro)
summary(bioAgro)
fhg

odel1=lm(Leaf_mass_g~Tmt + Spp + Interaction + Block,data=bioAgro)
summary(model1)
anova(model1)

library(nlme)
model1.gls =gls(Leaf_mass_g ~ Tmt + Spp + Interaction + factor(Block),
method=”REML”,data=bioAgro)
summary(model1.gls)
#Addressing heterogeneity in Species level
model1.gls1 =gls(Leaf_mass_g ~ Tmt + Spp + Interaction + Block,
weights=varIdent(form=~1|Spp),data=bioAgro)
model1.gls2=gls(Leaf_mass_g ~ Tmt + Spp + Interaction + factor(Block),
weights=varIdent(form=~1|Tmt),data=bioAgro)
model1.gls3=gls(Leaf_mass_g ~ Tmt + Spp + Interaction + factor(Block),
weights=varIdent(form=~1|Interaction),data=bioAgro)
model1.gls4=gls(Leaf_mass_g ~ Tmt + Spp + Interaction + Block,
weights=varIdent(form=~1|Spp*Tmt),method=”REML”,data=bioAgro)
summary(model1.gls4)
anova(model1.gls,model1.gls4)

boxplot(predict(model1.gls4)~ Tmt + Spp ,
data =bioAgro, cex.axis=0.75,xlab=”Treatment.Species.Interaction”,
ylab=”Dry Leaf Mass (g)”)

boxplot(predict(model1)~Tmt + Spp+ Interaction,data=bioAgro)

library(multcomp)
## Run the following lines. These introduce methods for ‘gls’ objects.
model.matrix.gls <- function(object, …) {
model.matrix(terms(object), data = getData(object), …)
}
model.frame.gls <- function(object, …) {
model.frame(formula(object), data = getData(object), …)
}
terms.gls <- function(object, …) {
terms(model.frame(object), …)
}
# The line below gives pairwise comparisons now. Note that the above
# performs t-tests for all pairwise differences.

multCompTukey <- glht(glsModel1, linfct = mcp(site = “Tukey”))
a=glht(model1.gls4,linfct=mcp(Spp=”Tukey”))
summary(a)
#Model1.gls4 is so far the best with AIC = 81.39

model1.gls5=gls(Leaf_mass_g ~ Tmt + Spp + Interaction + factor(Block),
weights=varIdent(form=~1|Spp*Interaction),data=bioAgro)
model1.gls6=gls(Leaf_mass_g ~ Tmt + Spp + Interaction + factor(Block),
weights=varIdent(form=~1|Spp*Block),data=bioAgro)
model1.gls7=gls(Leaf_mass_g ~ Tmt + Spp + Interaction + Block,
weights=varIdent(form=~1|Spp*Tmt*Block),method=”ML”,
data=bioAgro)

summary(list(bioAgro$Leaf_mass_g,bioAgro$Tmt,bioAgro$Spp,bioAgro$Block,
bioAgro$Interaction))
anova(model1.gls4,model1.gls1,model1)
AIC(model1.gls1,model1.gls2,model1.gls3,model1.gls4,model1.gls5,
model1.gls6,model1.gls7)

model1.gls =gls(Leaf_mass_g ~ Tmt + Spp + Interaction + factor(Block),
weights=varIdent(form=~1|Spp*Tmt),data=bioAgro)
summary(model1.gls2)

modell.gls3=gls(Leaf_mass_g ~ Tmt + Spp + Interaction + factor(Block),
weights=varIdent(form=~1|Spp*Interaction),data=bioAgro)
summary(modell.gls3)

modell.gls4=gls(Leaf_mass_g ~ Tmt + Spp + Interaction + Block,
weights=varIdent(form=~1|Spp*Interaction*Tmt),data=bioAgro)
summary(modell.gls4)

summary(model1.gls1)
aov(model1, model1.gls1)
AIC(model1,model1.gls,model1.gls1)

op<-par(mfrow=c(2,3),mar=c(4,4,1,2))
plot(model1.gls4,col=2,add.smooth=F, which=1)
residual_1=resid(model1.gls4,type=”normalized”)
hist(residual_1,xlab=”Residuals”,main=””)
plot(bioAgro$Tmt,residual_1,xlab=”Treatment”,ylab=”Residuals”)
plot(bioAgro$Spp,residual_1,xlab=”Species”,ylab=”Residuals”)
plot(factor(bioAgro$Block),residual_1,xlab=”Block”,ylab=”Residuals”)
plot(bioAgro$Interaction,residual_1,xlab=”Interaction”,ylab=”Residuals”)
par(op)
plot(bioAgro$Tmt, bioAgro$Leaf_mass_g)
plot(bioAgro$Spp, bioAgro$Leaf_mass_g)
plot(bioAgro$Block, bioAgro$Leaf_mass_g)
plot(bioAgro$Interaction, bioAgro$Leaf_mass_g)
bartlett.test(residual_1,bioAgro$Tmt)
bartlett.test(residual_1,bioAgro$Spp)
bartlett.test(residual_1,factor(bioAgro$Block))
bartlett.test(residual_1,bioAgro$Interaction)
hist(residual_1[bioAgro$Tmt])

About the Author:

Leave A Comment