###################################################################################################################
# Title: r-scripts to compute size-dependent growth 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
data["growth.mm"] <- data$Growthrate*10 

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

# model for suvival as function of Diameter.1
power.function <- function(diam,b,c) { b*(diam^c) }

# alternative models for growth as function of Diameter.1 and soil pH
pH.pow.function <- function(diam,a,b,c,pH) { a*(diam^b)*(pH)^c }

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

#  PDFs
my.dnorm <- function(x,mean,sd) {dnorm(x,mean,sd,log=T)}

# define an alternative to the normal PDF in which the variance is a power function of the mean
normal.with.power.variance <- function(x,mean,sigma)
{  sd <- mean^sigma
   dnorm(x,mean,sd,log=T)
}
########################################################################################################################
##############  Model growth rate as a power function of tree diameter #################################################

attach(data)
working.data <- na.omit(data.frame(growth.mm,Diameter.1,ANNUALMIN,PH,exceedence,GDD0C))
detach(data)

##############  Set up the optimization using anneal  ##############
#  Set up the parameter list - b and c are defined as vectors of appropriate length
par <- list(b = 1, c = 1, sd = 1)

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

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

var$mean <- "predicted"
var$x <- "growth.mm"

##  now call the annealing algorithm
results1<-anneal(power.function,par,var,working.data,par_lo,par_hi,my.dnorm,"growth.mm",max_iter=10000,hessian = TRUE)

## now write the results to a file...
write_results(results1,"Growth as a function of size - no intercept.txt")

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

#  plot the predicted shape of the function
x <- seq(0,max(working.data$Diameter.1, na.rm = TRUE),0.25)
y <- power.function(x,results1$best_pars$b,results1$best_pars$c)
windows()
plot(x,y, ylim=c(0,max(y)),type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Growth Rate (mm/yr)")
savePlot(file="Predicted growth as a function of size - without intercept.png",type="png")

###############################################################################################################################
################# model growth as a function of diameter and pH ###############################################################

attach(data)
working.data <- na.omit(data.frame(growth.mm,Diameter.1,ANNUALMIN,PH,exceedence,GDD0C))
detach(data)
# Note that to compare alternative models with AIC values one needs to use the exact same date to fit both models

par <- list(a = 1, b = 1, c = 0.5, sd = 1)

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

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

var$mean <- "predicted"
var$x <- "growth.mm"

##  now call the annealing algorithm
results1<-anneal(pH.pow.function,par,var,working.data,par_lo,par_hi,my.dnorm,"growth.mm",max_iter=10000,hessian = TRUE)

## now write the results to a file...
write_results(results1,"Growth as power function of size and pH.txt")

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

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

par <- list(a = 1, b = 1, c = 0.5, sd = 1)

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

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

var$mean <- "predicted"
var$x <- "growth.mm"

##  now call the annealing algorithm
results1<-anneal(pH.exp.function,par,var,working.data,par_lo,par_hi,my.dnorm,"growth.mm",max_iter=10000,hessian = TRUE)

## now write the results to a file...
write_results(results1,"Growth as exponential function of size and pH.txt")

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

########################################################################################################################
######################## individual species size-dependent growth functions ############################################

#setwd("/Users/cmarks/Documents/Data/demography/growth")
data$Species <- as.factor(data$Species)
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)

for (i in 1:sppnum) {

#take data for common species individually:
spp <- common.species[i]
data2 <- subset(data, data$Species == spp)
attach(data2)
working.data <- na.omit(data.frame(growth.mm,Diameter.1))
detach(data2)

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

#  Set up the parameter list - a,b, and c are defined as vectors of appropriate length
par <- list(b = 1, c = 1, sd = 1)

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

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

var$mean <- "predicted"
var$x <- "growth.mm"

##  now call the annealing algorithm
results1<-anneal(power.function,par,var,working.data,par_lo,par_hi,my.dnorm,"growth.mm",max_iter=10000,hessian = TRUE)

## now write the results to a file...
write_results(results1,paste("Growth as a function of size without intercept for ", as.character(spp), ".txt", sep = ""))

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

#  plot the predicted shape of the function

x <- seq(0,max(working.data$Diameter.1, na.rm = TRUE),0.25)
y <- power.function(x,results1$best_pars$b,results1$best_pars$c)
windows()
plot(x,y, ylim=c(0,max(y)),type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Growth Rate (mm/yr)")
savePlot(file=paste("Predicted growth as a function of size without intercept for ", as.character(spp), ".png",sep = ""),type="png")
dev.off()

}
###############################################################################################################################