################################################################################
#  Solution of the exercise 3-2 on binomial models                             #
#  
# Last revised: Spring 2019
#
# Copyright © 2018 by Rodrigo Labouriau
################################################################################
# rm(list = ls()) # This remove ALL the objecs previously defined 
                  # (there is no way to retrieve them again!)

load("DataFramesStatisticalModelling.RData")

# str(Ex3.2)
attach(Ex3.2)

#------------------------------------------------------------------------------#
# Item a)
#------------------------------------------------------------------------------#

# Calculating the proportions

Prop <- n.gram/n.grains

plot(Field+runif(length(Field))/10-0.05, Prop,col=Trap)


# Ploting the proportions

prop<-tapply(n.gram, list(Trap,Field), sum)/tapply(n.grains, list(Trap,Field), sum)

points(1:4, prop[1, ], col=1, pch=19, cex=2, type="b")
points(1:4, prop[2, ], col=2, pch=19, cex=2, type="b")


# b) 

# Calculating the lodds (i.e. the logit transformed proportions)
Logit <-function(p){
 log(p/(1-p))
}

Logit(prop)

Lodds <- Logit(Prop)
plot(Field+runif(length(Field))/10-0.05, Lodds,col=Trap)

lodds <- Logit(prop)


points(1:4, lodds[1, ], col=1, pch=19, cex=2, type="b")
points(1:4, lodds[2, ], col=2, pch=19, cex=2, type="b")

# Ploting the lodds
plot(1:4,Logit(prop[1,]) ,type='b', pch='1', ylim=c(-3,2), ylab="Lodds", xlab="Field")
points(1:4,Logit(prop[2,]) ,type='b', pch='2')

# c) 

resp <- cbind(n.gram, n.grains-n.gram)

Trap <- factor(Trap)
Field <- factor(Field)

saturated <- glm(resp~factor(Obs.num), family=binomial)
deviance(saturated)


two.way <- glm(resp~Trap*Field-Trap-0, family=binomial)

summary(two.way)

ILogit<-function(l){exp(l)/(1+exp(l))}

coef(two.way)

ILogit(coef(two.way))

prop

deviance(two.way)

# d)

additive <- glm(resp~Trap+Field, family=binomial)

coef(additive)
coef(two.way)

deviance(additive)
deviance(two.way)

anova(additive, two.way, test='Chisq')

NoTrap <- glm(resp~Field, family=binomial)

anova(NoTrap,additive, test='Chisq')

NoField <- glm(resp~Trap, family=binomial)

anova(NoField,additive, test='Chisq')

Null <-  glm(resp~1, family=binomial)

anova(Null,additive, test='Chisq')

deviance(additive)

pchisq(deviance(additive), df=235, lower.tail=F)

# Cleaning up
detach(Ex3.2)
ls()

# The next call (commented out) cleans everything from the environment, 
# i.e. removes all the objects listed by the call "ls()"

# rm(list = ls())