dat<-read.table("C:\\Users\\Grassi\\Desktop\\TABLE_5_10.txt",header=TRUE)

# Stima del modello di regressione logistica Tabella 5.10 pag 198
Logit_5_10 <- glm(y~x,family = binomial(link=logit),data=dat)
summary(Logit_5_10)

# punto a, valori previsti nei termini di probabilità:

B<-coefficients(Logit_5_10)
X <- as.matrix( data.frame(1,dat[,1]) )

logit_hat<-X%*%B
prob_hat<-exp(logit_hat)/(1+exp(logit_hat)) # ANTILOGIT

fitted.values(Logit_5_10)
predict(Logit_5_10)

# punto b, valori previsti (x=26):

logit_26<-c(1,26)%*%B
prob_26<-exp(logit_26)/(1+exp(logit_26))

# punto c, derivata prima locale della funzione logistica
X[1,] # x ==> 26
B[2]*prob_hat[1]*(1-prob_hat[1]) # =  0.009177602 
B[2]*prob_26[1]*(1-prob_26[1]) # =  0.03621476



# punto d
logit_14<-c(1,14)%*%B
prob_14<-exp(logit_14)/(1+exp(logit_14))

logit_28<-c(1,28)%*%B
prob_28<-exp(logit_28)/(1+exp(logit_28))

diff(c(prob_28,prob_14))

# punto e

exp(B[1]);exp(B[2])

logit_1<-c(1,1)%*%B
exp(logit_1)
prob_1<-exp(logit_1)/(1+exp(logit_1))
prob_1/(1-prob_1)

exp(logit_1)/
exp(B[1])



# punto fg, CI95% odds ratio

O.R.<-exp(B[2])

B[2]+c( qnorm(0.025),- qnorm(0.025))*0.05933801
confint(Logit_5_10)[2,]
exp(confint(Logit_5_10)[2,])

Wald.CI95.l<-exp( B[2]-1.96*sqrt(0.003521)  )
Wald.CI95.u<-exp( B[2]+1.96*sqrt(0.003521)  )

round(Wald.CI95.l,3)
round(Wald.CI95.u,3)

# Wald test

z1=B[1]/sqrt(1.900616)
z1;z1^2

pnorm(z1)*2

CovM<-vcov(Logit_5_10)
B[1]/sqrt(CovM[1,1])

z2<-B[2]/sqrt(CovM[2,2])
pnorm(z2,lower.tail=FALSE)*2 #pvalore

#Chi - squares

v1<-c(1,0)
v2<-c(0,1)

Chisq_1<-(v1%*%B)%*%solve(v1%*%CovM%*%v1)%*%(v1%*%B)

Chisq_2<-(v2%*%B)%*%solve(v2%*%CovM%*%v2)%*%(v2%*%B)

1-pchisq(c(Chisq_1,Chisq_2),df=1)

# punto h, Lik. R. Test: -2*log.LikH1 e -2*log.LikH0
 names(Logit_5_10)
 Logit_5_10$deviance
 Logit_5_10$null.deviance
 
 Logit_5_10$df.null
 Logit_5_10$df.residual

Chisq_test<- Logit_5_10$null.deviance - Logit_5_10$deviance
df_test<- 26-25

Chisq_test
1-pchisq(Chisq_test,df=df_test) # area coda destra Chisq, 1 gdl
anova(Logit_5_10,test="Chisq")

# confidence interval for fit LI=8 
#Chi - squares
logit_8<-c(1,8)%*%B
exp(logit_8)
prob_8<-exp(logit_8)/(1+exp(logit_8))
prob_8

v<-c(1,8)
s.err<-sqrt(v%*%CovM%*%v)

deltamethod.CI95.l<- logit_8-1.96*sqrt( v%*%CovM%*%v )  
deltamethod.CI95.u<- logit_8+1.96*sqrt( v%*%CovM%*%v )  

# --------------
prob_8
# -- lower 95% -
exp(deltamethod.CI95.l)/(1+exp(deltamethod.CI95.l))
# -- upper 95% -
exp(deltamethod.CI95.u)/(1+exp(deltamethod.CI95.u))