#################################################
   #                                               #
   # exemple de régression logistique              #
   #                                               #
   #         (gH)  gilles.hunault@univ-angers.fr   #
   #                                               #
   #################################################



   # 1. lecture des données pour lesquels on connait le groupe
   # ---------------------------------------------------------

   pg <- read.table("pg.dar",head=TRUE,row.names=1)

   print(pg)
   TAILLE <- pg$TAILLE
   GROUPE <- pg$GROUPE


   # 2. tracé des données
   # ---------------------

   png(file="reglogi1.png",width=800,height=600)
   plot(TAILLE,GROUPE,col=1+GROUPE,pch=16,xlim=c(130,200))
   rect(120,-0.1,210,1.1,col="lightblue1")
   par(new=TRUE)
   plot(TAILLE,GROUPE,col=1+GROUPE,pch=16,xlim=c(130,200))
   rug(TAILLE[GROUPE==0],side=1,ticksize=0.01)
   rug(TAILLE[GROUPE==1],side=3,ticksize=0.01)
   dev.off()


   # 3. Fonctions de prédiction
   # ---------------------------

   # 3.10 résultat de la régression (valeurs fournies par SAS)

   a <-  (   0.1812 ) # "SLOPE"
   b <-  ( -27.2103 ) # "INTERCEPT"

   # 3.11 Valeurs fournies par R

   rr <- glm(GROUPE~TAILLE,family=binomial)
   a  <- rr$coefficients[2]
   b  <- rr$coefficients[1]

   # 3.2 fonction de prédiction non arrondie

   pred <- function(x) {
      g <- a*x+b
      y <- exp(g)/(1+exp(g))
      return(y)
   } ; # fin de fonction pred

   # 3.3 fonction de prédiction arrondie

   preda <- function(x) { return( round(pred(x)) ) }


   # 4. prédiction des nouveaux individus
   # ------------------------------------

   # 4.1 lecture des données

   npg <- read.table("pginc.dar",head=TRUE)
   nx <- npg$TAILLE

   # 4.2 affectation et affichage

   nm <- matrix(nrow=3,ncol=2)
   rownames(nm) <- npg[,1]
   colnames(nm) <- c("TAILLE","GROUPE")
   nm[,1] <- nx
   nm[,2] <- preda(nx)
   print(nm)

   # ... mais bien sûr R dispose d'une fonction de prédiction
   #     et il aurait suffi d'écrire
   #
   #     nm[,2] <- predict(rr,newdata=data.frame(TAILLE=nx),type=c("response"))
   #

   # 4.3 rapport de cote ("odds ratio")

   odr <- exp( a)
   cat("\n le rapport de cote ou \"odds ratio\" vaut ",odr,"\n\n")

   # 5. Vérification graphique du modèle
   # -----------------------------------

   # 5.1 tracés des groupes, du prédicteur et des valeurs pour
   #     données sans groupe

   png(file="reglogi2.png",width=800,height=600)
   plot(TAILLE,GROUPE,col=1+GROUPE,pch=16,xlim=c(130,200))
   rect(120,-0.1,210,1.1,col="lightblue1")
   par(new=TRUE)
   plot(TAILLE,GROUPE,col=1+GROUPE,pch=16,xlim=c(130,200))
   yp <- pred(TAILLE)
   points(TAILLE,yp,col="green3",pch=16)
   points(nx,pred(nx),col="blue2",pch=16,lw=2)
   allTAILLE <- c(TAILLE,nx)
   xr <- seq(min(allTAILLE),max(allTAILLE),by=0.1)
   rug(TAILLE[GROUPE==0],side=1,ticksize=0.01)
   rug(TAILLE[GROUPE==1],side=3,ticksize=0.01)
   lines(xr,pred(xr),col="blue2")
   abline(h=0.5,col="orange",lw=2)
   dev.off()


   # 6. qu'aurait-donné une régression linéaire ?
   # ---------------------------------------------

   # 6.1 reprise des anciens graphiques

   png(file="reglogi3.png",width=800,height=600)
   plot(TAILLE,GROUPE,col=1+GROUPE,pch=16,xlim=c(130,200),ylim=c(-0.05,1.55))
   rect(120,-0.2,210,1.65,col=rgb(red=180/256,green=210/256,blue=190/256))
   par(new=TRUE)
   plot(TAILLE,GROUPE,col=1+GROUPE,pch=16,xlim=c(130,200),ylim=c(-0.05,1.55))
   abline(h=0,col="red",lw=2)
   abline(h=1,col="red",lw=2)
   points(TAILLE,yp,col="green3")
   lines(xr,pred(xr),col="blue2")
   abline(h=0.5,col="orange",lw=1)
   points(nx,pred(nx),col="blue2",pch=16,lw=2)
   rug(nx,side=1)
   rug(TAILLE[GROUPE==0],side=1,ticksize=0.01)
   rug(TAILLE[GROUPE==1],side=1,ticksize=0.01)

   # 6.2 régression linéaire et tracé :
   #     on voit bien qu'on sort de l'intervale [0,1]

   rl <- lm(GROUPE~TAILLE)$coefficients
   na <- as.real(rl[2])
   nb <- as.real(rl[1])
   yh <- na*allTAILLE + nb
   points(allTAILLE,yh,col="black",pch=16,lw=2)
   lines(allTAILLE,yh,col="black",pch=16,lw=1)

   nyh <- na*nx + nb

   points(nx,nyh,col="green",pch=16,lw=2)
   dev.off()