Here is the bootstrapping code to calculate confidence intervals for the variance components, as done in Messier et al. 2010. Ecology Letters (DOI: 10.1111/j.1461-0248.2010.01476). I have received numerous requests to share, so have decided to make it available publicly.

*** Note that***:

– The nesting structure of the model fits the sampling design of this study. Different studies using different sampling designs will have to adjust the structure of the model accordingly.

– This code is provided is as a courtesy. It was written in 2009, with the versions of the R packages available at the time. I cannot guarantee that it will still work as the packages are updated.

============================================

d<-read.csv(‘C:/Documents and Settings/U/Desktop/MSc/Data/Mydata/lft.csv’,h=T,na.strings=”NA”)

#Adjust path as necessary

library(ape)

library(nlme)

# VARCOMP

varcomp(lme(log(LMA)~1,random=~1|Site/Plot/Species/Tree/Strata,data=d,na.action=na.omit),1) # This is a fully nested model. Used in Messier et al. 2010. Ecology Letters.

# CONFIDENCE INTERVALS

# LMA –

Site<-rep(NA,500) ## create empty vectors in which to store data

Plot<-rep(NA,500)

Species<-rep(NA,500)

Tree<-rep(NA,500)

Strata<-rep(NA,500)

Leaf<-rep(NA,500)

for (i in 1:500){ ## Do 500 bootstrap runs

d1<-d[sample(1:1952,1952,replace=T),] ## Randomize dataset

try(varcomp.LMA<-varcomp(lme(log(LMA)~1,random=~1|Site/Plot/Species/Tree/Strata,data=d1,na.action=na.omit),1),silent=T)

## store output of each bootstrap runs in these vectors

if (is.numeric(varcomp.LMA)){

Site[i]<-varcomp.LMA[1]

Plot[i]<-varcomp.LMA[2]

Species[i]<-varcomp.LMA[3]

Tree[i]<-varcomp.LMA[4]

Strata[i]<-varcomp.LMA[5]

Leaf[i]<-varcomp.LMA[6]

}

}

## use ‘try’ and ‘if (is.numeric)’ because the bootstrap runs that fail produce an error

## message, which cannot be stored in the vectors. These bits of code avoid the whole

## analysis aborting if one of the runs fails.

(CI.Site<-quantile(Site,c(0.025,0.975),na.rm=T))

(CI.Plot<-quantile(Plot,c(0.025,0.975),na.rm=T))

(CI.Species<-quantile(Species,c(0.025,0.975),na.rm=T))

(CI.Tree<-quantile(Tree,c(0.025,0.975),na.rm=T))

(CI.Strata<-quantile(Strata,c(0.025,0.975),na.rm=T))

(CI.Leaf<-quantile(Leaf,c(0.025,0.975),na.rm=T))

### Use na.rm=T because the bootstrap runs that fail will print NAs in the vector, which ### will affect the calculation of the quantiles.

====================================================

### Like this:

Like Loading...