# RegLog.r

library(Design)
library(mlogit)

##############
# Exemple GPEM

load("GPEM.rda")

RegLog_GPEM <- glm(relevel(M.STATUS,2)~relevel(EMS,2)+relevel(PMS,2)+relevel(EMS,2):relevel(PMS,2),data=GPEM,family=binomial)
summary(RegLog_GPEM)


# Remarque: Il est possible de spécifier plus simplement le même modèle de la manière suivante:
summary(glm(relevel(M.STATUS,2)~relevel(EMS,2)*relevel(PMS,2),data=GPEM,family=binomial))

# Codages
# Sum
contrasts(GPEM$EMS) = contr.sum(2)
contrasts(GPEM$PMS) = contr.sum(2)
summary(glm(relevel(M.STATUS,2)~EMS*PMS,data=GPEM,family=binomial))

# Helmert
contr.helmert(3)
contr.helmert(4)


# Test de la déviance
1-pchisq(1331.3,1032)


# Utilisation de lrm()
RegLog_GPEM_2 <- lrm(relevel(M.STATUS,2)~relevel(EMS,2)*relevel(PMS,2),data=GPEM)
RegLog_GPEM_2


# Odds ratio
exp(RegLog_GPEM$coefficients)
exp(RegLog_GPEM_2$coefficients)


###########################################
# Exemple Revenu, opinion politique, région

load("RevOpReg.rda")

summary(RevOpReg)

contrasts(RevOpReg$Politique) = contr.SAS(3)
contrasts(RevOpReg$Region) = contr.SAS(2)

# Modèle de départ avec interaction
ROR_1 <- glm(Revenu~Politique*Region,data=RevOpReg,family=binomial)
summary(ROR_1)
1-pchisq(827.68,924)

# Suppression de l'interaction
ROR_2 <- glm(Revenu~Politique+Region,data=RevOpReg,family=binomial)
summary(ROR_2)
1-pchisq(827.79,926)

# Suppression de Politique
ROR_3 <- glm(Revenu~Region,data=RevOpReg,family=binomial)
summary(ROR_3)
1-pchisq(1220.7,928)
# --> Ce modèle est rejeté. Le meilleur est ROR_2.

# Suppression de Region (ne figure pas dans le polycopié, mais devrait!!!)
ROR_3b <- glm(Revenu~Politique,data=RevOpReg,family=binomial)
summary(ROR_3b)
1-pchisq(837.82,927)
# --> Ce modèle est acceptable.

# Procédure automatique avec BIC
RevopReg_BIC <- step(ROR_1,direction="backward",k=log(930))
# --> Meilleur modèle: Politique + Region.


#################
# Exemple Travail

# Data
load("Travail.rda")

summary(Travail)

# Modèle avec les 7 variables explicatives
Travail_1 <- glm(THISYR~LASTYR+BLACK+CHILD1+CHILD2+EDUC+AGE+HUBINC,data=Travail,family=binomial)
summary(Travail_1)

# Sélection basée sur les p-valeurs
Travail_2 <- glm(THISYR~LASTYR+BLACK+CHILD1+CHILD2+EDUC+AGE,data=Travail,family=binomial)
summary(Travail_2)

Travail_3 <- glm(THISYR~LASTYR+CHILD1+CHILD2+EDUC+AGE,data=Travail,family=binomial)
summary(Travail_3)

Travail_4 <- glm(THISYR~LASTYR+CHILD1+EDUC+AGE,data=Travail,family=binomial)
summary(Travail_4)

Travail_5 <- glm(THISYR~LASTYR+EDUC+AGE,data=Travail,family=binomial)
summary(Travail_5)


# Procédures de sélection automatique basées sur AIC et BIC
Travail_AIC <- step(Travail_1,direction="backward")
summary(Travail_AIC)

Travail_BIC <- step(Travail_1,direction="backward",k=log(200))
summary(Travail_BIC)


#########################################
# Graphique distributions logit et probit

x <- seq(-5, 5, length.out=100)
y1 <- plogis(.x, location=0, scale=1)
y2 <- pnorm(.x, mean=0, sd=1)

win.graph(width=8, height=4)

plot(x,y1,lwd=2,type="l",xlab="",ylab="")

abline(h=0, col="gray")
abline(h=1, col="gray")
abline(v=0, col="gray")
lines(x,y2,col="black",lwd=2,lty=2)

legend(-4.5,0.8, c("Logit", "Probit"),lwd=2, lty = c(1,2))


#########################
# Régression multinomiale

# Exemple Revenu, opinion politique, région

load("RevOpReg.rda")

summary(RevOpReg)

contrasts(RevOpReg$Revenu) = contr.SAS(2)
contrasts(RevOpReg$Region) = contr.SAS(2)

# Vraie régression multinomiale (les 2 régressions sont estimées simultanément)
MD <- mlogit.data(RevOpReg, varying=NULL, choice="Politique", shape="wide")
ROR_1 <- mlogit(Politique~1|Revenu+Region,data=MD,reflevel="Gauche")
summary(ROR_1)

# Régression multinomiale utilisant deux régressions estimées séparément
ROR_1a <- mlogit(Politique~1|Revenu+Region,data=MD,reflevel="Gauche",alt.subset=c("Gauche","Centre"))
summary(ROR_1a)
ROR_1b <- mlogit(Politique~1|Revenu+Region,data=MD,reflevel="Gauche",alt.subset=c("Gauche","Droite"))
summary(ROR_1b)


#####################
# Régression ordinale

# Exemple Revenu, opinion politique, région

load("RevOpReg.rda")

summary(RevOpReg)

contrasts(RevOpReg$Revenu) = contr.SAS(2)
contrasts(RevOpReg$Region) = contr.SAS(2)

ROR_ord <- lrm(Politique~Revenu+Region,data=RevOpReg)
ROR_ord