###################################################################################################################
# Title: r-scripts to compute survival functions using a maximum likelihood approach 
# Notes: uses the likelihood package  
# Authors: Charles D. Canham (Cary Institute of Ecosystem Studies) and Christian O. Marks (The Nature Conservancy)
# Date: January 22, 2015
####################################################################################################################

library(likelihood)

# read in and format the data:
#setwd("/Users/cmarks/Documents/Data/demography")
data <- read.csv(file="CT_River.csv",header=T,as.is=T)
data$Date.1 <- as.Date(data$Date.1)
data$Date.2 <- as.Date(data$Date.2)
data$Species <- as.factor(data$Species)
data["remper"] <- as.numeric(data$Date.2 - data$Date.1)/365 # calc remeasurement period

####################################################################################################################
# define functions:

# model for suvival as function of Diameter.1
double.exponential.size.function <- function(diam,remper,a,b,c,d)
{   ((1-a*exp(b*diam/100)) * (exp(c*((diam/100)^d))))^remper }

# alternative model for survival as a constant (i.e. independent of diameter)
constant.function2 <- function(remper,a){(1-a)^remper}

# alternative model for threshold effect of Growing Degree Days (GDD)
step.function <- function(GDD,remper,a,b)
{   ifelse(GDD<a, 1, b^remper)}

# alternative models for survival as a function of size and flooding
flood.pow.function <- function(diam,flood,remper,a,b,c)
{   ((1-a*exp(b*diam/100)) * (flood/100)^c )^remper }

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

# pdf:
loglikelihood <- function(pred,observed)
{ ifelse(observed == 1, log(pred), log(1-pred)) }

####################################################################################################################
# fit the survival function to the data for all species combined:

data["Survive"] <- ifelse(data$mortality == FALSE, 1, 0)
attach(data)
working.data <- na.omit(data.frame(Survive,Diameter.1,remper))
detach(data)

##############  Set up the optimization using anneal  ##############

#  Set up the parameter list - a,b,c and d are defined as vectors of appropriate length
par <- list(a = 0.15, b = -5, c = -0.05, d = 4)

#  Set up the variable list
var <- list(diam = "Diameter.1", remper = "remper")

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(a = 0, b = -500, c = -100, d = 1)
par_hi <- list(a = 1, b = 0, c = 0, d = 100)

#  NOTE:  there are no parameters for the PDF (since the scientific model is probabilistic)
#   "predicted" is a reserved word in anneal, which takes the model and applies it to the
#   data and generates a new column in the data frame labeled "predicted"
var$pred <- "predicted"
var$observed <- "Survive"

##  now call the annealing algorithm
results<-anneal(double.exponential.size.function,par,var,working.data,par_lo,par_hi,loglikelihood,"Survive",max_iter=15000,hessian = TRUE)

## now write the results to a file...
write_results(results,"Survival as a function of size incl fall and broken.txt")

## display some of the results in the console
results$best_pars
results$max_likeli
results$aic_corr
results$slope
results$R2

## plot the fitted function
x <- seq(0,max(working.data$Diameter.1, na.rm = TRUE),0.25)
y <- double.exponential.size.function(x,1,results$best_pars$a,results$best_pars$b,results$best_pars$c,results$best_pars$d)
windows()
plot(x,y, type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Probability of Survival (/yr)")
savePlot(file="Predicted survival as a function of size including fallen and broken trees.png",type="png")

###############################################################################################################
######################## fit individual species survival functions ############################################

#take data for common species individually:
data$Species <- replace(data$Species, data$Species=="ACNISA", "ACSA3") #combine sugar maple with black maple hybrids as single species
sppN<-data.frame(N=summary(data$Species))
common.species <- c(row.names(subset(sppN,sppN$N>=50)))
sppnum<-length(common.species)

#setwd("/Users/cmarks/Documents/Data/demography/survival")

#Select data and set up functions to fit to the data:

for (i in 1:sppnum) {

spp <- common.species[i]
data2 <- subset(data, data$Species == spp)
attach(data2)
working.data <- na.omit(data.frame(Survive,Diameter.1,remper))
detach(data2)

##############  Set up the optimization using anneal  ##############

par <- list(a = 0.15, b = -5, c = -0.05, d = 4)

#  Set up the variable list
var <- list(diam = "Diameter.1", remper = "remper")

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(a = 0, b = -500, c = -100, d = 1)
par_hi <- list(a = 1, b = 0, c = 0, d = 100)

var$pred <- "predicted"
var$observed <- "Survive"

##  now call the annealing algorithm
results<-anneal(double.exponential.size.function,par,var,working.data,par_lo,par_hi,loglikelihood,"Survive",max_iter=15000,hessian = TRUE)

## now write the results to a file...
write_results(results,paste("Survival function for ", as.character(spp), ".txt", sep = ""))

## display some of the results in the console
results$best_pars
results$max_likeli
results$aic_corr
results$slope
results$R2

#  plot the predicted shape of the function
x <- seq(0,max(working.data$Diameter.1, na.rm = TRUE),0.25)
y <- double.exponential.size.function(x,1,results$best_pars$a,results$best_pars$b,results$best_pars$c,results$best_pars$d)
windows()
plot(x,y, type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Probability of Survival (/yr)")
savePlot(file=paste("Survival function for ", as.character(spp), ".png",sep = ""),type="png")
dev.off()

}
###############################################################################################################
######################### model beaver cutting for individual species  ########################################

#setwd("/Users/cmarks/Documents/Data/demography/mortspecies")
data["Survive"] <- ifelse(data$beavermort == FALSE, 1, 0)

#take data for common species individually:
data$Species <- replace(data$Species, data$Species=="ACNISA", "ACSA3") #combine sugar maple with black maple hybrids as single species
sppN<-data.frame(N=summary(data$Species))
common.species <- c(row.names(subset(sppN,sppN$N>=50)))
sppnum<-length(common.species)

#Select data and set up functions to fit to the data:

for (i in 1:sppnum) {

spp <- common.species[i]
data2 <- subset(data, data$Species == spp)
attach(data2)
working.data <- na.omit(data.frame(Survive,remper))
detach(data2)

##############  Set up the optimization using anneal  ##############

par <- list(a = 0.01)

#  Set up the variable list
var <- list(remper = "remper")

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(a = 0)
par_hi <- list(a = 1)

var$pred <- "predicted"
var$observed <- "Survive"

##  now call the annealing algorithm
results<-anneal(constant.function2,par,var,working.data,par_lo,par_hi,loglikelihood,"Survive",max_iter=15000,hessian = TRUE)

## now write the results to a file...
write_results(results,paste("Mortality caused by beaver cutting for ", as.character(spp), ".txt", sep = ""))

}

###########################################################################################################################
######################## model effect of GDD0C on bittersweet caused mortality ############################################

# select the data for the analysis
data["Survive"] <- ifelse(is.na(data$bitmort), NA, ifelse(data$bitmort == FALSE, 1, 0))
attach(data)
working.data <- na.omit(data.frame(Survive,GDD0C,remper))
detach(data)

##############  Set up the optimization using anneal  ##############

par <- list(a = 2500, b = 0.98)

#  Set up the variable list
var <- list(GDD = "GDD0C", remper = "remper")

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(a = 1, b = 0)
par_hi <- list(a = 10000, b = 1)

var$pred <- "predicted"
var$observed <- "Survive"

##  now call the annealing algorithm
results<-anneal(step.function,par,var,working.data,par_lo,par_hi,loglikelihood,"Survive",max_iter=10000,hessian = TRUE)

## now write the results to a file...
write_results(results,"Survival as step function of GDD.txt")

## display some of the results in the console
results$best_pars
results$max_likeli
results$aic_corr
results$slope
results$R2

#  plot the predicted shape of the function

x <- seq(min(working.data$GDD0C, na.rm = TRUE),max(working.data$GDD0C, na.rm = TRUE),0.25)
y <- step.function(x,1,results$best_pars$a,results$best_pars$b)
windows()
plot(x,y, type="l",lwd=2,col="red",xlab="GDD (base 0C)",ylab="Probability of Survival (/yr)")
savePlot(file="Survival as step function of GDD.png")

#######################################################################################################################
############## model beaver caused mortality as a function of size and flooding #######################################

#select the data for the analysis
#setwd("/Users/cmarks/Documents/Data/demography/env")
data["Survive"] <- ifelse(data$beavermort == FALSE, 1, 0)
attach(data)
data2 <- na.omit(data.frame(Survive,Diameter.1,ANNUALMIN,PH,exceedence,GDD0C,remper))
detach(data)
working.data <- data.frame(Survive=data2$Survive, Diameter.1=data2$Diameter.1, remper=data2$remper, exceedence=data2$exceedence)

##############  Set up the optimization using anneal for exponential function model ##############

par <- list(a = 0.2, b = -25, c = -0.0075)

#  Set up the variable list
var <- list(diam = "Diameter.1", flood = "exceedence", remper = "remper")

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(a = 0, b = -500, c = -100)
par_hi <- list(a = 1, b = 0, c = 0)

var$pred <- "predicted"
var$observed <- "Survive"

##  now call the annealing algorithm
results<-anneal(flood.exp.function,par,var,working.data,par_lo,par_hi,loglikelihood,"Survive",max_iter=10000,hessian = TRUE)

## now write the results to a file...
write_results(results,"Beaver Survival as exponential function of size and flooding.txt")

## display some of the results in the console
results$best_pars
results$max_likeli
results$aic_corr
results$slope
results$R2

##############  Set up the optimization using anneal for power function model ##############

par <- list(a = 0.2, b = -2.5, c = 0.01)

#  Set up the variable list
var <- list(diam = "Diameter.1", flood = "exceedence", remper = "remper")

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(a = 0, b = -500, c = 0)
par_hi <- list(a = 1, b = 0, c = 100)

var$pred <- "predicted"
var$observed <- "Survive"

##  now call the annealing algorithm
results<-anneal(flood.pow.function,par,var,working.data,par_lo,par_hi,loglikelihood,"Survive",max_iter=10000,hessian = TRUE)

## now write the results to a file...
write_results(results,"Beaver Survival as power function of size and flooding.txt")

## display some of the results in the console
results$best_pars
results$max_likeli
results$aic_corr
results$slope
results$R2

###################################################################################################