#GLM 020 LOGISTIC REGRESSION
library(nlme)
setwd("c:/DATA/Models")
#========================================================== 
D=read.table("KNNL1401.txt",header=T)
D
FM1=glm(Y~X,data=D,family=binomial)
summary(FM1)

Results=cbind(X=D$X,Y=D$Y,Fitted=fitted(FM1))
Results
Residuals=cbind(X=D$X,Y=D$Y,
	Pearson=resid(FM1,type="pearson"),
	Deviance=resid(FM1,type="deviance"))
Residuals

#PLOTTING CURVE:
M=data.frame(X=D$X) #USING ORIGINAL INDEPENDENT VALUES
M=na.omit(M)
Pred=predict(FM1,newdata=M,type="response")
plot(x=D$X,y=D$Y,xlab="X",ylab="Y",pch=20)
points(M$X,Pred,pch=20,col="red")

summary(FM1)
anova(FM1,test="Chisq")
#==========================================================
B=read.table("Boar.txt",header=T)
B
FM2=glm(Tb~LengthCT,family=binomial(link="logit"),data=B) #(link="logit") DEFAULT
FM3=glm(Tb~LengthCT,family=binomial(link="probit"),data=B)
FM4=glm(Tb~LengthCT,family=binomial(link="cloglog"),data=B)

#PLOTTING CURVES:
M=data.frame(LengthCT=B$LengthCT) #USING ORIGINAL INDEPENDENT VALUES
M=na.omit(M)
#M=data.frame(LengthCT=seq(from=46.5,to=165,by=1)) #USING NEW SEQUENCE
Pred2=predict(FM2,newdata=M,type="response")
Pred3=predict(FM3,newdata=M,type="response")
Pred4=predict(FM4,newdata=M,type="response")
plot(x=B$LengthCT,y=B$Tb,xlab="Length",ylab="Tb",pch=20)
lines(M$LengthCT,Pred2,col="red")
lines(M$LengthCT,Pred3,col="blue")
lines(M$LengthCT,Pred4,col="green")
legend("topleft",legend=c("logit","probit","clog-log"),
	lty=c(1,1,1),col=c("red","blue","green"),cex=0.8)
 
#==========================================================
C=read.table("ParasiteCod.txt",header=T)
C=na.omit(C)
C$fArea=factor(C$Area)
C$fYear=factor(C$Year)
FM5=glm(Prevalence~fArea*fYear+Length,family=binomial,data=C)
summary(FM5)

C
anova(FM5,test="Chisq")
drop1(FM5,test="Chisq")

RM1=glm(Prevalence~fArea+fYear+Length,family=binomial,data=C)
anova(RM1,FM5,test="Chisq")

RM2=glm(Prevalence~fArea*fYear,family=binomial,data=C)
anova(RM2,FM5,test="Chisq")

summary(FM5)



























#==========================================================
I=read.table("Isolation.txt",header=T)
I
FM6=glm(incidence~area*isolation,family=binomial,data=I)
summary(FM6,test="Chisq")
anova(FM6,test="Chisq")

RM3=glm(incidence~area+isolation,family=binomial,data=I)
options(digits=9)
anova(RM3,FM6,test="Chisq")
	
#==========================================================	
D=read.table("KNNLt1403.txt", header=T)
D

FM7=glm(disease~se1+se2+sector+Age,family=binomial,data=D)	
summary(FM7)	
RM4=glm(disease~se1+se2+sector,family=binomial,data=D)
anova(RM4,FM7,test="Chisq")