#GLM 020 LOGISTIC REGRESSION
library(nlme)
library(MASS)
setwd("c:/DATA/Models")

#========================================================== 
R=read.table("Roadkills.txt",header=T)
R
killed=R$TOT.N
distance=R$D.PARK
FM1=glm(killed~distance,family=poisson)
FM1a=glm.nb(killed~distance)
summary(FM1)

fit=fitted(FM1,type="response")
Pearson=resid(FM1,type="pearson")
Deviance=resid(FM1,type="deviance")
cbind(killed,distance,fit,Pearson,Deviance)

#PLOTTING FITS:
plot(distance,fitted(FM1,type="response"),col='red',
	ylab="Number Killed",type='l')
lines(distance,fitted(FM1a,type="response"),col='blue')
points(distance,killed,col='black',pch=20)

#========================================================== 
phi=sum(Pearson^2)/50
phi

FM2=update(FM1,family=quasipoisson)

#PLOTTING FITS & 95% CI:
PredFM1=predict(FM1,type="link",se=T)
PredFM2=predict(FM2,type="link",se=T)

U.FM1=exp(PredFM1$fit+1.96*PredFM1$se.fit)	
L.FM1=exp(PredFM1$fit-1.96*PredFM1$se.fit)
U.FM2=exp(PredFM2$fit+1.96*PredFM2$se.fit)	
L.FM2=exp(PredFM2$fit-1.96*PredFM2$se.fit)

plot(distance,killed,col='black',pch=20)
lines(distance,fitted(FM2,type="response"),col='red',
	ylab="Number Killed")
lines(distance,U.FM1,lty=2,col='red')	
lines(distance,L.FM1,lty=2,col='red')
lines(distance,U.FM2,lty=2,col='blue')	
lines(distance,L.FM2,lty=2,col='blue')

summary(FM2)

#========================================================== 
S=read.table("species.txt",header=T)
S

FM3=glm(Species~Biomass*pH,family=poisson,data=S)
summary(FM3)

RM1=glm(Species~Biomass+pH,family=poisson,data=S)
anova(FM3,RM1,test="Chisq")
drop1(FM3,test="Chisq")

#PLOTTING FITS:
plot(S$Biomass,S$Species,pch=20,col='black',
	ylab="Number of Species",
	xlab="Biomass")
points(S$Biomass[S$pH=="high"],S$Species[S$pH=="high"],pch=20,col="violet")
points(S$Biomass[S$pH=="mid"],S$Species[S$pH=="mid"],pch=20,col="blue")
points(S$Biomass[S$pH=="low"],S$Species[S$pH=="low"],pch=20,col="green")

FM1a=glm(Species[pH=="high"]~Biomass[pH=="high"],family=poisson,data=S)
FM1b=glm(Species[pH=="mid"]~Biomass[pH=="mid"],family=poisson,data=S)
FM1c=glm(Species[pH=="low"]~Biomass[pH=="low"],family=poisson,data=S)

points(S$Biomass[S$pH=="high"],fitted(FM1a),col="violet",pch=4)
points(S$Biomass[S$pH=="mid"],fitted(FM1b),col="blue",pch=4)
points(S$Biomass[S$pH=="low"],fitted(FM1c),col="green",pch=4)

#========================================================== 
openl=R$OPEN.L
monts=R$MONT.S
polic=sqrt(R$POLIC)
dpark=distance
shrub=sqrt(R$SHRUB)
watres=sqrt(R$WAT.RES)
lwatc=R$L.WAT.C
lproad=sqrt(R$L.P.ROAD)
dwatcour=sqrt(R$D.WAT.COUR)

T=data.frame(killed,openl,monts,polic,shrub,watres,lwatc,lproad,dwatcour,dpark)
T
FM4=glm(killed~openl+monts+polic+shrub+watres+lwatc+lproad+dwatcour+dpark,
	family=poisson,data=T)
summary(FM4)

FM5=update(FM4,family=quasipoisson)
summary(FM5)

RM2=glm(killed~dpark,family=poisson,data=T)

anova(FM5,RM2,test="F")
drop1(FM5,test="F")

#---------------------------------------------------------
anova(FM4)
R=residuals(FM4,type="response")
SSE=as.numeric(t(R)%*%R)
SSE
dfE=42
MSE=SSE/dfE
X=model.matrix(FM4)
Vb=MSE*solve(t(X)%*%X)
sd=sqrt(diag(Vb))
RESULTS=cbind(coefficients(FM4),sd)
RESULTS
summary(FM4)