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

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

################################################################################
# Ex3.6 : Logistic regression for studying fish                                #
# Dead fishes, 6 doses antibiotic treatment, 2 Ecotypes                        #
################################################################################

attach(Ex3.6)

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

str(Ex3.6)

table(Doses, Ecotype)

prop <- tapply(DeadFishes, list(factor(Doses), factor(Ecotype)),sum)/
	tapply(Nfishes, list(factor(Doses), factor(Ecotype)),sum) ; prop

plot(levels(factor(Doses)), prop[,1], type='b', pch='1', ylim=range(prop),
		 ylab='Proportion of dead fishes', xlab='Dosis of antibiotics', col='red')
points(levels(factor(Doses)), prop[,2], type='b', pch='2', col='blue')
legend(40,0.8, legend=c('Ecotype 1','Ecotype 2'), col=c('red','blue'),
			 pch=c('1','2'), bty='n')

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

lodds <- Logit(prop) ;lodds

plot(levels(factor(Doses)), lodds[,1], type='b', pch='1', ylim=range(lodds),
		 ylab='Lodds of dead', xlab='Dosis of antibiotics', col='red')
points(levels(factor(Doses)), lodds[,2], type='b', pch='2')
legend(40,2, legend=c('Ecotype 1','Ecotype 2'), col=c('red','blue'),
			 pch=c('1','2'), bty='n')

#------------------------------------------------------------------------------#
# Modeling the fish mortality via logistic regressions                         #
#------------------------------------------------------------------------------#

resp <- cbind(DeadFishes, Nfishes-DeadFishes)
free.curves <- glm(resp ~ factor(Doses)*factor(Ecotype), family=binomial(link='logit'))
deviance(free.curves)
pchisq(deviance(free.curves), df=length(Doses)-12, lower.tail=F)

two.logistic.curves <- glm(resp ~ Doses*factor(Ecotype),
													 family=binomial(link='logit'))
summary(two.logistic.curves)
anova(two.logistic.curves, free.curves, test='Chisq')

parallel <- glm(resp ~ Doses+factor(Ecotype), family=binomial(link='logit'))
anova(parallel, two.logistic.curves, test='Chisq')

NoDoses <- glm(resp ~ factor(Ecotype), family=binomial(link='logit'))
anova(NoDoses, parallel, test='Chisq')

NoEcotype <- glm(resp ~ Doses, family=binomial(link='logit'))
anova(NoEcotype , parallel, test='Chisq')

detach(Ex3.6)

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

# End!