#GLM 030 LOGISTIC REGRESSION WITH PROPORTIONS
library(nlme)
setwd("c:/DATA/Models")

#========================================================== 
S=read.table("sexratio.txt",header=T)
S
Y=cbind(S$males,S$females)
logdensity=log(S$density)
FM1=glm(Y~logdensity,family=binomial)
summary(FM1)
Results=cbind(Y=Y,logdensity=logdensity,Fitted=fitted(FM1,type="response"))
Results
Pmales=S$males/(S$males+S$females)
Residuals=cbind(logdensity=logdensity,Pmales=Pmales,
	Pearson=resid(FM1,type="pearson"),
	Deviance=resid(FM1,type="deviance"))
Residuals
sum((resid(FM1,type='deviance')^2))

#PLOTTING REGRESSION FIT:
Pmales=S$males/(S$males+S$females)
plot(logdensity,Pmales,pch=20,col="black")
points(logdensity,predict(FM1,type="response"),pch=20,col='red')

#========================================================== 
#ZUUR ET AL. 10.3 DEER TB DATA:
T=read.table("TBdeer.txt",header=T)
T

T$fFenced=factor(T$Fenced)
pos=T$DeerPosCervi
neg=T$DeerSampledCervi-T$DeerPosCervi

FM2=glm(cbind(pos,neg)~
       OpenLand+ScrubLand+QuercusPlants+QuercusTrees+
       ReedDeerIndex+EstateSize+fFenced,
       family=binomial,data=T)
summary(FM2)

FM3=update(FM2,family=quasibinomial)
summary(FM3)
drop1(FM3,test="F")
anova(FM3,test="F")

RM1=glm(cbind(pos,neg)~OpenLand,family=quasibinomial,data=T)
summary(RM1)
drop1(RM1,test="F")

#PLOTTING FIT:
M=data.frame(OpenLand=seq(from = min(T$OpenLand),to = max(T$OpenLand),by=0.01))

Pred=predict(RM1,newdata=M,type="response",se=TRUE)
plot(M$OpenLand,Pred$fit,
     type="l",ylim=c(0,1),col='red',
     xlab="Percentage open land",
     ylab="E. cervi Probability")
PROPORTION=T$DeerPosCervi/T$DeerSampledCervi
points(T$OpenLand,PROPORTION,pch=20)
lines(M$OpenLand,Pred$fit+1.96*Pred$se.fit,lty=2,col='blue')
lines(M$OpenLand,Pred$fit-1.96*Pred$se.fit,lty=2,col='blue')

#==========================================================