Category: Programming Language

R routines for modeling coarse woody debris dynamics

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

R-ggplot2: Adding both vertical and horizontal error bars in each point of a scatter plot

Sometimes it might be necessary to show both horizontal and vertical error bars/point in a scatterplot if both the dependent (Y) and independent (X) variables are averages of some other values. This is my personal note to keep things organized and ease of future use. I’d be happy if anyone is benefitted.

Suppose, we have collected data of biomass and height of some plants and recorded the species names. Now we want to see how biomass is related to plant’s height across species. To do so lets plot the average biomass in Y-axis and average height in X-axis in a scatter plot and fit a linear regression line across the species. We also want to draw error bars for each mean biomass and mean height pair. Here is how we can do:


# Create the dataset
biomass = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
height = c(0.15,.22,.34,.41,.5,.62,.71,.81,.92,1.01,1.1,1.25,1.32,1.41,1.52)
species = c('a','a','a','a','a','b','b','b','b','b','c','c','c','c','c')

dummyData = data.frame(biomass,height,species)
dummyData

library(Rmisc)

# Calculating SE (standard error) which will be used for error bars
bmassStat = summarySE(dummyData,measurevar=”biomass”,groupvars = “species”)
bmassStat
heightStat = summarySE(data=dummyData,measurevar = “height”,groupvars = “species”)
heightStat

# Combining variables and statistics
bmassHt = data.frame(bmassStat$biomass,bmassStat$se,heightStat$height,heightStat$se,heightStat$species)
names(bmassHt)=c(“avgBiomass”,”biomassSE”,”avgHeight”,”heightSE”,”species”)
bmassHt

library(ggplot2)
ggplot(data=bmassHt,aes(x=avgHeight,y=avgBiomass,color=species))+
geom_smooth(method=’lm’,se=FALSE,color=”black”)+
geom_point()+
geom_errorbar(aes(ymax=avgBiomass+biomassSE,ymin=avgBiomass-biomassSE,width=0))+
geom_errorbarh(aes(xmax=avgHeight+heightSE,xmin=avgHeight-heightSE))

# geom_errobarh creates the horizontal error bar for the variable in the X-axis (in this case avgHeight)

So, this is the final graph:
Horizontal & Vertical Error Bars

Please let me know if you have any questions or suggestions.

Sometimes it might be necessary to show both horizontal and vertical error bars/point in a scatterplot if both the dependent (Y) and independent (X) variables are averages of some other values. This is my personal note to keep things organized and ease of future use. I’d be happy if anyone is benefitted.

Suppose, we have collected data of biomass and height of some plants and recorded the species names. Now we want to see how biomass is related to plant’s height across species. To do so lets plot the average biomass in Y-axis and average height in X-axis in a scatter plot and fit a linear regression line across the species. We also want to draw error bars for each mean biomass and mean height pair. Here is how we can do:


# Create the dataset
biomass = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
height = c(0.15,.22,.34,.41,.5,.62,.71,.81,.92,1.01,1.1,1.25,1.32,1.41,1.52)
species = c('a','a','a','a','a','b','b','b','b','b','c','c','c','c','c')

dummyData = data.frame(biomass,height,species)
dummyData

library(Rmisc)

# Calculating SE (standard error) which will be used for error bars
bmassStat = summarySE(dummyData,measurevar=”biomass”,groupvars = “species”)
bmassStat
heightStat = summarySE(data=dummyData,measurevar = “height”,groupvars = “species”)
heightStat

# Combining variables and statistics
bmassHt = data.frame(bmassStat$biomass,bmassStat$se,heightStat$height,heightStat$se,heightStat$species)
names(bmassHt)=c(“avgBiomass”,”biomassSE”,”avgHeight”,”heightSE”,”species”)
bmassHt

library(ggplot2)
ggplot(data=bmassHt,aes(x=avgHeight,y=avgBiomass,color=species))+
geom_smooth(method=’lm’,se=FALSE,color=”black”)+
geom_point()+
geom_errorbar(aes(ymax=avgBiomass+biomassSE,ymin=avgBiomass-biomassSE,width=0))+
geom_errorbarh(aes(xmax=avgHeight+heightSE,xmin=avgHeight-heightSE))

# geom_errobarh creates the horizontal error bar for the variable in the X-axis (in this case avgHeight)

So, this is the final graph:
Horizontal & Vertical Error Bars

Please let me know if you have any questions or suggestions.

Sometimes it might be necessary to show both horizontal and vertical error bars/point in a scatterplot if both the dependent (Y) and independent (X) variables are averages of some other values. This is my personal note to keep things organized and ease of future use. I’d be happy if anyone is benefitted.

Suppose, we have collected data of biomass and height of some plants and recorded the species names. Now we want to see how biomass is related to plant’s height across species. To do so lets plot the average biomass in Y-axis and average height in X-axis in a scatter plot and fit a linear regression line across the species. We also want to draw error bars for each mean biomass and mean height pair. Here is how we can do:


# Create the dataset
biomass = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
height = c(0.15,.22,.34,.41,.5,.62,.71,.81,.92,1.01,1.1,1.25,1.32,1.41,1.52)
species = c('a','a','a','a','a','b','b','b','b','b','c','c','c','c','c')

dummyData = data.frame(biomass,height,species)
dummyData

library(Rmisc)

# Calculating SE (standard error) which will be used for error bars
bmassStat = summarySE(dummyData,measurevar=”biomass”,groupvars = “species”)
bmassStat
heightStat = summarySE(data=dummyData,measurevar = “height”,groupvars = “species”)
heightStat

# Combining variables and statistics
bmassHt = data.frame(bmassStat$biomass,bmassStat$se,heightStat$height,heightStat$se,heightStat$species)
names(bmassHt)=c(“avgBiomass”,”biomassSE”,”avgHeight”,”heightSE”,”species”)
bmassHt

library(ggplot2)
ggplot(data=bmassHt,aes(x=avgHeight,y=avgBiomass,color=species))+
geom_smooth(method=’lm’,se=FALSE,color=”black”)+
geom_point()+
geom_errorbar(aes(ymax=avgBiomass+biomassSE,ymin=avgBiomass-biomassSE,width=0))+
geom_errorbarh(aes(xmax=avgHeight+heightSE,xmin=avgHeight-heightSE))

# geom_errobarh creates the horizontal error bar for the variable in the X-axis (in this case avgHeight)

So, this is the final graph:
Horizontal & Vertical Error Bars

Please let me know if you have any questions or suggestions.