cuse <- data.frame(matrix(c(
    1,   3,  61, 232,
    2,  80, 137, 400,
    3, 216, 131, 301,
    4, 268,  76, 203,
    5, 197,  50, 188,
    6, 150,  24, 164,
    7,  91,  10, 183), 7, 4, byrow=TRUE))

names(cuse) <- c("a","ster","other","none")
ages <- paste(seq(15,45,5),seq(19,49,5),sep="-")
cuse$ageg <- ages[cuse$a]
cuse
cuse$Y <- as.matrix(cuse[,c("none","ster","other")])

cuse$age <- seq(17.5,47.5,5)[cuse$a]
cuse$agesq <- cuse$age^2


# Pour ajuster un modèle logistique (Y à 2 modalités) : glm 
# Dans ce cas, il faut créer un variable Ynew à 2 modalités pour chaque (sous) modèle

# 1. Aucune/tout type de contraception
Ynew = cbind(apply(cuse$Y[,2:3],1,sum),cuse$Y[,1])
colnames(Ynew) = c("use","none")

cuse$Ynew = Ynew
mod = glm(Ynew~1,data=cuse,family=binomial)
summary(mod)
deviance(mod)
qchisq(.95,df=6)
1-pchisq(deviance(mod),df=df.residual(mod))
# pour faire plus joli
res.mod = c(deviance(mod),qchisq(.95,df=6),1-pchisq(deviance(mod),df=df.residual(mod)))
names(res.mod) <- c("Deviance","Seuil chi2","p-value")
res.mod


mod1 = glm(Ynew~age,data=cuse,family=binomial) # age continu
summary(mod1)

# Test qualité d'ajustement
# H_0  :  "le modèle est bon" | H_1  :  "le modèle n'est pas bon"
df = df.residual(mod1)
res.mod1 = c(deviance(mod1),qchisq(.95,df=df),1-pchisq(deviance(mod1),df=df)) 
names(res.mod1) <- c("Deviance","Seuil chi2","p-value")
res.mod1 # -> on conclut que le modèle n'est pas bon


# Amélioration  d'ajustement
# H_0 : le modèle 1 n'est pas meilleur que le modèle 0 | H_1 : le modèle 1 est meilleur que le modèle 0
D = deviance(mod)-deviance(mod1)
1-pchisq(D,df=1) # on conclut que le modèle 1 améliore la prédiction par rapport au modèle 0

AIC(mod) # Attention "aic" donnent des résultats étranges -> utiliser AIC
AIC(mod1)

plot(cuse$age,cuse$Ynew[,1]/(cuse$Ynew[,1]+cuse$Ynew[,2]),pch=20,xlab="age",ylab="P(Y=use)")
points(cuse$age,mod1$fitted.values,pch=18,col="red")
l = coef(mod1)[1]+coef(mod1)[2]*cuse$age
lines(cuse$age,exp(l)/(1+exp(l)),col="red")
#lines(cuse$age,mod1$fitted.values,col="red")

mod1g = glm(Ynew~ageg,data=cuse,family=binomial)
summary(mod1g)

deviance(mod1)
1-pchisq(deviance(mod1),df=5)


deviance(mod1g)

mod2 = glm(Ynew~age+agesq,data=cuse,family=binomial)
summary(mod2)
# Amélioration  d'ajustement
D = deviance(mod1)-deviance(mod2)
1-pchisq(D,df=1)
D0 = deviance(mod)-deviance(mod2)
1-pchisq(D,df=2)


# Test qualité d'ajustement
# H_0  :  "le modèle est bon" | H_1  :  "le modèle n'est pas bon"
deviance(mod2)
[1] 6.119465
1-pchisq(6.119465,df=4) # -> on conclut que le modèle est bon

anova(mod2,test="Chisq") # permet de faire les tests d'"amélioration" d'ajustement - Exemple  H_0 : le modèle 1 n'est pas meilleur que le modèle 0 | H_1 : le modèle 1 est meilleur que le modèle 0
# A chaque ligne on entre une nouvelle variable et on regarde si l'ajoût de cette variable améliore significativement le modèle ou non. 

AIC(mod) 
AIC(mod1)
AIC(mod2)

plot(cuse$age,cuse$Ynew[,1]/(cuse$Ynew[,1]+cuse$Ynew[,2]),pch=20,xlab="age",ylab="P(Y=use)")
lines(cuse$age,mod1$fitted.values,lty=3,col="red")
points(cuse$age,mod2$fitted.values,pch=18,col="red")
l = coef(mod2)[1]+coef(mod2)[2]*cuse$age+coef(mod2)[3]*cuse$agesq
lines(cuse$age,exp(l)/(1+exp(l)),col="red")

# 2. Contraception long terme/contraception court terme
Ynew = cbind(cuse$Y[,2],cuse$Y[,3])
colnames(Ynew) = c("long","short")
cuse$Ynew = Ynew

mod1 = glm(Ynew~age,data=cuse,family=binomial) 
mod2 = glm(Ynew~age+agesq,data=cuse,family=binomial) 
anova(mod2,test="Chisq")
plot(cuse$age,cuse$Ynew[,1]/(cuse$Ynew[,1]+cuse$Ynew[,2]),pch=20,xlab="age",ylab="P(Y=ster.)")
lines(cuse$age,mod1$fitted.values,lty=3,col="red")
points(cuse$age,mod2$fitted.values,pch=18,col="red")
l = coef(mod2)[1]+coef(mod2)[2]*cuse$age+coef(mod2)[3]*cuse$agesq
lines(cuse$age,exp(l)/(1+exp(l)),col="red")

c(AIC(mod) , AIC(mod1), AIC(mod2))

cuse$jeune = cuse$age<25
mod3 = glm(Ynew~age+agesq+jeune*age,data=cuse,family=binomial)
deviance(mod3)
plot(cuse$age,cuse$Ynew[,1]/(cuse$Ynew[,1]+cuse$Ynew[,2]),pch=20,xlab="age",ylab="P(Y=ster.)")
lines(cuse$age,mod1$fitted.values,lty=3,col="red")
points(cuse$age,mod2$fitted.values,pch=18,col="red")
l = coef(mod2)[1]+coef(mod2)[2]*cuse$age+coef(mod2)[3]*cuse$agesq
lines(cuse$age,mod2$fitted.values,col="red")
lines(cuse$age,mod3$fitted.values,col="blue")


cuse$agec = cuse$age-mean(cuse$age)
cuse$agec_abs = abs(cuse$agec)
mod2 = glm(Ynew~agec+agec_abs,data=cuse,family=binomial) 
plot(cuse$agec,cuse$Ynew[,1]/(cuse$Ynew[,1]+cuse$Ynew[,2]),pch=20,xlab="age",ylab="P(Y=ster.)")
lines(cuse$agec,mod1$fitted.values,lty=3,col="red")
points(cuse$agec,mod2$fitted.values,pch=18,col="red")
lines(cuse$agec,mod2$fitted.values,col="red")

dev.new()
plot(cuse$agec,log(cuse$Ynew[,1]/(cuse$Ynew[,1]+cuse$Ynew[,2])/(1-cuse$Ynew[,1]/(cuse$Ynew[,1]+cuse$Ynew[,2]))),pch=20,xlab="age",ylab="P(Y=ster.)")
lines(cuse$agec,log(mod1$fitted.values/(1-mod1$fitted.values)),lty=3,col="red")
points(cuse$agec,mod2$fitted.values,pch=18,col="red")
lines(cuse$agec,mod2$fitted.values,col="red")



# Pour ajuster un modèle multinomial (Y à plus de 2 modalités) : multinom (package nnet)
#install.packages("nnet")
library(nnet)
msat <- multinom(Y ~ ageg, data=cuse); msat
ml <- multinom(Y ~ age , data=cuse)
mlq <- multinom(Y ~ age + agesq, data=cuse)


obs <- log(cuse$Y[,2:3]/cuse$Y[,1])
age <- seq(15,50,0.5)
fit <- cbind(1,age,age^2) %*% t(coef(mlq))
color <- c("red","green")
method <- c("ster","other")
plot(cuse$age, obs[,1], type="n", 
      xlab="age", ylab="log-odds", xlim=c(15,55))
for (j in 1:2) {
   points(cuse$age, obs[,j], pch=19, col=color[j])
   lines(age, fit[,j], col=color[j])
   text(52, obs[7,j], label=method[j], col=color[j])
 }


# Deviance, comparaison modèle quadratique-modèle nul 
deviance(multinom(Y~1, data=cuse)) - deviance(mlq)

#The graph suggests that most of the lack of fit comes from overestimation of the relative odds of being sterilized compared to using no method at ages 15-19. Adding a dummy for this age group confirms the result:
mlqt <- multinom(Y ~ age + agesq + I(ageg=="15-19"), data=cuse)
x2 <- deviance(mlqt) - deviance(msat); x2
1-pchisq(x2, 6)