################################################################################
#                                                                              #
# Solution to the exercises 3.3                                                #
#                                                                              #
# Last revised: Spring 2019
#
# Copyright © 2018 by Rodrigo Labouriau                              #
################################################################################

load("DataFramesStatisticalModelling.RData")

attach(Ex3.3)
str(Ex3.3)

# Exploratory analyses
table(P,field)

Proportions <- tapply(mycor, P, sum)/tapply(nsamples, P, sum)
Proportions

# Ploting the proportions against the P
plot(c(0,100,200,300,400,500,600), Proportions, type="b")

logit <- function(p) log(p/(1-p))

logit(Proportions)

# Ploting the lodds against the P
plot(c(0,100,200,300,400,500,600), logit(Proportions), type="b")
# Drawing a line (roughtly) adjusted wilth least square (this is primitive!)
LSM <- lm(logit(Proportions)~c(0,100,200,300,400,500,600))
coef(LSM)
abline(coef(LSM), lty=2)

# Preparing the response matrix for fitting a series of binomial models
Resp <- cbind(mycor, nsamples-mycor)
Resp[1:10, ]

FreeCurve <- glm(Resp~factor(P), family=binomial(link="logit"))
summary(FreeCurve)

deviance(FreeCurve)
length(coef(FreeCurve )) # Number of parameters
length(P)                # Number of observations
pchisq(deviance(FreeCurve), df=140-7, lower.tail = FALSE) # Test of homogeneity

LogisticModel <- glm(Resp~P, family=binomial(link="logit"))
summary(LogisticModel)

# Testing whether the logistic regression curve is adequate
anova(LogisticModel, FreeCurve, test="Chisq")

coef(LogisticModel)
summary(LogisticModel)

# Testing homogeneity and adequacy of the logistic regression SIMULTANEOUSLY
deviance(LogisticModel)
length(coef(LogisticModel)) # Number of parameters
length(P)                   # Number of observations
pchisq(deviance(LogisticModel), df=140-2, lower.tail = FALSE)

NullModel <- glm(Resp~1, family=binomial(link="logit"))
coef(NullModel)

# Testing whether the response to P is not 0 (under a logistic regression model)
anova(NullModel, LogisticModel, test="Chisq")

# Summing up
summary(LogisticModel)
confint(LogisticModel)

########################################################################
# Optional
# Performing some model control
#
########################################################################
RawResiduals <- residuals(FreeCurve, "response")
PearsonResiduals <- residuals(FreeCurve, "pearson")
Fitted <- fitted(FreeCurve)

plot(Fitted, RawResiduals)
abline(h=0)

plot(Fitted, PearsonResiduals)
abline(h=0)

# Checking the logistic regression curve
RawResiduals <- residuals(LogisticModel, "response")

plot(P, RawResiduals)
abline(h=0)
LL <- loess (RawResiduals~P)
Grid <- 1:600
Pred <- predict(LL, data.frame(P=Grid))
points(Grid, Pred, type="l", lwd=4, col="red")

detach(Ex3.3)

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