################Sosk##############
data=matrix(scan(),byrow=T,ncol=3)
1.691 59 6 
1.724 60 13
1.755 62 18
1.784 56 28 
1.811 63 52
1.837 59 53
1.861 62 61
1.884 60 60

dim(data)
colnames(data)<-c("x","n","y")
data<-data.frame(data)
y=data$y
n=data$n
x=data$x

combdat=cbind(y,nminusy=n-y)
model1=glm(combdat~x,family=binomial(link=logit),data=data)
print(summary(model1))
residuals(model1) 
fitted.values(model1) 
predict.glm(model1) 
anova(model1)

model2=glm(combdat~x,family=binomial(link=probit),data=data) 
print(summary(model2)) 


model3=glm(combdat~x,family=binomial(link=cloglog),data=data) 
print(summary(model3)) 


R<-resid(model1)
R
shapiro.test(R). ## normality test
qqnorm(R)
qqline(R)

plot(fitted(model1),R)
library(lmtest)  
bptest(model1)  ## constant variance of residuals

ts.plot(R,type="b") ## 
abline(h=0,lty=2,col=2)

######################mice####################
data=matrix(scan(),byrow=T,ncol=4)
1.5 96 44 120
1.5 168 37 80
1.5 336 43 80
1.5 504 35 60
3.5 0.5 29 100
3.5 1 53 200
3.5 2 13 40 
3.5 3 75 200
3.5 5 23 40
3.5 7 152 280
3.5 14 55 80
3.5 24 98 140
3.5 48 121 160
7 0.5 52 120
7 1 62 120
7 1.5 61 120
7 2 86 120
 
dim(data)

colnames(data)<-c("no2","time","y","n")
data<-data.frame(data)
y=data$y
n=data$n
no2=data$no2
time=data$time
logno2=log(no2) 
logtime=log(time) 
logno2=(logno2-mean(logno2))/sqrt(var(logno2))
logtime=(logtime-mean(logtime))/sqrt(var(logtime)) 
combdat=cbind(y,nminusy=n-y)
model1=glm(combdat~logno2+logtime,family=binomial(link=logit),data=data)
print(summary(model1))
residuals(model1) 
rstudent(model1)
fitted.values(model1) ####pihat######
n*fitted.values(model1) ###yhat###
predict.glm(model1) ####hat-logit######
anova(model1)
model2=glm(combdat~logno2*logtime,family=binomial(link=logit),data=data) 
print(summary(model2)) 

######################snoring####################
y = c(24,35,21,30)
n = c(1379,638,213,254)
snore = c(0,2,4,5)

snoring<-data.frame(snore,y,n)
M.logit<-glm(cbind(y,n-y)~snore,family=binomial(link=logit),data=snoring)
M.logit
coef(M.logit)
confint(M.logit)
summary(M.logit)
fitted(M.logit)


M.probit<-glm(cbind(y,n-y)~snore,family=binomial(link=probit),data=snoring)
M.probit
fitted(M.probit)


M.lm<-glm(y/n~snore,weight=n,data=snoring)
M.lm
fitted(M.lm)

beta0<-coef(M.logit)[1]
beta1<-coef(M.logit)[2]
plot(snore,y/n,type="p",pch=16)
abline(coef(M.lm),lty=2,col=2)
curve(exp(beta0+beta1*x)/(1+exp(beta0+beta1*x)),0,5,add=T,col=4)
beta01<-coef(M.probit)[1]
beta11<-coef(M.probit)[2]
curve(pnorm(beta01+beta11*x),0,5,add=T,col=8)
legend("topleft",legend=c(NA,"lm","logit","probit"),lty=c(NA,2,1,1),col=c(NA,2,4,8),bty="n")

fitted(M.logit)
fitted(M.probit)
fitted(M.lm)

#####################Snoring Example #####################
snoring <-matrix(c(24,1355,35,603,21,192,30,224), ncol=2, byrow=TRUE)
dimnames(snoring) <-list(snore=c("never","sometimes","often","always"),heartdisease=c("yes","No"))
scores.a <- c(0,2,4,5) 
scores.b <- c(0,2,4,6) 
scores.c <- 0:3 
scores.d <- 1:4 
snoring.lg.a <- glm( snoring ~ scores.a, family=binomial() ) 
summary(snoring.lg.a)
snoring.lg.b <- glm( snoring ~ scores.b, family=binomial() ) 
summary(snoring.lg.b)
snoring.lg.c <- glm( snoring ~ scores.c, family=binomial() ) 
summary(snoring.lg.c)
snoring.lg.d <- glm( snoring ~ scores.d, family=binomial ()) 
summary(snoring.lg.d)
coef(snoring.lg.a)
coef(snoring.lg.b)
coef(snoring.lg.c )
coef(snoring.lg.d) 
predict(snoring.lg.a, type="response" )
predict(snoring.lg.b, type="response" )
predict(snoring.lg.c, type="response" )
predict(snoring.lg.d, type="response" )
snoring.logit <-glm( snoring ~ scores.a, family=binomial(link="logit"))
summary(snoring.logit)
predict(snoring.logit, type="response")


#################################################################
## glm 
snore.f<-factor(1:4,labels=c("never","sometime","almost","all"))
snore.f

M<-glm(cbind(y,n-y)~snore.f,family=binomial(link=logit))
M
summary(M)
exp(coef(M))
exp(confint(M))


anova(M,test="Chisq")
anova(M,test="LRT")
anova(M.logit,M,test="Chisq")
anova(M.probit,M,test="Chisq")
anova(M.logit,M.probit,test="Rao")

######################poisson regression####################
data<-read.table(file="/Users/macbook/Documents/Fatemeh/Fati-Old-PC/examine/categgorical data/R-program/poisson regression/data.txt")
dim(data)
y=data$V5
x1=data$V2
x2=data$V3
x3=data$V4

x2<-as.factor(x2)

model1=glm(y~x1+x2+x3,family=poisson(link=log))
summary(model1)

c.x2<-coef(model1)[3]

exp(c.x2) # odds ration 
1/exp(c.x2)

exp(confint(model1))



x22<-relevel(x2,ref="small")
x22

model2=glm(y~x1+x22+x3,family=poisson(link=log))
summary(model2)