aaronberdanier

time <- 100
N <- Nhat <- dN <- K <- rep(NA,time)
P <- rlnorm(100,7,0.7)
r <- 0.05
Kmax <- 500
alpha <- mean(P)
N[1] <- 100
Nhat[1] <- rpois(1,N[1])
K[1] <- NA

aaronberdanier

# Calculating water balance based on algorithms in Willmott, Rowe, and Mintz (1985) J. Clim. 5:589-606
# Aaron Berdanier and Matt Kwit (2012)
# Questions or comments to: aaron.berdanier@duke.edu

# Packages:
library(reshape)

# Load the function:
waterbalance <- function(rawdat,nh,wstar){

aaronberdanier

# 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

aaronberdanier

# calculating water balance based on Willmott, Rowe, and Mintz 1985 J. Clim.

Temp <- read.csv(file.choose()) # monthly_temp.csv
Precip <- read.csv(file.choose()) # monthly_precip.csv
Temp <- Temp[Temp[,1]>1905,]
Precip <- Precip[Precip[,1]>1905,]

T <- t(Temp)[-1,]
colnames(T) <- Temp[,1]

aaronberdanier

t <- 10 # number of times
n <- 100 # number of individuals 'n'
x <- rnorm(t,0,1) # environmental variable 'x'
mu <- 4 # average survival response to 'x'
s <- 2 # variation among individuals in survival response to 'x'
b0 <- 1.5
beta <- rnorm(n,mu,s) # survival response for individual to 'x'

# Deterministic survival probabilty for each individual at each time
theta <- exp(b0+x%*%t(beta))/(1+exp(b0+x%*%t(beta))) 

aaronberdanier

response <- cbind(matrix(y,ncol=1),rep(1,n*t)-matrix(y,ncol=1)) # response matrix for GLM

# y ~ Bin(1,p)
# logit(p) = xB
model <- summary(glm(response~rep(x,n),binomial(link=logit)))

aaronberdanier

t <- 10 # number of times
n <- 100 # number of individuals 'n'
x <- rnorm(t,0,1) # environmental variable 'x'
mu <- 4 # average survival response to 'x'
s <- 2 # variation among individuals in survival response to 'x'
b0 <- 1.5
beta <- rnorm(n,mu,s) # survival response for individual to 'x'


# deterministic survival probabilty for each individual at each time

aaronberdanier

xdat <- seq(0,10,.1)
ydat <- sin(xdat)
ndat <- length(xdat)
err <- rnorm(ndat,1,0.25)
err_hi <- ydat + err
err_lo <- ydat - err

plot(ydat~xdat, type="n", ylim=c(min(err_lo),max(err_hi)))
# add the error polygon
# it draws around the shape, so one half must be read backwards

aaronberdanier

ng <- 1000

# create the progress bar
prog <- txtProgressBar(min=0, max=ng, char="*", style=3)

for(g in 1:ng){
 # do stuff
 z <- numeric(10000)
 z <- rgamma(10000, 1, 1)
 setTxtProgressBar(prog, g) # update the progress bar

aaronberdanier

td <- numeric(0)
ng <- 1000

### WITH A PROGRESS BAR
# create the progress bar
prog <- txtProgressBar(min=0, max=ng, char="*", style=3)
time1 <- Sys.time() # check system time
for(g in 1:ng){
 # do stuff
 z <- numeric(10000)