################################################################################
# Ex4.3                                                                        #
#
# Last revised: Spring 2019
#
# Copyright © 2018 by Rodrigo Labouriau
################################################################################

load("DataFramesStatisticalModelling.RData")

attach(Ex4.3)

# -----------------------------------------------------------------------------#
# Exploratory analysis                                                         #
# -----------------------------------------------------------------------------#

str(Ex4.3)

table(Temperature)

# Studying the dependency of the mean on the temperature
plot(Temperature, Ninsects)
mmean <- tapply(Ninsects, factor(Temperature),mean)
points(levels(factor(Temperature)), mmean, type='b', pch=19, cex=1.5)

plot(Temperature, log(Ninsects))
points(levels(factor(Temperature)), log(mmean), pch=19, cex=1.5)
abline(coef(lm(log(Ninsects)~ Temperature), col='red'))

# Studying the relationship between the mean and the variance
vvar <- tapply(Ninsects, factor(Temperature),var)
plot(mmean, vvar)
abline(lm(vvar~ mmean ))
summary(lm(vvar~ mmean))

# -----------------------------------------------------------------------------#
# Fiting several Poisson models                                                #
# -----------------------------------------------------------------------------#

# Free-curve model

free.curve <- glm(Ninsects~factor(Temperature), family=poisson(link='log'))
summary(free.curve)
deviance(free.curve)
pchisq(deviance(free.curve), df=length(Ninsects)-5, lower.tail=F)

scatter.smooth (fitted.values(free.curve),residuals(free.curve,'response'),
     ylab='Raw residuals',
     xlab='Fitted values',
     main='Free curve'
)
abline(h=0)

plot(fitted.values(free.curve),residuals(free.curve,'pearson'),
     ylab='Pearson residuals',
     xlab='Fitted values',
     main='Free curve'
)
abline(h=0)

# Exponential response to the temperature

exponential.curve <- glm(Ninsects~Temperature, family=poisson(link='log'))
summary(exponential.curve)

plot(fitted.values(exponential.curve),residuals(exponential.curve,'response'),
     ylab='Raw residuals',
     xlab='Fitted values',
     main='Exponential curve'
)
abline(h=0)

plot(fitted.values(exponential.curve),residuals(exponential.curve,'pearson'),
     ylab='Pearson residuals',
     xlab='Fitted values',
     main='Exponential curve'
)

anova(exponential.curve,free.curve, test='Chisq')

deviance(exponential.curve)
pchisq(deviance(exponential.curve), df=length(Ninsects)-2, lower.tail=F)

# Null model

null.model <-  glm(Ninsects~ 1, family=poisson(link='log'))
summary(null.model)
anova(null.model, exponential.curve , test='Chisq')

scatter.smooth(Temperature,residuals(null.model,'response'),
     ylab='Raw residuals',
     xlab='Temperature',
     main='Null-model'
)
abline(h=0)


# Non-scaled exponential

non.scaled <-  glm(Ninsects~ Temperature +0, family=poisson(link='log'))
summary(non.scaled)
anova(non.scaled, exponential.curve , test='Chisq')

################################################################################
#                                                                              #
#  Alternative solution to the exercise Ex.4.3                                 #
#                                                                              #
################################################################################

# Here we give an alternative model to the temperature dependence.
# The model is inspired in a classic thermodynamic model, the Arrhenius model.
#
# According to this model the rate of a chemical reaction is given by the
# Arrhenius equation:  
# k = A exp(-Ea/ R T),
# where, k is the rate of the reaction, Ea is the activation energy, R is
# the universal gas constant and T is the absolute temperature.
# Taking logarithms, we see that the Arrhenius equation is equivalent to
# log(k) = log(A) - Ea/R 1/T,
# that is, the logarithm of the rate is an affine (i.e. linear) transformation 
# of the inverse (i.e. one divided by) of the temperature.

# First, I check whether such a model would give a reasonable approximation
# for our data by drawing the so called Arrhenius plot. If such a plot show
# a linear pattern then the model describes well the data. 

# Calculating the negative of the inverse of the temperatures
Atemp <- -1/(Temperature + 273.15)

# Drawing the Arrhenius plot
plot(Atemp, log(Ninsects), xlab="Negative inverse of the absolute temperature",
     main="Arrhenius Plot")
abline(coef(lm(log(Ninsects)~Atemp)))

Arrhenius.model <- glm(Ninsects ~ Atemp, family=poisson(link="log"))

free.curve <- glm(Ninsects~factor(Atemp), family=poisson(link='log'))

# Testing homogeneity 
deviance (free.curve)
pchisq(deviance (free.curve), df=95, lower.tail=F)

# So, we do not have evidences of lack of homogeneity

# Testing linearity of the Arrhenius plot, i.e. adecquacy of the Arrhenius model
anova(Arrhenius.model, free.curve, test="Chisq")

# So, there is no evidence against the Arrhenius model!

plot(Atemp, residuals(Arrhenius.model, method="pearson"))
abline(h=0)

coef(Arrhenius.model)

detach(Ex4.3)