###################################################################################################################
# Title: r-scripts to fit density functions to tree size distributions 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$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)

###############################################################
#### bin the data by initial diameters (dbh)

bins <- 172 #number of bins specified in bp - use 1cm wide bins
#bp <- c(10^(seq(from = 0.475, to = 2.25, length.out = bins+1)))#log scale bins
bp <- c(seq(from = 3, to = 175, length.out = bins+1))#linear scale bins
data$Bin <- cut(data$Diameter.1, breaks=bp, include.lowest = TRUE)
mp <- c(rep(NA, bins))#calculate bin midpoints
for(j in 1:bins) { 
mp[j] <- (bp[j]+bp[j+1])/2
}
###############################################################
#### function to calculate density of trees in diameter size bin (number trees/area):

mvden <- function(tdata){
totarea <- 23.77952 #total area sampled in hectares
den <- dim(tdata)[1]/totarea
return(den)
}
###############################################################
#Calculate tree size distribution for all species combined and for common species individually:

#create empty data frames for results:
sizedata <- data.frame(matrix(NA, nrow = bins, ncol = 2+sppnum))
col.names <- c("Diameter", "Density", common.species)
names(sizedata) <- col.names

for(j in 1:bins) { 
cbin <- levels(data$Bin)[j]
bindat <- subset(data, data$Bin == cbin)
sizedata[j,"Diameter"] <- mp[j]
sizedata[j,"Density"] <- mvden(bindat)
for (i in 1:sppnum) {
spp <- common.species[i]
spdat <- subset(bindat , bindat$Species == spp)
sizedata[j,as.character(spp)] <- mvden(spdat)
}
}
####################################################################################################################
# define functions to model tree density as a function of diameter:

weibull.function <- function(diam,k,L,G)
{    (k*(diam)^(G-1))*(exp(-L*(diam)^G))  }

negexp.function <- function(diam,k,L)
{    k*exp(-L*diam)   }

power.function <- function(diam,k,L)
{    k*(diam^(-L))    }

#  PDFs

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

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

attach(sizedata)
working.data <- na.omit(data.frame(Density,Diameter))
detach(sizedata)

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

#  Set up the parameter list - k and L are defined as vectors of appropriate length
par <- list(k = 400, L = 0.05, sd = 1)

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

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(k = 0, L = 0, sd=0)
par_hi <- list(k = 1000, L = 1, sd=1000)

var$mean <- "predicted"
var$x <- "Density"

##  now call the annealing algorithm
results<-anneal(negexp.function,par,var,working.data,par_lo,par_hi,my.dnorm,"Density",max_iter=25000,hessian = TRUE)

## now write the results to a file...
write_results(results,"Density as a function of size modelled with negative exponential.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
xx<-working.data$Diameter
yy<-working.data$Density
x <- seq(0,max(working.data$Diameter, na.rm = TRUE),0.25)
y <- negexp.function(x,k=results$best_pars$k,L=results$best_pars$L)

windows()
plot(x,y, ylim=c(0,max(yy)), type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Population density (/ha)")
points(xx,yy,lwd=2,col="blue")
savePlot(file="Modelled density as a negative exponential function of size.png",type="png")

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

#  Set up the parameter list - k and L are defined as vectors of appropriate length
par <- list(k = 200, L = 0.8, sd = 1)

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

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

var$mean <- "predicted"
var$x <- "Density"

##  now call the annealing algorithm
results.2<-anneal(power.function,par,var,working.data,par_lo,par_hi,my.dnorm,"Density",max_iter=25000,hessian = TRUE)

## now write the results to a file...
write_results(results.2,"Density as a function of size modelled with negative power function.txt")

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

#  plot the predicted shape of the function
xx<-working.data$Diameter
yy<-working.data$Density
x <- seq(0,max(working.data$Diameter, na.rm = TRUE),0.25)
y <- power.function(x,k=results.2$best_pars$k,L=results.2$best_pars$L)

windows()
plot(x,y, ylim=c(0,max(yy)), type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Population density (/ha)")
points(xx,yy,lwd=2,col="blue")
savePlot(file="Modelled density as a function of size with negative power function.png",type="png")

##############  Set up the optimization using anneal for Weibull function ##############

#  Set up the parameter list - k, G and L are defined as vectors of appropriate length
par <- list(k = 95, L = 0.15, G = 0.7, sd = 5)

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

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(k = 0, L = 0, G=0, sd=0)
par_hi <- list(k = 100000, L = 1000, G=100, sd=100)

var$mean <- "predicted"
var$x <- "Density"

##  now call the annealing algorithm
results.3<-anneal(weibull.function,par,var,working.data,par_lo,par_hi,my.dnorm,"Density",max_iter=25000,hessian = TRUE)

## now write the results to a file...
write_results(results.3,"Density as a function of size modelled with Weibull function.txt")

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

#  plot the predicted shape of the function
xx<-working.data$Diameter
yy<-working.data$Density
x <- seq(0,max(working.data$Diameter, na.rm = TRUE),0.25)
y <- weibull.function(x,k=results.3$best_pars$k,L=results.3$best_pars$L,G=results.3$best_pars$G)

windows()
plot(x,y, ylim=c(0,max(yy)), type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Population density (/ha)")
points(xx,yy,lwd=2,col="blue")
savePlot(file="Modelled density as a function of size with Weibull function.png",type="png")

############################################################################################
############################################################################################
#  fit model distributions to each species' size data:

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

#Select data for each species and set up functions to fit to the data:

for (i in 1:sppnum) {

spp <- common.species[i]
attach(sizedata)
working.data <- na.omit(data.frame(Density=get(as.character(spp)),Diameter))
detach(sizedata)

##############  Set up the optimization using anneal for Weibull function ##############

#  Set up the parameter list - k, G and L are defined as vectors of appropriate length
par <- list(k = 95, L = 0.15, G = 0.7, sd = 5)

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

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(k = 0, L = 0, G=0, sd=0)
par_hi <- list(k = 1000000, L = 1000, G=100, sd=100)

var$mean <- "predicted"
var$x <- "Density"

##  now call the annealing algorithm
results.4<-anneal(weibull.function,par,var,working.data,par_lo,par_hi,my.dnorm,"Density",max_iter=25000,hessian = TRUE)

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

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

#  plot the predicted shape of the function
xx<-working.data$Diameter
yy<-working.data$Density
x <- seq(0,max(working.data$Diameter, na.rm = TRUE),0.25)
y <- weibull.function(x,k=results.4$best_pars$k,L=results.4$best_pars$L,G=results.4$best_pars$G)

windows()
plot(x,y, ylim=c(0,max(yy)), type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Population density (/ha)")
points(xx,yy,lwd=2,col="blue")
savePlot(file=paste("Weibull density function for ", as.character(spp), ".png",sep = ""),type="png")
dev.off()

##############  Set up the optimization using anneal for negative exponential ##############

#  Set up the parameter list - k and L are defined as vectors of appropriate length
par <- list(k = 400, L = 0.05, sd = 1)

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

#  Set up bounds and initial search steps for parameters to optimize
par_lo <- list(k = 0, L = 0, sd=0)
par_hi <- list(k = 1000, L = 1, sd=1000)

var$mean <- "predicted"
var$x <- "Density"

##  now call the annealing algorithm
results<-anneal(negexp.function,par,var,working.data,par_lo,par_hi,my.dnorm,"Density",max_iter=25000,hessian = TRUE)

## now write the results to a file...
write_results(results,paste("Negative exponential 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
xx<-working.data$Diameter
yy<-working.data$Density
x <- seq(0,max(working.data$Diameter, na.rm = TRUE),0.25)
y <- negexp.function(x,k=results$best_pars$k,L=results$best_pars$L)

windows()
plot(x,y, ylim=c(0,max(yy)), type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Population density (/ha)")
points(xx,yy,lwd=2,col="blue")
savePlot(file=paste("Negative exponential function for ", as.character(spp), ".png",sep = ""),type="png")
dev.off()

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

#  Set up the parameter list - k and L are defined as vectors of appropriate length
par <- list(k = 200, L = 0.8, sd = 1)

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

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

var$mean <- "predicted"
var$x <- "Density"

##  now call the annealing algorithm
results.2<-anneal(power.function,par,var,working.data,par_lo,par_hi,my.dnorm,"Density",max_iter=25000,hessian = TRUE)

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

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

#  plot the predicted shape of the function
xx<-working.data$Diameter
yy<-working.data$Density
x <- seq(0,max(working.data$Diameter, na.rm = TRUE),0.25)
y <- power.function(x,k=results.2$best_pars$k,L=results.2$best_pars$L)

windows()
plot(x,y, ylim=c(0,max(yy)), type="l",lwd=2,col="red",xlab="DBH (cm)",ylab="Population density (/ha)")
points(xx,yy,lwd=2,col="blue")
savePlot(file=paste("Negative power function for ", as.character(spp), ".png",sep = ""),type="png")
dev.off()

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