########################################################################################################
# Title: r-script to plot growth and mortality relationships for different environmental factors 
# Notes: includes ploting support intervals
# Author: Christian O. Marks (The Nature Conservancy)
# Date: January 20, 2015
########################################################################################################
library(gtools)

setwd("/Users/cmarks/Documents/Data/demography/env")
growth <- read.csv(file="EnvCoeffGrowth.csv",header=T)
mort <- read.csv(file="EnvCoeffMort.csv",header=T)
load("selfthinning.Rda")
load("equilibrium.Rda")
diam <- 30 # set the diameter for the comparison in cm

########################################################################################################
# prepare functions for computing support intervals:

#function to replace high (1) and low (0) from permutation with actual upper and lower support interval fitted values
valrep <- function(v,l,h){  
v <- ifelse(v==0, l, h)  
return(v) 
}

#function to get all permutations of upper and lower support interval values for fitted coefficients:
permhighlow <- function(r){ # r is the number of coefficients in equation to permute
bv <- permutations(n=2,r=r,v=c(0, 1),repeats.allowed=T)
return(bv)
}

# define demographic functions:

flood.mort.function <- function(diam,flood,a,b,c)
{   1 - ((1-a*exp(b*diam/100)) * exp(c*(flood/100)))}

gdd.mort.function <- function(diam,gdd,a,b,c)
{   1 - ((1-a*exp(b*diam/100)) * exp(c*(gdd/1000)))}

Tmin.mort.function <- function(diam,t,a,b,c)
{   1 - ((1-a*exp(b*diam/100)) * exp(c*(t+10)))}

pH.mort.function <- function(diam,pH,a,b,c)
{   1 - ((1-a*exp(b*diam/100)) * exp(c*(pH)))}

flood.grow.function <- function(diam,flood,a,b,c) { a*(diam^b)*exp(c*(flood/100)) }

gdd.grow.function <- function(diam,gdd,a,b,c) { a*(diam^b)*(gdd/1000)^c }

Tmin.grow.function <- function(diam,t,a,b,c) { a*(diam^b)*exp(c*(t+10)) }

pH.grow.function <- function(diam,pH,a,b,c) { a*(diam^b)*(pH)^c }

####################################################
# find fitted equation parameters for flooding:
sp1 <- "flood"
flood <- c(0,12,25,38,50)
a1 <- subset(growth$a,growth$Factor==sp1)
b1 <- subset(growth$b,growth$Factor==sp1)
c1 <- subset(growth$c,growth$Factor==sp1)
al1 <- subset(growth$al,growth$Factor==sp1)
bl1 <- subset(growth$bl,growth$Factor==sp1)
cl1 <- subset(growth$cl,growth$Factor==sp1)
au1 <- subset(growth$au,growth$Factor==sp1)
bu1 <- subset(growth$bu,growth$Factor==sp1)
cu1 <- subset(growth$cu,growth$Factor==sp1)
a <- subset(mort$a, mort$Factor==sp1)
b <- subset(mort$b, mort$Factor==sp1)
c <- subset(mort$c, mort$Factor==sp1)
al <- subset(mort$al, mort$Factor==sp1)
bl <- subset(mort$bl, mort$Factor==sp1)
cl <- subset(mort$cl, mort$Factor==sp1)
au <- subset(mort$au, mort$Factor==sp1)
bu <- subset(mort$bu, mort$Factor==sp1)
cu <- subset(mort$cu, mort$Factor==sp1)

len <- length(flood)
results.flood <- data.frame(factor=flood,growth=rep(NA,len),growth.l=rep(NA,len),growth.u=rep(NA,len),mort=rep(NA,len),mort.l=rep(NA,len),mort.u=rep(NA,len))

#compute observed flood mortality with support intervals: 
y1 <- flood.mort.function(diam,flood,a,b,c)
bv <- permhighlow(3)
coefs <- data.frame(a=valrep(bv[,1],al,au),b=valrep(bv[,2],bl,bu),c=valrep(bv[,3],cl,cu)) 
vals <- data.frame(matrix( rep(NA, len), nrow=len, ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- flood.mort.function(diam,flood,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl1 <- apply(vals,1,min)
yu1 <- apply(vals,1,max)
results.flood$mort <- y1
results.flood$mort.l <- yl1
results.flood$mort.u <- yu1

#compute observed flood growth with support intervals:
y0 <- flood.grow.function(diam,flood,a1,b1,c1) 
bv <- permhighlow(3)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),c1=valrep(bv[,3],cl1,cu1)) 
vals <- data.frame(matrix( rep(NA, len), nrow=len, ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- flood.grow.function(diam,flood,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results.flood$growth <- y0/10
results.flood$growth.l <- yl0/10
results.flood$growth.u <- yu0/10

# find fitted equation parameters for pH:
pH <- c(4.3,5,5.7,6.4,7.1)
sp1 <- "pH"
a1 <- subset(growth$a,growth$Factor==sp1)
b1 <- subset(growth$b,growth$Factor==sp1)
c1 <- subset(growth$c,growth$Factor==sp1)
al1 <- subset(growth$al,growth$Factor==sp1)
bl1 <- subset(growth$bl,growth$Factor==sp1)
cl1 <- subset(growth$cl,growth$Factor==sp1)
au1 <- subset(growth$au,growth$Factor==sp1)
bu1 <- subset(growth$bu,growth$Factor==sp1)
cu1 <- subset(growth$cu,growth$Factor==sp1)
a <- subset(mort$a, mort$Factor==sp1)
b <- subset(mort$b, mort$Factor==sp1)
c <- subset(mort$c, mort$Factor==sp1)
al <- subset(mort$al, mort$Factor==sp1)
bl <- subset(mort$bl, mort$Factor==sp1)
cl <- subset(mort$cl, mort$Factor==sp1)
au <- subset(mort$au, mort$Factor==sp1)
bu <- subset(mort$bu, mort$Factor==sp1)
cu <- subset(mort$cu, mort$Factor==sp1)

len <- length(pH)
results.pH <- data.frame(factor=pH,growth=rep(NA,len),growth.l=rep(NA,len),growth.u=rep(NA,len),mort=rep(NA,len),mort.l=rep(NA,len),mort.u=rep(NA,len))

#compute observed pH mortality with support intervals:
y1 <- pH.mort.function(diam,pH,a,b,c)
bv <- permhighlow(3)
coefs <- data.frame(a=valrep(bv[,1],al,au),b=valrep(bv[,2],bl,bu),c=valrep(bv[,3],cl,cu)) 
vals <- data.frame(matrix( rep(NA, len), nrow=len, ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- pH.mort.function(diam,pH,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl1 <- apply(vals,1,min)
yu1 <- apply(vals,1,max)
results.pH$mort <- y1
results.pH$mort.l <- yl1
results.pH$mort.u <- yu1

#compute observed pH growth with support intervals:
y0 <- pH.grow.function(diam,pH,a1,b1,c1) 
bv <- permhighlow(3)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),c1=valrep(bv[,3],cl1,cu1)) 
vals <- data.frame(matrix( rep(NA, len), nrow=len, ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- pH.grow.function(diam,pH,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results.pH$growth <- y0/10
results.pH$growth.l <- yl0/10
results.pH$growth.u <- yu0/10

# find fitted equation parameters for gdd:
gdd <- c(4200,3850,3500,3150,2800)
sp1 <- "gdd"
a1 <- subset(growth$a,growth$Factor==sp1)
b1 <- subset(growth$b,growth$Factor==sp1)
c1 <- subset(growth$c,growth$Factor==sp1)
al1 <- subset(growth$al,growth$Factor==sp1)
bl1 <- subset(growth$bl,growth$Factor==sp1)
cl1 <- subset(growth$cl,growth$Factor==sp1)
au1 <- subset(growth$au,growth$Factor==sp1)
bu1 <- subset(growth$bu,growth$Factor==sp1)
cu1 <- subset(growth$cu,growth$Factor==sp1)
a <- subset(mort$a, mort$Factor==sp1)
b <- subset(mort$b, mort$Factor==sp1)
c <- subset(mort$c, mort$Factor==sp1)
al <- subset(mort$al, mort$Factor==sp1)
bl <- subset(mort$bl, mort$Factor==sp1)
cl <- subset(mort$cl, mort$Factor==sp1)
au <- subset(mort$au, mort$Factor==sp1)
bu <- subset(mort$bu, mort$Factor==sp1)
cu <- subset(mort$cu, mort$Factor==sp1)

len <- length(gdd)
results.gdd <- data.frame(factor=gdd,growth=rep(NA,len),growth.l=rep(NA,len),growth.u=rep(NA,len),mort=rep(NA,len),mort.l=rep(NA,len),mort.u=rep(NA,len))

#compute observed GDD mortality with support intervals:
y1 <- gdd.mort.function(diam,gdd,a,b,c)
bv <- permhighlow(3)
coefs <- data.frame(a=valrep(bv[,1],al,au),b=valrep(bv[,2],bl,bu),c=valrep(bv[,3],cl,cu)) 
vals <- data.frame(matrix( rep(NA, len), nrow=len, ncol=nrow(bv)))  
for (i in 1:nrow(bv)){
vals[,i] <- gdd.mort.function(diam,gdd,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl1 <- apply(vals,1,min)
yu1 <- apply(vals,1,max)
results.gdd$mort <- y1
results.gdd$mort.l <- yl1
results.gdd$mort.u <- yu1

#compute observed GDD growth with support intervals:
y0 <- gdd.grow.function(diam,gdd,a1,b1,c1) 
bv <- permhighlow(3)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),c1=valrep(bv[,3],cl1,cu1)) 
vals <- data.frame(matrix( rep(NA, len), nrow=len, ncol=nrow(bv)))  
for (i in 1:nrow(bv)){
vals[,i] <- gdd.grow.function(diam,gdd,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results.gdd$growth <- y0/10
results.gdd$growth.l <- yl0/10
results.gdd$growth.u <- yu0/10

# find fitted equation parameters for Tmin:
Tmin <- c(2.9,1.05,-0.8,-2.65,-4.5)
sp1 <- "Tmin"
a1 <- subset(growth$a,growth$Factor==sp1)
b1 <- subset(growth$b,growth$Factor==sp1)
c1 <- subset(growth$c,growth$Factor==sp1)
al1 <- subset(growth$al,growth$Factor==sp1)
bl1 <- subset(growth$bl,growth$Factor==sp1)
cl1 <- subset(growth$cl,growth$Factor==sp1)
au1 <- subset(growth$au,growth$Factor==sp1)
bu1 <- subset(growth$bu,growth$Factor==sp1)
cu1 <- subset(growth$cu,growth$Factor==sp1)
a <- subset(mort$a, mort$Factor==sp1)
b <- subset(mort$b, mort$Factor==sp1)
c <- subset(mort$c, mort$Factor==sp1)
al <- subset(mort$al, mort$Factor==sp1)
bl <- subset(mort$bl, mort$Factor==sp1)
cl <- subset(mort$cl, mort$Factor==sp1)
au <- subset(mort$au, mort$Factor==sp1)
bu <- subset(mort$bu, mort$Factor==sp1)
cu <- subset(mort$cu, mort$Factor==sp1)

len <- length(Tmin)
results.Tmin <- data.frame(factor=Tmin,growth=rep(NA,len),growth.l=rep(NA,len),growth.u=rep(NA,len),mort=rep(NA,len),mort.l=rep(NA,len),mort.u=rep(NA,len))

#compute observed Tmin mortality with support intervals:
y1 <- Tmin.mort.function(diam,t=Tmin,a,b,c)
bv <- permhighlow(3)
coefs <- data.frame(a=valrep(bv[,1],al,au),b=valrep(bv[,2],bl,bu),c=valrep(bv[,3],cl,cu)) 
vals <- data.frame(matrix( rep(NA, len), nrow=len, ncol=nrow(bv)))  
for (i in 1:nrow(bv)){
vals[,i] <- Tmin.mort.function(diam,t=Tmin,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl1 <- apply(vals,1,min)
yu1 <- apply(vals,1,max)
results.Tmin$mort <- y1
results.Tmin$mort.l <- yl1
results.Tmin$mort.u <- yu1

#compute observed Tmin growth with support intervals:
y0 <- Tmin.grow.function(diam,t=Tmin,a1,b1,c1) 
bv <- permhighlow(3)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),c1=valrep(bv[,3],cl1,cu1)) 
vals <- data.frame(matrix( rep(NA, len), nrow=len, ncol=nrow(bv)))  
for (i in 1:nrow(bv)){
vals[,i] <- Tmin.grow.function(diam,t=Tmin,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results.Tmin$growth <- y0/10
results.Tmin$growth.l <- yl0/10
results.Tmin$growth.u <- yu0/10

############## create plots of mort versus growth at different env values ##################

#set up plotting area:
par(mfrow=c(2,2))
par(mar = c(0, 0, 0, 0), oma = c(4, 4, 0.5, 0.5))
ming <- 0.25
maxg <- 0.75

######### make plot for pH ###############

plot(results.pH$growth, results.pH$mort, xlim=c(ming,maxg), ylim=c(0,0.1), xaxt="n", cex.axis=1.0,cex.lab=1.25, pch=16)
arrows(results.pH[,'growth'], results.pH[,'mort.l'], results.pH[,'growth'], results.pH[,'mort.u'], col="grey", length=0.02, angle=90, code=3)
arrows(results.pH[,'growth.l'], results.pH[,'mort'], results.pH[,'growth.u'], results.pH[,'mort'],  col="grey", length=0.02, angle=90, code=3)
text(results.pH$growth, results.pH$mort.u, paste("pH =",as.character(round(results.pH$factor, digits = 2))), cex=0.8, pos=3, col=1) 
points(results.pH[,'growth'], results.pH[,'mort'], pch=16)
lines(equilthin$x,equilthin$y,lwd=2, lty=3)
lines(selfthin$x,selfthin$y,lwd=2, lty=3)

######### make plot for GDD0C ###############

plot(results.gdd$growth, results.gdd$mort, xlim=c(ming,maxg), ylim=c(0,0.1), xaxt="n",yaxt="n", cex.axis=1.0,cex.lab=1.25, pch=16)
arrows(results.gdd[,'growth'], results.gdd[,'mort.l'], results.gdd[,'growth'], results.gdd[,'mort.u'], col="grey", length=0.02, angle=90, code=3)
arrows(results.gdd[,'growth.l'], results.gdd[,'mort'], results.gdd[,'growth.u'], results.gdd[,'mort'],  col="grey", length=0.02, angle=90, code=3)
text(results.gdd$growth.u, results.gdd$mort, paste("GDD =",as.character(round(results.gdd$factor, digits = 2))), cex=0.8, pos=4, col=1) 
points(results.gdd[,'growth'], results.gdd[,'mort'], pch=16)
lines(equilthin$x,equilthin$y,lwd=2, lty=3)
lines(selfthin$x,selfthin$y,lwd=2, lty=3)

######### make plot for flood exceedence ###############

plot(results.flood$growth, results.flood$mort, xlim=c(ming,maxg), ylim=c(0,0.1), cex.axis=1.0,cex.lab=1.25, pch=16)
arrows(results.flood[,'growth'], results.flood[,'mort.l'], results.flood[,'growth'], results.flood[,'mort.u'], col="grey", length=0.02, angle=90, code=3)
arrows(results.flood[,'growth.l'], results.flood[,'mort'], results.flood[,'growth.u'], results.flood[,'mort'],  col="grey", length=0.02, angle=90, code=3)
text(results.flood$growth.u, results.flood$mort, paste("Flood =",as.character(round(results.flood$factor, digits = 2))), cex=0.8, pos=4, col=1)
points(results.flood[,'growth'], results.flood[,'mort'], pch=16)
lines(equilthin$x,equilthin$y,lwd=2, lty=3)
lines(selfthin$x,selfthin$y,lwd=2, lty=3)

######### make plot for minimum daily temperature averaged over the year ###############

plot(results.Tmin$growth, results.Tmin$mort, xlim=c(ming,maxg), ylim=c(0,0.1), yaxt="n", cex.axis=1.0,cex.lab=1.25, pch=16)
arrows(results.Tmin[,'growth'], results.Tmin[,'mort.l'], results.Tmin[,'growth'], results.Tmin[,'mort.u'], col="grey", length=0.02, angle=90, code=3)
arrows(results.Tmin[,'growth.l'], results.Tmin[,'mort'], results.Tmin[,'growth.u'], results.Tmin[,'mort'],  col="grey", length=0.02, angle=90, code=3)
text(results.Tmin$growth.u, results.Tmin$mort, paste("Tmin =",as.character(round(results.Tmin$factor, digits = 2))), cex=0.8, pos=4, col=1) 
points(results.Tmin[,'growth'], results.Tmin[,'mort'], pch=16)
lines(equilthin$x,equilthin$y,lwd=2, lty=3)
lines(selfthin$x,selfthin$y,lwd=2, lty=3)

mtext("DBH growth rate (cm/year)", side = 1, outer = TRUE, cex = 1.25, line = 2.3,col = "black")
mtext("Probability of mortality (/year)", side = 2, outer = TRUE, cex = 1.25, line = 2.3,col = "black")

savePlot(file="Effect of environment on demographic rates for all species combined.png",type="png")