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

# Loading the datasets required
load("/DataFramesStatisticalModelling.RData")

################################################################################
# Ex3.5 : Binomial models for studying the flora composition                   #
# Number a Weed out of Nplants sampled with rings in four places (Place)       #
# in two topographic positions (Position, 1 = low, 2 = high)                   #
################################################################################

attach(Ex3.5)

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

str(Ex3.5)

table(Place,Position)

prop <- tapply(Weed, list(factor(Place), factor(Position)),sum)/
  tapply(Nplants, list(factor(Place), factor(Position)),sum)  ; prop

plot(1:4, prop[,1], type='b', pch='1', ylim=range(prop), col='red',
     ylab='Proportion of weeds',
     xlab='Place'
)
points(1:4, prop[,2], type='b', pch='2', col='blue')
legend(1,0.8, legend=c('Low position','High position'), col=c('red','blue'),
       pch=c('1','2'), bty='n')


Logit <- function(p){
  return(log(p/(1-p)))
}

lodds <- Logit(prop) ;lodds

plot(1:4, lodds[,1], type='b', pch='1', ylim=range(lodds), col='red',
     ylab='Lodds of weeds',
     xlab='Place'
)
points(1:4, lodds[,2], type='b', pch='2', col='blue')
legend(1,0.8, legend=c('Low position','High position'), col=c('red','blue'),
       pch=c('1','2'), bty='n')

#------------------------------------------------------------------------------#
# Analysis using binomial models                                               #
#------------------------------------------------------------------------------#

resp <- cbind(Weed, Nplants-Weed)

# Model with interaction (full)
full <- glm(resp ~ factor(Place)*factor(Position) , family=binomial(link = "logit"))
summary(full)
deviance(full)
pchisq(deviance(full), df=length(Weed)-8, lower.tail=F)

# Additive model
additive <- glm(resp ~ factor(Place)+factor(Position) , family=binomial(link = "logit"))
anova(additive, full, test='Chisq')

# No effect of  Place
NOsoil <- glm(resp ~ factor(Position) , family=binomial(link = "logit"))
anova(NOsoil, additive, test='Chisq')

NOposition <- glm(resp ~ factor(Place) , family=binomial(link = "logit"))
anova(NOposition, additive, test='Chisq')

Null <- glm(resp ~ 1 , family=binomial(link = "logit"))
anova(Null, additive, test='Chisq')

# Some model control

RawResiduals <- residuals(additive, "response")
PearsonResiduals <- residuals(additive, "pearson")
Fitted <- fitted(additive)

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

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

detach(Ex3.5)
 
# The End!