Ecology with R
From BiodivBorneo09
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)
