###############################################################
#
#         2011 Likelihood Course - Lab 8
#
#  Nonlinear Quantile Regression using nlrq in "quantreg"
#
#     adapted from the example in the nlrq help file
#
###############################################################

library(quantreg)

###########  Load the sapling growth data from Lab#2 on Day 1 ##############

BCdata <- read.table("BC Sapling Growth Data.txt",header=T,sep="\t",as.is=T)

# create a working dataset for one of the species
#  in one of the subzones

working.data <- subset(BCdata, spp == "HW" & subzone == "ICHMC2")

# drop any observations containing missing values for the independent or
#  dependent variables
#  NOTE:  by first "attaching" the dataframe, you don't have to preface
#    all of the variable names by the dataframe name...

attach(working.data)
data <- na.omit(data.frame(radius,global,meanrg))
detach(working.data)

#  create a variable with size (radius in mm) at the beginning
#    of the measurement period

data$size <- data$radius - 5*data$meanrg

# examine scatterplots

plot(data$global,data$meanrg)
windows()
plot(data$size,data$meanrg)


############  Define the scientific models you want to try ###############
#

# define a basic Michaelis Menton model

MM.model <- function(a,s,X1)
{ (a* X1)/((a/s)+X1) }

# a size-dependent Michaelis Menton model

size.MM.model <- function(a,s,d,X1,X2)
{  (a* X1*X2^d)/((a/s)+X1) }

#  add more as interested...


##############################################################################
#
# Now fit a series of quantile regression models for different quantiles (tau)
#
##############################################################################


#  SENSITIVITY TO STARTING VALUES
#  parameter estimates in nlrq appear to be very sensitive to the starting values
#  for the search
#  NOTE: for some reason, I get an error unless I use the start= option with a named
#    vector of parameters

#  First, try some default values of the parameters

initial <- list(a=1,s=1)

#  Now fit a model for tau = 0.5 (median)

res.q50 <- nlrq(meanrg ~ MM.model(a,s,global), data=data, start = initial, tau = 0.5, trace = T)

# display output
print(res.q50)

#  whoa!  even this simple model with fairly innocuous starting values gives utter nonsense for a
#     solution

#  I think you'll find that regardless of what quantile you are trying to fit, the model
#  converges on nonsense unless the starting values are fairly close to what they
#  should be

# start with more reasonable initial values
initial <- list(a=10,s=0.1)

#  Now fit a model for tau = 0.5 (median)

res.q50 <- nlrq(meanrg ~ MM.model(a,s,global), data=data, start = initial, tau = 0.5, trace = T)

# display output
#  NOTE:  print gives shorter output, but includes deviance
print(res.q50)
summary(res.q50)

#  this should make more sense.
#  Note:  the print output reports a "deviance": this is -2 * log likelihood
#    you'll notice that the reasonable output has much lower deviance than the first (garbage) output

#  BUT MORE IMPORTANTLY - you can use the reported deviance to compute AIC for the model, and 
#    then use model comparison to evaluate alternative models...
#
#  Remember:  AIC = -2*loglikelihood + 2*(number of parameters)

#  Plot the data and the fitted line
windows()
plot(data$global,data$meanrg)

#  nlrq has a "predict" method that takes the estimated parameters
#    and generates a vector of predicted values based on data passed
#    to the "predict" method - you can use this with the R "lines" function
#    to add the predicted line to your scatter plot
#    NOTE: that in "predict", I had to create a new (temporary) data list that
#      had a variable named "global" so that "predict" knew what to plug
#      into the model

lines(1:100, predict(res.q50,newdata=list(global=1:100)), col = 1)

# nlrq has a "fitted" method that stores the predicted values along the
#   quantile regression line, so you could always display the shape of the
#   line just using the "points" function, too

#  points(data$global,fitted(res.q50),col=2)

###########  Now, do other quartiles...   #################################

# 25% and 75% Quartile
res.q25 <- nlrq(meanrg ~ MM.model(a,s,global), data=data, start = initial, tau = 0.25, trace = T)

lines(1:100, predict(res.q25,newdata=list(global=1:100)), col = 2)
summary(res.q25)

res.q75 <- nlrq(meanrg ~ MM.model(a,s,global), data=data, start = initial, tau = 0.75, trace = T)

lines(1:100, predict(res.q75,newdata=list(global=1:100)), col = 2)
summary(res.q75)


# 10% and 90% Quantile
res.q10 <- nlrq(meanrg ~ MM.model(a,s,global), data=data, start = initial, tau = 0.10, trace = T)

lines(1:100, predict(res.q10,newdata=list(global=1:100)), col = 3)
summary(res.q10)

res.q90 <- nlrq(meanrg ~ MM.model(a,s,global), data=data, start = initial, tau = 0.90, trace = T)

lines(1:100, predict(res.q90,newdata=list(global=1:100)), col = 3)
summary(res.q90)

#  add a legend...
leg <- c("median","25% and 75% quartiles","10% and 90% quantiles")
legend(x="topleft", legend=leg, lty=1, col=1:3)


#########  Exercise  ##############

####  Try fitting a median model using a simple power function for growth (i.e.  growth = a * light ^s)

####  Is the model better than the Michaelis-Menton function?