This example is intended to present the methods develloped for the paper Blarquez *et al.* (2013). We will use the macroremain record of the Lac du Loup a small subalpine lake with continuous macoremain counts for the last 11 750 years. More details on this site can be found in Blarquez et al. (2010).

First load the macroremain data:

dataMacro=read.csv("http://blarquez.com/public/code/dataMacro.csv")
# Macroremain counts
paramMacro=read.csv("http://blarquez.com/public/code/paramMacro.csv")
# Associated depths, ages and sample volumes
dataMacro[is.na(dataMacro)]=0
# Set missing values to zero

We now caluculates the median resolution of the record and the interpolated ages at which the influx will be calculated:

resMed=median(diff(paramMacro[,3]))
Ages=seq(-50, 11750, resMed)

We reconstruct the influx matrix using the pretreatment function which is the R implementation of the CharAnalysis CharPretreatment.m function originally develloped by P. Higuera and available at https://sites.google.com/site/charanalysis

Requires a matrix as input with the following columns:

CmTop CmBot AgeTop AgeBot Volume (params)

A serie from which to calculate accumulation rates (serie)

source("http://blarquez.com/public/code/pretreatment.R")
infMacro=matrix(ncol=length(dataMacro[1,]),nrow=length(Ages))
for (i in 1:length(dataMacro[1,])){
+ infMacro[,i]=c(pretreatment(params=paramMacro,serie=dataMacro[,i],
+ yrInterp=resMed,first=-50,last=11750,
+ plotit=F,Int=F)$accI)
+ }

Calculate samples macroremain sums and the minimun influx sum:

S=rowSums(infMacro)
nMin=min(na.omit(S[1:length(S)-1]))

We exclude the last sample (not calculated because of unknown accumulation rate) and samples with NA (because of vol==0):

del=which(S==0 | is.na(S))
infMacro=infMacro[-del,]
Ages=Ages[-del]
S=S[-del]

We calculate a proportion matrix which is used to replicate the rarefaction procedure:

S_=matrix(ncol=ncol(infMacro),nrow=nrow(infMacro))
for (i in 1:length(infMacro[1,])){
+ S_[,i]=c(S)
+ }
propMat=infMacro/S_

Replicate rarefaction with min influx n=1,2,…,500.

We will use the rarefy function from the package vegan:

library(vegan)
rare_=matrix(nrow=nrow(infMacro),ncol=500)
for (n in 1:500){
+ infN=ceiling(S_*n/nMin*propMat)
+ rare_[,n]=c(rarefy(infN,sample=n))
+ }

And finnally calculates the 90 percent confidence intervals and plot the rarefied richness using a 500 years locally weighted scatter plot smoother:

CI_=t(apply(rare_, 1, quantile, probs = c(0.05, 0.5, .95), na.rm = TRUE))
span=500/resMed/length(Ages)
plot(Ages,lowess(CI_[,2], f =span)$y,type="l",ylab="E[Tn]",ylim=c(min(CI_),
+ max(CI_)),xlab="Age",main="Rarefied Richness",xlim=c(12000,-100))
lines(Ages,lowess(CI_[,1], f =span)$y,type="l",lty=2)
lines(Ages,lowess(CI_[,3], f =span)$y,type="l",lty=2)

References:

Blarquez O., Finsinger W. and C. Carcaillet. 2013. Assessing paleo-biodiversity using low proxy influx. PLoS ONE 8(4): e65852

Blarquez O., Carcaillet C., Mourier B., Bremond L. and O. Radakovitch. 2010. Trees in the subalpine belt since 11 700 cal. BP : origin, expansion and alteration of the modern forest. The Holocene 20 : 139-146. DOI : 10.1177/0959683609348857

Download pdf vignette