Ecology with R

From BiodivBorneo09

Jump to: navigation, search

Download the data file here.

Contents

Analysis of tree plot data

Enter our trees, look at growth rates; OBSERVATIONAL

# Load our data objects
load("Rday2.RData")
ls()

# The 6-year old data for our 20x20 subplot
lamplot[1:10, ]

# A test set is in the set
# tmp <- lamplot[100:120,c(-3,-4)]
# myplot <- data.frame(tmp$tag, as.numeric(tmp$diam)+1, tmp$taxon)
# colnames(myplot) <- c("tag", "diam", "taxon")

# NB import your data as ourplot, e.g.
ourplot <- read.table("ourplot.csv", sep=",", header=T)
# use "tag", "diam", "taxon" as the headings

# we need to merge with the 6-year-old data
m <- merge(lamplot, myplot, by ="tag")
# or m <- merge(lamplot, myplot, by.x ="tag", by.y="MyTagHeader") if diff

# now get a difference:
ddiam <- m$diam.y - m$diam.x
ddiam
hist(ddiam)
# the outliers and negative growth rates!
mean(ddiam , na.rm=T)

# now relative growth rate:
rgr <- (m$diam.y - m$diam.x) / m$diam.x

Analysis of litter arthropod data, counts

# Our litter data set in a single file
litter

# let's split it into taxa and environment:
dim(litter)
littax <- litter[,8:30]
rownames(littax) <- litter$sample
litenv <- litter[ , c(1,2,4,5,6,7) ]
rownames(litenv) <- litter$sample

# total indivs per sample
rowSums(littax)

# differnce in sampling method:
plot(rowSums(littax) ~ litenv$type)
t.test(rowSums(littax) ~ litenv$type)
wilcox.test(rowSums(littax) ~ litenv$type)

# but these are counts of individuals in space, so use an appropriate model
lit.glm <- glm(rowSums(littax) ~ litenv$type, family=poisson)
summary(lit.glm)

# two factors
lit.glm <- glm(rowSums(littax) ~ litenv$type + litenv$gap, family=poisson)
summary(lit.glm)

# but hang on, let's look first. Litter traps:
plot(rowSums(littax)[1:6] ~ litenv$gap[1:6])

# an effect in one direction:
lit.glm <- glm(rowSums(littax)[1:6] ~ litenv$gap[1:6], family=poisson)
summary(lit.glm)

# and another for the pitfalls
plot(rowSums(littax)[7:12] ~ litenv$gap[7:12])
lit.glm <- glm(rowSums(littax)[7:12] ~ litenv$gap[7:12] , family=poisson)
summary(lit.glm)

# so we see a strong interaction effect:    
lit.glm <- glm(rowSums(littax) ~ litenv$type + litenv$gap +  litenv$type:litenv$gap, family=poisson)
summary(lit.glm)

Other statistical/ecological tests

A simple classification of tests, based on input data class:

independent variable: Continuous Nominal
(Y/N, live/die...)
dependent variable:
Continuous scatterplots, cor.test, lm, glm box plots, t.tests,
wilcox.test, lm, glm
Nominal logistic regression (glm, family=binomial) contingency tables, chisq.test, fisher.test


  • Nominal vs. nominal:
comp
table(comp)
chisq.test(comp$sp1, comp$sp2)
fisher.test(comp$sp1, comp$sp2)
  • Continuous as independent variable, nominal as dependent:
mort
plot( mort$surv ~ mort$turb)
summary( glm( mort$surv ~ mort$turb, family=binomial))

Community ecology

# Let's look at species.  Collapse sampling methods and convert to presence/absence
litsp <- littax[1:6, ] + littax[7:12,]
litsp > 0
litsp[litsp > 0]
litsp[litsp > 0] <- 1
rowSums(litsp)
wilcox.test( rowSums(litsp) ~ litenv$gap[1:6] )

The vegan package

# Now some compositional analysis
library(vegan)

# demo the Rank Abundance curve
data(BCI)
BCI[1:10,1:3]
plot(rad.lognormal(BCI[1,]))

# diversity indices
diversity(BCI, "shannon")
barplot(sort(diversity(BCI, "shannon")))

# demo clustering
data(dune)
data(dune.env)

plot(hclust(vegdist(dune) method="ward"))
plot(hclust(vegdist(dune) method="ward"), labels=dune.env$Management)

Analysis of litter arthropod data, species composition

# For our data
plot( rad.lognormal( colSums(littax[1:6,]) ) )
plot( rad.lognormal( colSums(littax[7:12,]) ) )

# make the distance matrix
vegdist(littax, method="bray")

# clustergram of taxon composition
plot( hclust( vegdist(littax, method="bray") ) )
plot( hclust( vegdist(littax, method="bray") ), label = litenv$type )
plot( hclust( vegdist(littax, method="bray"), method="ward" ) )
# dissatisfying difference in results?  But remember the cluster is unrooted.

# Permution analysis of similarity
summary( anosim (vegdist(littax, method="bray"), litenv$type) )
plot( anosim (vegdist(littax, method="bray"), litenv$type) )
summary( anosim (vegdist(littax, method="bray"), litenv$gap) )
plot( anosim (vegdist(littax, method="bray"), litenv$gap) )

Miscellaneous tools

Compare two regression slopes with:

novar.test(case1.lm, case2.lm)

Save your data to a binary file:

save(lamplot, myplot, litter, file="Rday2.R")

Fix your data:

fix(comp)