Code for calculating psi_crit for B. gracilis in Colorado

BerdanierKwit_psicrit.R

# Estimates of psi_crit for Colorado were made with published data of water potential and leaf conductance
# fitted to an exponential-sigmoid model for the loss of conductance.
# The psi_crit parameter was estimated to match the criteria used by Craine et al. (2012),
# where psi_crit = water potential when conductance is 5% of maximum (or, 95% loss of conductance)
# Aaron Berdanier and Matt Kwit (2012)
# Questions or comments to: aaron.berdanier@duke.edu

# samples 1:14 are from Sala et al. (1981) Oecologia 28:327-331
# 15:17 are from Sala and Lauenroth (1982) Oecologia 53:301-304
# 18:19 are from Sala et al. (1982) J. Appl. Ecol 19:647-657
# Mid-day water potential (MPa)
psi <- c(-2.37,-2.92,-2.99,-2.65,-3.39,-3.48,-3.50,-3.32,-3.01,-3.87,-3.58,-3.78,-4.86,-5.00,-6.19,-5.17,-3.34,-3.00,-3.40)
# Leaf conductance (mm/s)
cond <- c(1.294,1.246,1.279,1.116,1.123,0.855,0.750,0.735,0.510,0.643,0.625,0.482,0.307,0.295,0.700,0.915,1.105,0.800,1.250)
# Maximum leaf conductance, from Sala et al. (1991)
condmax <- 1.6
# Percent loss of conductance
plc <- 1-cond/condmax

plot(psi,plc,xlim=c(min(psi),0),ylim=c(0,1))

# Exponential-sigmoid model
model <- function(b) 1/(1+exp(b[1]*(psi-b[2])))
y <- plc
# Likelihood function
likelihood <- function(par){
 -sum(dnorm(y,model(par[1:2]),sqrt(par[3]),log=T))
}
# Maximum likelihood parameter estimates
params <- nlminb(c(0.5,-4,0.05),likelihood,lower=c(0,min(psi),0.001),upper=c(5,max(psi),0.1))$par

# Parametric bootstrap to estimate psi_crit
nboot <- 1000
pboot <- matrix(NA,nboot,3)
for(i in 1:nboot){
 y <- rnorm(length(plc),plc,sqrt(params[3]))
 y[y>0.98] <- 0.98 # censor extreme values
 pars <- nlminb(params,likelihood,lower=c(0,min(psi),0.001),upper=c(5,max(psi),0.1))$par
 pboot[i,] <- pars
}
psicrits <- log(1/0.95-1)/pboot[,1]+pboot[,2] # bootstrapped psi_crit, the solution to the exponential-sigmoid model when plc=0.95
# To make psi_crit estimates more conservative,
# censor extremely negative values.
psicrits[psicrits<as.numeric(quantile(psicrits,0.05))] <- NA 

# Mean psi_crit
mu_psicrit <- mean(psicrits,na.rm=TRUE)
# Standard deviation psi_crit
sd_psicrit <- sd(psicrits,na.rm=TRUE)