##################################################################################
#
#   2009 Likelihood Course - Lab 7 - Logistic Regression Example
#
#   Parameter estimation for Adirondack windthrow data, using simulated annealing
#
#   DATA FROM:
#   Canham, C. D., Papaik, M. J., and Latty, E. F.  2001.  Interspecific variation
#      in susceptibility to windthrow as a function of tree size and storm severity
#      for northern temperate tree species.
#      Canadian Journal of Forest Research 31:1-10.
#
##################################################################################

#  load the "likelihood" library... (version 1.0)
library(likelihood)

# load the windthrow data file
load("Damagedata.Rdata")

################################################################
#
# Functions for logistic regression of Adirondack windthrow data
#
################################################################

###############   Scientific Model function #########################
#
#  Simple nonlinear logistic regression for probability of windthrow as
#    a function of species, size, and storm severity.
#    {a,b,c} are species-specific parameters
#    See Canham et al. (2001, CJFR) for details

logistic <- function(a,b,c,S,DBH,Plot,Spp)
#  {a,b,c} are set up as vectors of length(Spp);  Plot is a vector of length(Plot)
{
  logit = exp(a[Spp] + c[Spp]*S[Plot]*(DBH^b[Spp]))
  logit/(1+logit)   # this returns the predicted probability of windthrow
}

################## Likelihood function ##########################

  # Since the scientific model is already expressed as a probability
  #   of observing a given (binary) result, the likelihood of
  #   observing a "yes" (windthrown) event is simply the predicted
  #   probability from the scientific model,
  #   and the likelihood of observing a "no" (not windthrown) is
  #   simply 1- predicted
  #   So, the log-likelihood function is simply:

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


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

#  Set up the parameter list - a,b, and c are defined as vectors of appropriate length
#  Note:  species and plots have been numbered sequentially starting at 1, so I
#            can use the maximum number as the number of species and plots...

par <- list(a = rep(0,max(Damagedata$Spp)),b = rep(1,max(Damagedata$Spp)),
               c = rep(0,max(Damagedata$Spp)), S = rep(0.5,max(Damagedata$Plot)))

#  Set up the variable list
var <- list(DBH = "DBH", Plot = "Plot", Spp = "Spp")

#  Set up bounds and initial search steps for parameters to optimize
par_min <- list(a = rep(-20,max(Damagedata$Spp)),b = rep(-4,max(Damagedata$Spp)),
               c = rep(-20,max(Damagedata$Spp)), S = rep(0,max(Damagedata$Plot)))
par_max <- list(a = rep(20,max(Damagedata$Spp)),b = rep(4,max(Damagedata$Spp)),
               c = rep(20,max(Damagedata$Spp)), S = rep(1,max(Damagedata$Plot)))

#  NOTE:  there are no parameters for the PDF (since the scientific model is probabilistic)
var$pred <- "predicted"
var$observed <- "Mortality"


##  now call the annealing algorithm
results<-anneal(logistic,par,var,Damagedata,par_min,par_max,loglikelihood,"Mortality",max_iter=10000,hessian = TRUE)

## now write the results to a file...
write_results(results,"sample windthrow logistic regression output.txt")

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