Valid XHTML     Valid CSS2    

Logiciel R, séances de perfectionnement

    en 4 demi-journées

Séance 3 : Modélisation (régression linéaire, logistique, simple, multiple)

                     gilles.hunault "at" univ-angers.fr

 

Table des matières cliquable

  1. Notion de régression (simple, multiple, linéaire, logistique)

  2. Régression linéaire simple

  3. Régression linéaire multiple

  4. Régression logistique simple

  5. Régression logistique multiple

 

Il est possible d'afficher toutes les solutions via ?solutions=1 et de toutes les masquer avec via ?solutions=0.

 

1. Notion de régression (simple, multiple, linéaire, logistique)

Comment effectue-t-on une régression en R ? Quelles sont les principales fonctions génériques ?

Solution :  

Pour réaliser une régression linéaire en R, on utilise la fonction lm avec l'opérateur ~ (tilde) alors qu'une régression logistique passe par la fonction glm.

La fonction lm renvoie une liste avec entre autres les variables suivantes :

    $coefficients         estimation des coefficients    
    $fitted.values         les valeurs calculées pour chaque y   
    $residuals         les résidus (différences entre valeurs originales et valeurs estimées)   

Si on écrit print(ml) où ml est le résultat d'une fonction lm ou glm, on n'obtient que les coefficients de la régression. Pour avoir le détail de la régression, il faut utiliser summary(ml) et confint(ml) pour avoir les intervalles de confiance des coefficients de la régression. Les fonctions correspondantes sont : summary.lm et summary.glm car summary est une fonction générique qui peut s'appliquer à de nombreux objets. On peut ensuite prédire de nouvelles valeurs avec la fonction predict qui est une aussi fonction générique qui se décline en predict.lm et predict.glm.

La partie graphique est très importante pour les régressions. Pour une régression linéaire, après un plot élémentaire des données, un abline du modèle trace la droite de régression. A la suite d'une régression, la fonction plot.lm produit 6 graphiques pour analyser la qualité de la régression, dont seuls 4 sont affichés par défaut, repérés ci-dessous par une étoile :

            x         y    
   1       * Fitted values       Residuals   
   2       * Theoretical quantiles       Standardized residuals   
   3       * Fitted values       Root of standardized residuals   
   4         Observation number       Cook's distance   
   5       * Leverage       Standardized residuals   
   6         Leverage hn       Cook's distance   

Enfin, signalons qu'avec predict et matlines il est possible de tracer l'intervalle de confiance de la droite et l'intervalle de confiance des prévisons.

 

2. Régression linéaire simple

Après avoir vu la documentation de la fonction lm() du package stats essayer de réaliser la modélisation par une relation linéaire la dépendance entre la variable consommation d'essence et la variable distance parcourue dans le jeu de données km.dar (sans transformation des données).

Dans quel ordre doit-on lire les résultats  ? Faut-il regarder en premier le R2, la p-value de F, les coefficients ?

Y a-t-il des représentations graphiques associées ? Pourquoi y a-t-il 4 graphiques produits avec plot(lm(modele)) alors qu'on pourrait en avoir 6 ? Comment obtenir les 6 tracés ? A quoi correspondent-ils ?

Faut-il envisager une relation de causalité ? Quelles sont les valeurs prédites par le modèle pour 100 et 250 km ? Quelle est l'augmentation de la consommation si la distance augmente de 1 km ?

Solution :  

Rappel : pour définir un modèle linéaire avec R, il faut utiliser la fonction lm() et la droite correspondante se trace directement avec abline(). Pour afficher le modèle, on peut utiliser les fonctions génériques summary() et coef().


     # chargement des fonctions (gH)     
          
     source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
          
     # lecture des données (enregistrées en local)     
          
     km  <- lit.dar("km.dar")     
          
     # lecture des données sur internet     
     # km <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/km.dar")     
          
     idx <- order(km$distance)     
          
     # tracé rapide     
          
     # surtout pas     
     #  with( km, plot(essence~distance,type="b") )     
          
     with( km, plot(essence[idx]~distance[idx],type="p",pch=19,col="red" ))     
          
     # ajout de la droite de régression et un titre     
          
     ml <- lm(essence ~ distance,data=km) # pas besoin de idx ici     
     abline(ml, col="blue",lwd=2)     
     title("Consommation d'essence en fonction de la distance")     
          
     # affichons le nom des points     
          
     text(x=km$distance,y=km$essence+2,label=1:nrow(km))     
          
     # oups, le point 4 déborde, ajoutons-le     
          
     text(x=400,y=33,label="4")     
          
     # description du modèle     
          
     cats("Modèle linéaire ESSENCE en fonction de DISTANCE")     
     print( summary(ml) )     
          
     cats("ANOVA de ce modèle")     
     print( anova(ml) )     
          
     # coefficients     
          
     cats("coefficients seulement")     
     print( coef(ml) )     
          
     # l'objet issu de lm est complexe     
          
     cats("Objet modèle linéaire :")     
     print( names(ml) )     
     print( str(ml) )     
          

          
     Modèle linéaire ESSENCE en fonction de DISTANCE     
     ===============================================     
          
          
     Call:     
     lm(formula = essence ~ distance, data = km)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -1.6824 -0.9326 -0.6783  0.0686  4.4836     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept) 3.350376   0.894150   3.747  0.00458 **     
     distance    0.083320   0.004985  16.714 4.39e-08 ***     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 1.835 on 9 degrees of freedom     
     Multiple R-squared:  0.9688,	Adjusted R-squared:  0.9653     
     F-statistic: 279.4 on 1 and 9 DF,  p-value: 4.393e-08     
          
          
     ANOVA de ce modèle     
     ==================     
          
     Analysis of Variance Table     
          
     Response: essence     
               Df Sum Sq Mean Sq F value    Pr(>F)     
     distance   1 940.61  940.61  279.36 4.393e-08 ***     
     Residuals  9  30.30    3.37     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     coefficients seulement     
     ======================     
          
     (Intercept)    distance     
      3.35037574  0.08331991     
          
     Objet modèle linéaire :     
     =======================     
          
      [1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values" "assign"        "qr"            "df.residual"   "xlevels"       "call"          "terms"     
     [12] "model"     
     List of 12     
      $ coefficients : Named num [1:2] 3.3504 0.0833     
       ..- attr(*, "names")= chr [1:2] "(Intercept)" "distance"     
      $ residuals    : Named num [1:11] -0.6824 -0.0144 1.6537 -0.6783 0.1516 ...     
       ..- attr(*, "names")= chr [1:11] "traj0001" "traj0002" "traj0003" "traj0004" ...     
      $ effects      : Named num [1:11] -50.051 -30.669 1.756 -0.609 0.304 ...     
       ..- attr(*, "names")= chr [1:11] "(Intercept)" "distance" "" "" ...     
      $ rank         : int 2     
      $ fitted.values: Named num [1:11] 11.7 20 28.3 36.7 15.8 ...     
       ..- attr(*, "names")= chr [1:11] "traj0001" "traj0002" "traj0003" "traj0004" ...     
      $ assign       : int [1:2] 0 1     
      $ qr           :List of 5     
       ..$ qr   : num [1:11, 1:2] -3.317 0.302 0.302 0.302 0.302 ...     
       .. ..- attr(*, "dimnames")=List of 2     
       .. .. ..$ : chr [1:11] "traj0001" "traj0002" "traj0003" "traj0004" ...     
       .. .. ..$ : chr [1:2] "(Intercept)" "distance"     
       .. ..- attr(*, "assign")= int [1:2] 0 1     
       ..$ qraux: num [1:2] 1.3 1.19     
       ..$ pivot: int [1:2] 1 2     
       ..$ tol  : num 1e-07     
       ..$ rank : int 2     
       ..- attr(*, "class")= chr "qr"     
      $ df.residual  : int 9     
      $ xlevels      : Named list()     
      $ call         : language lm(formula = essence ~ distance, data = km)     
      $ terms        :Classes 'terms', 'formula' length 3 essence ~ distance     
       .. ..- attr(*, "variables")= language list(essence, distance)     
       .. ..- attr(*, "factors")= int [1:2, 1] 0 1     
       .. .. ..- attr(*, "dimnames")=List of 2     
       .. .. .. ..$ : chr [1:2] "essence" "distance"     
       .. .. .. ..$ : chr "distance"     
       .. ..- attr(*, "term.labels")= chr "distance"     
       .. ..- attr(*, "order")= int 1     
       .. ..- attr(*, "intercept")= int 1     
       .. ..- attr(*, "response")= int 1     
       .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>     
       .. ..- attr(*, "predvars")= language list(essence, distance)     
       .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"     
       .. .. ..- attr(*, "names")= chr [1:2] "essence" "distance"     
      $ model        :'data.frame':	11 obs. of  2 variables:     
       ..$ essence : int [1:11] 11 20 30 36 16 12 7 8 5 10 ...     
       ..$ distance: int [1:11] 100 200 300 400 150 50 60 70 20 100 ...     
       ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 essence ~ distance     
       .. .. ..- attr(*, "variables")= language list(essence, distance)     
       .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1     
       .. .. .. ..- attr(*, "dimnames")=List of 2     
       .. .. .. .. ..$ : chr [1:2] "essence" "distance"     
       .. .. .. .. ..$ : chr "distance"     
       .. .. ..- attr(*, "term.labels")= chr "distance"     
       .. .. ..- attr(*, "order")= int 1     
       .. .. ..- attr(*, "intercept")= int 1     
       .. .. ..- attr(*, "response")= int 1     
       .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>     
       .. .. ..- attr(*, "predvars")= language list(essence, distance)     
       .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"     
       .. .. .. ..- attr(*, "names")= chr [1:2] "essence" "distance"     
      - attr(*, "class")= chr "lm"     
     NULL     
          

km.png

Il est d'usage de commencer par regarder si le modèle est significativement différent de zéro dans les sorties de summary() grâce à la p-value du test de Fisher. Si cette p-value est supérieure à 0,05 «le modèle ne vaut rien». Si le modèle est significativement différent de zéro, on peut s'intéresser directement au coefficient de détermination en régression linéaire simple. Par contre en régression linéaire multiple, il faut regarder les p-values des coefficients pour envisager une sélection de variables avant d'utiliser le coefficient de détermination ajusté.

Ici, le modèle est significativement différent de zéro (p=4.393e-08<<0,05) et on peut donc dire ici que la consommation est d'environ 0,08332 litre par km, ou, pour parler plus couramment, d'environ 8,332 litre "aux 100". Donc si la distance augmente d'un km, la consommation augmente de 0,008332 litre. Par contre il ne faut pas mésestimer le terme constant, soit ici 3,53 l qui doit correspondre au premier démarrage de la voiture (il s'agit de ma propre voiture, une XM et il ne faut pas rouler tant que le systéme hydraulique n'est pas opérationnel, donc tant que la voiture n'est pas en "position haute", ce qui dure au moins 30 secondes...) et à la consommation pendant l'arrêt. Cette constante sera d'autant plus importante que la distance est effectuée en ville, à cause des stops et des feux rouges, mais il est clair que ces interprétations ne sont pas statistiques et font appel à des connaissances contextuelles liées aux données. Le modèle est un "très bon modèle" puisque le coefficient de détermination vaut 0.97 (donc très proche de 1).

Dans la mesure où un modèle linéaire est [presque] toujours calculable, il est important de savoir si la modélisation est significativement non nulle, c'est-à-dire si on peut accepter le modèle. C'est pour cela qu'il faut utiliser le test de Fisher (de l'analyse de la variance) appliqué au modèle, et ensuite regarder le coefficient de corrélation linéaire. Il ne faut pas accorder une confiance aveugle à ces tests, qui peuvent être mis en défaut notamment pour de petits échantillons. On doit donc systématiquement réaliser une analyse des résidus, ce qui peut souvent se réduire à des tracés diagnostics que fournit la fonction plot.lm() du package stats. C'est ce qui explique que R fournit 6 tracés (dont 4 seulement par défaut) pour plot.lm(modele) qui est équivalent à plot(lm(modele)).

Pour obtenir les 6 tracés, il faut utiliser le paramètre which :


     # lecture des données     
          
     km  <- lit.dar("km.dar")     
          
     # tracé des 6 graphiques de diagnostics pour le modèle     
          
     ml <- lm(essence ~ distance,data=km)     
          
     for (idg in (1:6)) {     
       plot(ml,which=idg)     
     } # fin pour idg     
          
     # tracés et stockage des tracés dans des fichiers PNG     
          
     par(ask=FALSE)     
          
     for (idg in (1:6)) {     
       fgr <- paste("kmdiagno_",idg,".png",sep="")     
       gr(fgr)     
       plot(ml,which=idg)     
       dev.off()     
       cat("fichier ",fgr," généré.\n")     
     } # fin pour idg     
          
          
kmdiagno_1.png

          

kmdiagno_2.png

          

kmdiagno_3.png
kmdiagno_4.png

          

kmdiagno_5.png

          

kmdiagno_6.png

Les graphiques obtenus se nomment tracés de diagnostic du modèle et viennent analyser les résidus, les points-leviers, les outliers... Ici par exemple, les trois premiers graphiques montrent que les trajets 3, 6 et 10 semblent ne pas bien être en accord avec le modèle.

Avec la fonction predict() appliquée à un data frame avec les mêmes noms de variables, on peut utiliser le modèle linéaire pour prédire de nouvelles valeurs :


     # lecture des données     
          
       km  <- lit.dar("km.dar")     
          
     # régression linéaire     
          
       ml <- lm(essence ~ distance,data=km)     
          
     # prédictions     
          
       cats("valeurs prédites")     
          
       newx <- c(100,250)     
          
     # surtout pas     
     #   newy   <- predict(ml,newx)     
     # car cela produit l'erreur     
     #   Erreur dans eval(predvars, data, env) :     
     #   argument numérique 'envir' n'ayant une longueur unitaire     
          
     # on construit un nouveau data frame avec le même nom pour x     
          
       newxdf <- data.frame(distance=newx)     
       newy   <- predict(ml,newxdf)     
          
       print(rbind(newx,newy))     
          
     # tracé avec ajout des points prédits et des intervalles de confiance     
          
       co      <- round(coef(ml),3)     
       idx     <- order(km$distance)     
       formule <- paste("Consommation = ", 100*co[2]," x distance/100 + ",co[1],sep="")     
          
       gr("kmpreds.png")     
          
       par(omi=c(0,1,0,0)) # pour avoir plus place à gauche de l'axe     
          
       with( km,     
        plot(essence[idx]~distance[idx],type="p",     
          pch=19,col="blue",xlim=c(0,800),ylim=c(0,60),     
          xlab="distance (kilomètres)",ylab="essence (litres)",main=formule     
        ) # fin de plot     
       ) # fin de with     
          
       abline(ml,lwd=2,col="black") # droite de régression (moindres carrés ordinaires)     
          
     # points pour lesquels on veut une prédiction     
          
       points(cbind(newx,newy),cex=1.5,col="red",pch=19)     
       lines(rep(newx[1],2),c(-10,newy[1]),lty=3,col="blue",lwd=2)     
       lines(rep(newx[2],2),c(-10,newy[2]),lty=3,col="blue",lwd=2)     
       lines(c(-50,newx[1]),rep(newy[1],2),lty=3,col="blue",lwd=2)     
       lines(c(-50,newx[2]),rep(newy[2],2),lty=3,col="blue",lwd=2)     
          
     # tracé des intervalles de confiance     
          
       grillex <- data.frame(distance=1:800)     
       icConf  <- predict(ml,new=grillex,interval="confidence",level=0.95)     
       icPred  <- predict(ml,new=grillex,interval="prediction",level=0.95)     
       matlines(grillex,cbind(icConf,icPred[,-1]),lty=c(1,2,2,3,3),col=c("black","red","red","green","green"),lwd=2)     
          
     # un légende, c'est utile     
          
       legend(     
         x="topleft",     
         legend=c("valeurs originales","valeurs prédites","droite MCO","IC conf","IC pred"),     
         col=c("blue","red","black","red","green"),     
         pch=c(19,19,NA,NA,NA),     
         lty=c(0,0,1,3,3),     
         lwd=c(0,0,2,2,2)     
       ) # fin de legend     
          
     dev.off()     
          

     valeurs prédites     
     ================     
          
                  1         2     
     newx 100.00000 250.00000     
     newy  11.68237  24.18035     
          

          kmpreds

Il y a bien sûr ici une relation de causalité, la cause étant la distance et l'effet étant la consommation, ce que l'on peut résumer par plus on roule, plus on consomme.

 

3. Régression linéaire multiple

Effectuez une première RLM (Régression Linéaire Multiple de LN_NID en fonction de ALT, PENTE, HAUT et DIAM, puis une seconde en utilisant seulement ALT, PENTE et HAUT, puis une troisième avec ALT, PENTE et DIAM dans le dossier CHENILLES.

Y a-t-il une différence au niveau du signe des coefficients calculés ?

Quel est le meilleur de ces trois modèles ? Faut-il distinguer le meilleur modèle prédictif du meilleur modèle explicatif ?

Que peut-on en conclure sur l'influence de la colinéarité des variables dans une RLM ?

Comment demander à R de calculer le meilleur modèle pour LN_NID à partir des quatre variables ALT, PENTE, HAUT et DIAM ?

Solution :  

Voici le code R pour réaliser les régressions demandées :


     attach(chen)     
          
     # régression 1     
          
     cats("Régression numéro 1 ")     
          
     m1 <- lm( LN_NID ~ ALT + PENTE + HAUT + DIAM )     
     print(summary(m1))     
     fit1 <- m1$fitted.values     
          
     # régression 2     
          
     cats("Régression numéro 2 ")     
          
     m2 <- lm( LN_NID ~ ALT + PENTE + HAUT )     
     print(summary(m2))     
     fit2 <- m2$fitted.values     
          
     # régression 3     
          
     cats("Régression numéro 3 ")     
          
     m3 <- lm( LN_NID ~ ALT + PENTE + DIAM )     
     print(summary(m3))     
     fit3 <- m3$fitted.values     
          
     # tracés     
          
     gr("chen3modeles.png")     
          
     plot(LN_NID,LN_NID,pch=19,col="black")     
     abline(a=0,b=1)     
     points(LN_NID,fit1,pch=19,col="red",cex=0.7)     
     points(LN_NID,fit2,pch=19,col="green",cex=0.7)     
     points(LN_NID,fit3,pch=19,col="blue",cex=0.7)     
          
     legend(     
       x="topleft",     
       legend=c("points originaux","modèle 1","modèle 2","modèle 3"),     
       col=c("black","red","green","blue"),     
       pch=19     
     ) # fin de legend     
          
     dev.off()     
          
     detach(chen)     
          

Et les résultats correspondants :


          
     Régression numéro 1     
          
     =====================     
          
          
     Call:     
     lm(formula = LN_NID ~ ALT + PENTE + HAUT + DIAM)     
          
     Residuals:     
          Min       1Q   Median       3Q      Max     
     -2.02087 -0.25012  0.09002  0.35179  1.71056     
          
     Coefficients:     
                  Estimate Std. Error t value Pr(>|t|)     
     (Intercept)  7.732147   1.488585   5.194 1.63e-05 ***     
     ALT         -0.003924   0.001148  -3.419 0.001946 **     
     PENTE       -0.057343   0.019388  -2.958 0.006236 **     
     HAUT        -1.356138   0.319834  -4.240 0.000220 ***     
     DIAM         0.283059   0.076260   3.712 0.000905 ***     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 0.7906 on 28 degrees of freedom     
     Multiple R-squared:  0.6471,	Adjusted R-squared:  0.5967     
     F-statistic: 12.83 on 4 and 28 DF,  p-value: 4.677e-06     
          
          
     Régression numéro 2     
          
     =====================     
          
          
     Call:     
     lm(formula = LN_NID ~ ALT + PENTE + HAUT)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.0837 -0.5037 -0.0635  0.6903  1.5616     
          
     Coefficients:     
                  Estimate Std. Error t value Pr(>|t|)     
     (Intercept)  7.448568   1.784311   4.174 0.000249 ***     
     ALT         -0.003986   0.001377  -2.894 0.007153 **     
     PENTE       -0.058997   0.023264  -2.536 0.016860 *     
     HAUT        -0.293479   0.171126  -1.715 0.097016 .     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 0.949 on 29 degrees of freedom     
     Multiple R-squared:  0.4734,	Adjusted R-squared:  0.419     
     F-statistic: 8.691 on 3 and 29 DF,  p-value: 0.0002869     
          
          
     Régression numéro 3     
          
     =====================     
          
          
     Call:     
     lm(formula = LN_NID ~ ALT + PENTE + DIAM)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.2638 -0.7811  0.2125  0.7013  1.6313     
          
     Coefficients:     
                  Estimate Std. Error t value Pr(>|t|)     
     (Intercept)  7.237019   1.868584   3.873 0.000564 ***     
     ALT         -0.004659   0.001428  -3.262 0.002834 **     
     PENTE       -0.062853   0.024358  -2.580 0.015197 *     
     DIAM        -0.006389   0.042806  -0.149 0.882384     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 0.9955 on 29 degrees of freedom     
     Multiple R-squared:  0.4205,	Adjusted R-squared:  0.3605     
     F-statistic: 7.014 on 3 and 29 DF,  p-value: 0.001093     
          
          

Si on compare les modèles 1 et 3, on remarque le coefficient de DIAM est positif avec le modèle 1 (valeur 0.283) mais négatif avec le modèle 3 (valeur -0.006). Compte-tenu de la corrélation entre les varables (détaillée ci-dessous), le coefficient devrait être négatif. C'est sans doute la colinéarité entre HAUT et DIAM (r=0.905) qui cause cette "instabilité" de coefficient.

Le meilleur modèle prédictif de ces trois modèles est le premier, parce que son R2 ajusté est le plus élevé (0.5967 vs 0.4190 et 0.3605).

Malheureusement, ce modèle 1 est sans doute incorrect, à cause du signe de DIAM.

Le meilleur modèle explicatif est certainement le modèle 2.

          chen3modeles

L'étude de la corrélation entre les variables, y compris la "cible" LN_NID et la régression en fonction de chaque variable séparément peut s'obtenir par :


          
     # données     
          
     chen <- lit.dar("chenilles.dar")     
          
     # on retire NIDS (variable 11)     
          
     chen <- chen[,-11]     
          
     # régression de LN_NIDS en fonction chaque variable séparément     
          
     nbvar <- 1:(ncol(chen) - 1)     
     for (idv in nbvar) {     
      lavar <- names(chen)[idv]     
      cats(paste("Régression en fonction de",lavar,"seulement"))     
      fo <- paste("LN_NID ~",lavar)     
      ml <- lm( as.formula(fo),data=chen)     
      print(summary(ml))     
     } # fin pour idv     
          
     # matrice des corrélations     
          
     mdc(chen)     
          
     ## tracés de LN_NIDS en fonction chaque variable séparément     
          
     # avec la variable pairs() originale     
          
     gr("chenpairs1.png")     
     pairs(chen)     
     dev.off()     
          
     # avec notre fonction pairsi()     
          
     gr("chenpairs2.png")     
     pairsi(chen,opt=2)     
     dev.off()     
          

          
     Régression en fonction de ALT seulement     
     =======================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
          Min       1Q   Median       3Q      Max     
     -2.72074 -0.68616  0.00178  0.94230  1.74964     
          
     Coefficients:     
                  Estimate Std. Error t value Pr(>|t|)     
     (Intercept)  5.958414   1.936589   3.077  0.00435 **     
     ALT         -0.005148   0.001465  -3.513  0.00138 **     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 1.07 on 31 degrees of freedom     
     Multiple R-squared:  0.2847,	Adjusted R-squared:  0.2617     
     F-statistic: 12.34 on 1 and 31 DF,  p-value: 0.001384     
          
          
     Régression en fonction de PENTE seulement     
     =========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.2142 -0.9769  0.2790  0.9433  1.5437     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept)  1.31175    0.82689   1.586   0.1228     
     PENTE       -0.07320    0.02765  -2.648   0.0126 *     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 1.142 on 31 degrees of freedom     
     Multiple R-squared:  0.1844,	Adjusted R-squared:  0.1581     
     F-statistic:  7.01 on 1 and 31 DF,  p-value: 0.01263     
          
          
     Régression en fonction de PINS seulement     
     ========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.3184 -0.7781  0.2604  0.8380  1.8454     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept) -0.03901    0.29711  -0.131  0.89639     
     PINS        -0.06760    0.02006  -3.370  0.00203 **     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 1.082 on 31 degrees of freedom     
     Multiple R-squared:  0.2681,	Adjusted R-squared:  0.2445     
     F-statistic: 11.36 on 1 and 31 DF,  p-value: 0.002027     
          
          
     Régression en fonction de HAUT seulement     
     ========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.3125 -1.1291  0.2460  0.9391  1.8610     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept)   1.4513     0.8882   1.634   0.1124     
     HAUT         -0.5087     0.1944  -2.616   0.0136 *     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 1.145 on 31 degrees of freedom     
     Multiple R-squared:  0.1809,	Adjusted R-squared:  0.1545     
     F-statistic: 6.845 on 1 and 31 DF,  p-value: 0.01361     
          
          
     Régression en fonction de DIAM seulement     
     ========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.5451 -0.9392  0.2989  0.8394  1.7781     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept)  0.07346    0.80584   0.091    0.928     
     DIAM        -0.05814    0.05091  -1.142    0.262     
          
     Residual standard error: 1.239 on 31 degrees of freedom     
     Multiple R-squared:  0.04038,	Adjusted R-squared:  0.00942     
     F-statistic: 1.304 on 1 and 31 DF,  p-value: 0.2622     
          
          
     Régression en fonction de DENS seulement     
     ========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.2265 -0.8000  0.2660  0.8289  1.8401     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept)   0.8286     0.5095   1.626  0.11400     
     DENS         -0.9168     0.2647  -3.464  0.00158 **     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 1.074 on 31 degrees of freedom     
     Multiple R-squared:  0.2791,	Adjusted R-squared:  0.2558     
     F-statistic:    12 on 1 and 31 DF,  p-value: 0.001577     
          
          
     Régression en fonction de ORIEN seulement     
     =========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.3924 -0.8565  0.1141  0.9578  1.9767     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept)    1.719      1.939   0.887    0.382     
     ORIEN         -1.528      1.163  -1.314    0.199     
          
     Residual standard error: 1.231 on 31 degrees of freedom     
     Multiple R-squared:  0.05275,	Adjusted R-squared:  0.0222     
     F-statistic: 1.726 on 1 and 31 DF,  p-value: 0.1985     
          
          
     Régression en fonction de HAUTDOM seulement     
     ===========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -1.9251 -0.8830  0.4163  0.8720  1.5343     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept)  1.34730    0.63046   2.137  0.04060 *     
     HAUTDOM     -0.28657    0.07994  -3.585  0.00114 **     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 1.063 on 31 degrees of freedom     
     Multiple R-squared:  0.2931,	Adjusted R-squared:  0.2703     
     F-statistic: 12.85 on 1 and 31 DF,  p-value: 0.00114     
          
          
     Régression en fonction de STRAT seulement     
     =========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.0833 -0.8512  0.2844  0.9167  1.5247     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept)   1.7737     0.6537   2.713 0.010780 *     
     STRAT        -1.3054     0.3175  -4.111 0.000268 ***     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Residual standard error: 1.018 on 31 degrees of freedom     
     Multiple R-squared:  0.3528,	Adjusted R-squared:  0.3319     
     F-statistic:  16.9 on 1 and 31 DF,  p-value: 0.0002681     
          
          
     Régression en fonction de MELAN seulement     
     =========================================     
          
          
     Call:     
     lm(formula = as.formula(fo), data = chen)     
          
     Residuals:     
         Min      1Q  Median      3Q     Max     
     -2.6174 -0.8257  0.4276  0.9753  1.9877     
          
     Coefficients:     
                 Estimate Std. Error t value Pr(>|t|)     
     (Intercept)  -0.2787     1.5406  -0.181    0.858     
     MELAN        -0.3052     0.8706  -0.351    0.728     
          
     Residual standard error: 1.262 on 31 degrees of freedom     
     Multiple R-squared:  0.003949,	Adjusted R-squared:  -0.02818     
     F-statistic: 0.1229 on 1 and 31 DF,  p-value: 0.7283     
          
          
     Matrice des corrélations au sens de  Pearson pour  33  lignes et  11  colonnes     
          
                ALT  PENTE   PINS   HAUT   DIAM   DENS  ORIEN HAUTDOM  STRAT  MELAN LN_NID     
     ALT      1.000     
     PENTE    0.121  1.000     
     PINS     0.538  0.322  1.000     
     HAUT     0.321  0.137  0.414  1.000     
     DIAM     0.284  0.113  0.295  0.905  1.000     
     DENS     0.515  0.301  0.980  0.439  0.306  1.000     
     ORIEN    0.268 -0.152  0.128  0.058 -0.079  0.151  1.000     
     HAUTDOM  0.360  0.262  0.759  0.772  0.596  0.810  0.060   1.000     
     STRAT    0.364  0.326  0.877  0.460  0.267  0.909  0.063   0.854  1.000     
     MELAN   -0.100  0.129  0.206 -0.045 -0.025  0.130  0.138   0.054  0.175  1.000     
     LN_NID  -0.534 -0.429 -0.518 -0.425 -0.201 -0.528 -0.230  -0.541 -0.594 -0.063  1.000     
          

          chenpairs1          chenpairs2

Le dernier graphique montre une tendance négative pour la liaison entre DIAM et LN_NID, d'où notre remarque sur le modèle 1.

La fonction step() permet de déterminer le meilleur modèle par sélection ascendante et descendante de variables :


     # lecture des données     
          
     che <- lit.dar("chenilles.dar.txt")     
          
     # on se restreint aux variables indiquées     
          
     chen4v <- che[,c("ALT","PENTE","HAUT","DIAM","LN_NID")]     
          
     # modèle minimal     
          
     lmMin <- lm(LN_NID~1,data=chen4v)     
          
     # modèle maxnimal     
          
     lmMax <- lm(LN_NID~.,data=chen4v)     
          
     # sélection forward du meilleur modèle     
          
     cats("sélection ascendante")     
     step(lmMin,scope=list(upper=lmMax),trace=1,direction="forward")     
          
     # sélection backward du meilleur modèle     
          
     cats("sélection descendante")     
     step(lmMax,scope=list(lower=lmMin),trace=1,direction="backward")     
          
     # remarque : sélection sur l'ensemble des variables     
          
     catss("sélection globale")     
          
     lmMinG <- lm(LN_NID~1,data=che)     
     lmMaxG <- lm(LN_NID~.,data=che)     
     cats("sélection ascendante globale")     
     step(lmMinG,scope=list(upper=lmMaxG),trace=1,direction="forward")     
     cats("sélection descendante globale")     
     step(lmMaxG,scope=list(lower=lmMinG),trace=1,direction="backward")     
     cats("sélection descendante et ascendante")     
     step(lmMinG,scope=list(upper=lmMaxG),trace=1,direction="both")     
          
     # et bien sur sans la variable NIDS     
          
     catss("sélection globale sans NIDS")     
          
     cheSN  <- che[,-11]     
     lmMinG <- lm(LN_NID~1,data=cheSN)     
     lmMaxG <- lm(LN_NID~.,data=cheSN)     
     cats("sélection ascendante globale sans NIDS")     
     step(lmMinG,scope=list(upper=lmMaxG),trace=1,direction="forward")     
     cats("sélection descendante globale sans NIDS")     
     step(lmMaxG,scope=list(lower=lmMinG),trace=1,direction="backward")     
     cats("sélection descendante et ascendante sans NIDS")     
     step(lmMinG,scope=list(upper=lmMaxG),trace=1,direction="both")     
          

          
     sélection ascendante     
     ====================     
          
     Start:  AIC=15.44     
     LN_NID ~ 1     
          
             Df Sum of Sq    RSS     AIC     
     + ALT    1   14.1222 35.474  6.3856     
     + PENTE  1    9.1464 40.450 10.7172     
     + HAUT   1    8.9707 40.625 10.8602     
     <none>               49.596 15.4443     
     + DIAM   1    2.0025 47.594 16.0843     
          
     Step:  AIC=6.39     
     LN_NID ~ ALT     
          
             Df Sum of Sq    RSS    AIC     
     + PENTE  1    6.7095 28.764 1.4668     
     + HAUT   1    3.5668 31.907 4.8886     
     <none>               35.474 6.3856     
     + DIAM   1    0.1322 35.342 8.2623     
          
     Step:  AIC=1.47     
     LN_NID ~ ALT + PENTE     
          
            Df Sum of Sq    RSS    AIC     
     + HAUT  1   2.64866 26.116 0.2790     
     <none>              28.764 1.4668     
     + DIAM  1   0.02208 28.742 3.4414     
          
     Step:  AIC=0.28     
     LN_NID ~ ALT + PENTE + HAUT     
          
            Df Sum of Sq    RSS     AIC     
     + DIAM  1    8.6123 17.503 -10.926     
     <none>              26.116   0.279     
          
     Step:  AIC=-10.93     
     LN_NID ~ ALT + PENTE + HAUT + DIAM     
          
          
     Call:     
     lm(formula = LN_NID ~ ALT + PENTE + HAUT + DIAM, data = chen4v)     
          
     Coefficients:     
     (Intercept)          ALT        PENTE         HAUT         DIAM     
        7.732147    -0.003924    -0.057343    -1.356138     0.283059     
          
          
     sélection descendante     
     =====================     
          
     Start:  AIC=-10.93     
     LN_NID ~ ALT + PENTE + HAUT + DIAM     
          
             Df Sum of Sq    RSS      AIC     
     <none>               17.503 -10.9257     
     - PENTE  1    5.4683 22.972  -3.9540     
     - ALT    1    7.3067 24.810  -1.4135     
     - DIAM   1    8.6123 26.116   0.2790     
     - HAUT   1   11.2389 28.742   3.4414     
          
     Call:     
     lm(formula = LN_NID ~ ALT + PENTE + HAUT + DIAM, data = chen4v)     
          
     Coefficients:     
     (Intercept)          ALT        PENTE         HAUT         DIAM     
        7.732147    -0.003924    -0.057343    -1.356138     0.283059     
          
          
     + ================= +     
     +                   +     
     + sélection globale +     
     +                   +     
     + ================= +     
          
          
     sélection ascendante globale     
     ============================     
          
     Start:  AIC=15.44     
     LN_NID ~ 1     
          
               Df Sum of Sq    RSS      AIC     
     + NIDS     1    38.072 11.524 -30.7186     
     + STRAT    1    17.499 32.097   3.0848     
     + HAUTDOM  1    14.536 35.060   5.9986     
     + ALT      1    14.122 35.474   6.3856     
     + DENS     1    13.841 35.755   6.6459     
     + PINS     1    13.297 36.299   7.1445     
     + PENTE    1     9.146 40.450  10.7172     
     + HAUT     1     8.971 40.625  10.8602     
     <none>                 49.596  15.4443     
     + ORIEN    1     2.616 46.980  15.6559     
     + DIAM     1     2.002 47.594  16.0843     
     + MELAN    1     0.196 49.400  17.3138     
          
     Step:  AIC=-30.72     
     LN_NID ~ NIDS     
          
               Df Sum of Sq    RSS     AIC     
     + HAUT     1   0.70995 10.814 -30.817     
     <none>                 11.524 -30.719     
     + ALT      1   0.32905 11.195 -29.675     
     + HAUTDOM  1   0.24374 11.280 -29.424     
     + DIAM     1   0.19999 11.324 -29.296     
     + STRAT    1   0.11318 11.411 -29.044     
     + MELAN    1   0.10510 11.419 -29.021     
     + ORIEN    1   0.10125 11.423 -29.010     
     + DENS     1   0.06037 11.464 -28.892     
     + PENTE    1   0.05779 11.466 -28.884     
     + PINS     1   0.04093 11.483 -28.836     
          
     Step:  AIC=-30.82     
     LN_NID ~ NIDS + HAUT     
          
               Df Sum of Sq    RSS     AIC     
     + DIAM     1   0.69885 10.115 -31.021     
     <none>                 10.814 -30.817     
     + ALT      1   0.19360 10.620 -29.413     
     + ORIEN    1   0.11195 10.702 -29.160     
     + PENTE    1   0.07140 10.743 -29.035     
     + MELAN    1   0.05941 10.755 -28.999     
     + HAUTDOM  1   0.03580 10.778 -28.926     
     + STRAT    1   0.00474 10.809 -28.831     
     + PINS     1   0.00097 10.813 -28.820     
     + DENS     1   0.00018 10.814 -28.817     
          
     Step:  AIC=-31.02     
     LN_NID ~ NIDS + HAUT + DIAM     
          
               Df Sum of Sq     RSS     AIC     
     <none>                 10.1152 -31.021     
     + ALT      1   0.41371  9.7015 -30.400     
     + PENTE    1   0.19232  9.9229 -29.655     
     + HAUTDOM  1   0.15231  9.9629 -29.522     
     + MELAN    1   0.02261 10.0926 -29.095     
     + ORIEN    1   0.01665 10.0985 -29.076     
     + STRAT    1   0.01133 10.1038 -29.058     
     + DENS     1   0.00230 10.1129 -29.029     
     + PINS     1   0.00113 10.1140 -29.025     
          
     Call:     
     lm(formula = LN_NID ~ NIDS + HAUT + DIAM, data = che)     
          
     Coefficients:     
     (Intercept)         NIDS         HAUT         DIAM     
        -0.86290      1.19199     -0.50992      0.08869     
          
          
     sélection descendante globale     
     =============================     
          
     Start:  AIC=-19.42     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + DENS + ORIEN + HAUTDOM +     
         STRAT + MELAN + NIDS     
          
               Df Sum of Sq     RSS      AIC     
     - ORIEN    1    0.0001  8.8527 -21.4210     
     - MELAN    1    0.0155  8.8681 -21.3637     
     - DENS     1    0.0518  8.9043 -21.2292     
     - PINS     1    0.1159  8.9684 -20.9923     
     - STRAT    1    0.1525  9.0050 -20.8580     
     - HAUTDOM  1    0.3073  9.1598 -20.2954     
     - PENTE    1    0.4532  9.3057 -19.7741     
     <none>                  8.8525 -19.4216     
     - ALT      1    0.5751  9.4276 -19.3445     
     - DIAM     1    0.9891  9.8417 -17.9262     
     - HAUT     1    1.4373 10.2898 -16.4566     
     - NIDS     1    6.2773 15.1299  -3.7346     
          
     Step:  AIC=-21.42     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + DENS + HAUTDOM +     
         STRAT + MELAN + NIDS     
          
               Df Sum of Sq     RSS      AIC     
     - MELAN    1    0.0168  8.8695 -23.3583     
     - DENS     1    0.0589  8.9116 -23.2021     
     - PINS     1    0.1264  8.9791 -22.9531     
     - STRAT    1    0.1661  9.0188 -22.8076     
     - HAUTDOM  1    0.3145  9.1672 -22.2690     
     - PENTE    1    0.4696  9.3223 -21.7154     
     <none>                  8.8527 -21.4210     
     - ALT      1    0.5960  9.4487 -21.2707     
     - DIAM     1    1.1896 10.0423 -19.2602     
     - HAUT     1    1.6953 10.5480 -17.6389     
     - NIDS     1    6.3155 15.1682  -5.6511     
          
     Step:  AIC=-23.36     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + DENS + HAUTDOM +     
         STRAT + NIDS     
          
               Df Sum of Sq     RSS      AIC     
     - DENS     1    0.0436  8.9131 -25.1965     
     - PINS     1    0.1106  8.9801 -24.9492     
     - STRAT    1    0.1788  9.0483 -24.6996     
     - HAUTDOM  1    0.3103  9.1798 -24.2236     
     - PENTE    1    0.4568  9.3264 -23.7009     
     <none>                  8.8695 -23.3583     
     - ALT      1    0.5895  9.4591 -23.2346     
     - DIAM     1    1.1807 10.0502 -21.2342     
     - HAUT     1    1.6933 10.5629 -19.5924     
     - NIDS     1    6.6786 15.5481  -6.8348     
          
     Step:  AIC=-25.2     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + HAUTDOM + STRAT +     
         NIDS     
          
               Df Sum of Sq     RSS      AIC     
     - PINS     1    0.1018  9.0149 -26.8217     
     - STRAT    1    0.2258  9.1389 -26.3708     
     - HAUTDOM  1    0.2672  9.1803 -26.2219     
     - PENTE    1    0.4322  9.3453 -25.6338     
     <none>                  8.9131 -25.1965     
     - ALT      1    0.6467  9.5598 -24.8849     
     - DIAM     1    1.1528 10.0659 -23.1828     
     - HAUT     1    1.6602 10.5733 -21.5597     
     - NIDS     1    6.6380 15.5511  -8.8285     
          
     Step:  AIC=-26.82     
     LN_NID ~ ALT + PENTE + HAUT + DIAM + HAUTDOM + STRAT + NIDS     
          
               Df Sum of Sq     RSS      AIC     
     - STRAT    1    0.1252  9.1401 -28.3665     
     - HAUTDOM  1    0.3019  9.3168 -27.7348     
     - PENTE    1    0.3788  9.3937 -27.4634     
     - ALT      1    0.5590  9.5739 -26.8364     
     <none>                  9.0149 -26.8217     
     - DIAM     1    1.2879 10.3028 -24.4150     
     - HAUT     1    1.8765 10.8914 -22.5815     
     - NIDS     1    7.3061 16.3210  -9.2338     
          
     Step:  AIC=-28.37     
     LN_NID ~ ALT + PENTE + HAUT + DIAM + HAUTDOM + NIDS     
          
               Df Sum of Sq     RSS      AIC     
     - HAUTDOM  1    0.2170  9.3572 -29.5921     
     - PENTE    1    0.4028  9.5430 -28.9432     
     <none>                  9.1401 -28.3665     
     - ALT      1    0.5900  9.7301 -28.3023     
     - DIAM     1    1.3283 10.4685 -25.8886     
     - HAUT     1    1.7550 10.8951 -24.5703     
     - NIDS     1    8.3282 17.4683  -8.9919     
          
     Step:  AIC=-29.59     
     LN_NID ~ ALT + PENTE + HAUT + DIAM + NIDS     
          
             Df Sum of Sq     RSS     AIC     
     - PENTE  1    0.3443  9.7015 -30.400     
     - ALT    1    0.5657  9.9229 -29.655     
     <none>                9.3572 -29.592     
     - DIAM   1    1.1472 10.5044 -27.776     
     - HAUT   1    1.6400 10.9972 -26.263     
     - NIDS   1    8.1462 17.5034 -10.926     
          
     Step:  AIC=-30.4     
     LN_NID ~ ALT + HAUT + DIAM + NIDS     
          
            Df Sum of Sq     RSS     AIC     
     - ALT   1    0.4137 10.1152 -31.021     
     <none>               9.7015 -30.400     
     - DIAM  1    0.9190 10.6204 -29.413     
     - HAUT  1    1.3915 11.0929 -27.977     
     - NIDS  1   13.2703 22.9717  -3.954     
          
     Step:  AIC=-31.02     
     LN_NID ~ HAUT + DIAM + NIDS     
          
            Df Sum of Sq    RSS      AIC     
     <none>              10.115 -31.0215     
     - DIAM  1    0.6989 10.814 -30.8169     
     - HAUT  1    1.2088 11.324 -29.2963     
     - NIDS  1   21.2832 31.398   4.3582     
          
     Call:     
     lm(formula = LN_NID ~ HAUT + DIAM + NIDS, data = che)     
          
     Coefficients:     
     (Intercept)         HAUT         DIAM         NIDS     
        -0.86290     -0.50992      0.08869      1.19199     
          
          
     sélection descendante et ascendante     
     ===================================     
          
     Start:  AIC=15.44     
     LN_NID ~ 1     
          
               Df Sum of Sq    RSS      AIC     
     + NIDS     1    38.072 11.524 -30.7186     
     + STRAT    1    17.499 32.097   3.0848     
     + HAUTDOM  1    14.536 35.060   5.9986     
     + ALT      1    14.122 35.474   6.3856     
     + DENS     1    13.841 35.755   6.6459     
     + PINS     1    13.297 36.299   7.1445     
     + PENTE    1     9.146 40.450  10.7172     
     + HAUT     1     8.971 40.625  10.8602     
     <none>                 49.596  15.4443     
     + ORIEN    1     2.616 46.980  15.6559     
     + DIAM     1     2.002 47.594  16.0843     
     + MELAN    1     0.196 49.400  17.3138     
          
     Step:  AIC=-30.72     
     LN_NID ~ NIDS     
          
               Df Sum of Sq    RSS     AIC     
     + HAUT     1     0.710 10.814 -30.817     
     <none>                 11.524 -30.719     
     + ALT      1     0.329 11.195 -29.675     
     + HAUTDOM  1     0.244 11.280 -29.424     
     + DIAM     1     0.200 11.324 -29.296     
     + STRAT    1     0.113 11.411 -29.044     
     + MELAN    1     0.105 11.419 -29.021     
     + ORIEN    1     0.101 11.423 -29.010     
     + DENS     1     0.060 11.464 -28.892     
     + PENTE    1     0.058 11.466 -28.884     
     + PINS     1     0.041 11.483 -28.836     
     - NIDS     1    38.072 49.596  15.444     
          
     Step:  AIC=-30.82     
     LN_NID ~ NIDS + HAUT     
          
               Df Sum of Sq    RSS     AIC     
     + DIAM     1    0.6989 10.115 -31.021     
     <none>                 10.814 -30.817     
     - HAUT     1    0.7099 11.524 -30.719     
     + ALT      1    0.1936 10.620 -29.413     
     + ORIEN    1    0.1119 10.702 -29.160     
     + PENTE    1    0.0714 10.743 -29.035     
     + MELAN    1    0.0594 10.755 -28.999     
     + HAUTDOM  1    0.0358 10.778 -28.926     
     + STRAT    1    0.0047 10.809 -28.831     
     + PINS     1    0.0010 10.813 -28.820     
     + DENS     1    0.0002 10.814 -28.817     
     - NIDS     1   29.8113 40.625  10.860     
          
     Step:  AIC=-31.02     
     LN_NID ~ NIDS + HAUT + DIAM     
          
               Df Sum of Sq     RSS      AIC     
     <none>                 10.1152 -31.0215     
     - DIAM     1    0.6989 10.8140 -30.8169     
     + ALT      1    0.4137  9.7015 -30.3996     
     + PENTE    1    0.1923  9.9229 -29.6550     
     + HAUTDOM  1    0.1523  9.9629 -29.5222     
     - HAUT     1    1.2088 11.3240 -29.2963     
     + MELAN    1    0.0226 10.0926 -29.0954     
     + ORIEN    1    0.0167 10.0985 -29.0759     
     + STRAT    1    0.0113 10.1038 -29.0585     
     + DENS     1    0.0023 10.1129 -29.0291     
     + PINS     1    0.0011 10.1140 -29.0252     
     - NIDS     1   21.2832 31.3984   4.3582     
          
     Call:     
     lm(formula = LN_NID ~ NIDS + HAUT + DIAM, data = che)     
          
     Coefficients:     
     (Intercept)         NIDS         HAUT         DIAM     
        -0.86290      1.19199     -0.50992      0.08869     
          
          
     + =========================== +     
     +                             +     
     + sélection globale sans NIDS +     
     +                             +     
     + =========================== +     
          
          
     sélection ascendante globale sans NIDS     
     ======================================     
          
     Start:  AIC=15.44     
     LN_NID ~ 1     
          
               Df Sum of Sq    RSS     AIC     
     + STRAT    1   17.4987 32.097  3.0848     
     + HAUTDOM  1   14.5357 35.060  5.9986     
     + ALT      1   14.1222 35.474  6.3856     
     + DENS     1   13.8412 35.755  6.6459     
     + PINS     1   13.2969 36.299  7.1445     
     + PENTE    1    9.1464 40.450 10.7172     
     + HAUT     1    8.9707 40.625 10.8602     
     <none>                 49.596 15.4443     
     + ORIEN    1    2.6164 46.980 15.6559     
     + DIAM     1    2.0025 47.594 16.0843     
     + MELAN    1    0.1958 49.400 17.3138     
          
     Step:  AIC=3.08     
     LN_NID ~ STRAT     
          
               Df Sum of Sq    RSS     AIC     
     + ALT      1    5.7642 26.333 -1.4474     
     + PENTE    1    3.0899 29.008  1.7446     
     <none>                 32.097  3.0848     
     + ORIEN    1    1.8378 30.260  3.1390     
     + HAUT     1    1.4585 30.639  3.5502     
     + HAUTDOM  1    0.2153 31.882  4.8627     
     + DIAM     1    0.0945 32.003  4.9875     
     + MELAN    1    0.0853 32.012  4.9971     
     + DENS     1    0.0368 32.061  5.0470     
     + PINS     1    0.0020 32.095  5.0828     
          
     Step:  AIC=-1.45     
     LN_NID ~ STRAT + ALT     
          
               Df Sum of Sq    RSS     AIC     
     + PENTE    1   3.07012 23.263 -3.5382     
     + DENS     1   2.27340 24.060 -2.4269     
     + PINS     1   1.94158 24.392 -1.9749     
     <none>                 26.333 -1.4474     
     + HAUT     1   0.60003 25.733 -0.2080     
     + ORIEN    1   0.55986 25.773 -0.1565     
     + HAUTDOM  1   0.04815 26.285  0.4922     
     + DIAM     1   0.03828 26.295  0.5046     
     + MELAN    1   0.01903 26.314  0.5288     
          
     Step:  AIC=-3.54     
     LN_NID ~ STRAT + ALT + PENTE     
          
               Df Sum of Sq    RSS     AIC     
     + PINS     1   2.42797 20.835 -5.1757     
     + DENS     1   2.34020 20.923 -5.0370     
     <none>                 23.263 -3.5382     
     + ORIEN    1   1.21501 22.048 -3.3084     
     + HAUT     1   0.64494 22.618 -2.4660     
     + HAUTDOM  1   0.07706 23.186 -1.6477     
     + DIAM     1   0.06082 23.202 -1.6246     
     + MELAN    1   0.00000 23.263 -1.5382     
          
     Step:  AIC=-5.18     
     LN_NID ~ STRAT + ALT + PENTE + PINS     
          
               Df Sum of Sq    RSS     AIC     
     + ORIEN    1   1.38056 19.454 -5.4381     
     <none>                 20.835 -5.1757     
     + HAUT     1   0.47688 20.358 -3.9398     
     + MELAN    1   0.13013 20.705 -3.3824     
     + DENS     1   0.10288 20.732 -3.3390     
     + HAUTDOM  1   0.07171 20.763 -3.2895     
     + DIAM     1   0.03910 20.796 -3.2377     
          
     Step:  AIC=-5.44     
     LN_NID ~ STRAT + ALT + PENTE + PINS + ORIEN     
          
               Df Sum of Sq    RSS     AIC     
     <none>                 19.454 -5.4381     
     + HAUT     1   0.50574 18.949 -4.3073     
     + DENS     1   0.28304 19.172 -3.9218     
     + HAUTDOM  1   0.08666 19.368 -3.5855     
     + MELAN    1   0.01632 19.438 -3.4658     
     + DIAM     1   0.00004 19.454 -3.4382     
          
     Call:     
     lm(formula = LN_NID ~ STRAT + ALT + PENTE + PINS + ORIEN, data = cheSN)     
          
     Coefficients:     
     (Intercept)        STRAT          ALT        PENTE         PINS        ORIEN     
       11.107720    -1.745234    -0.004373    -0.054446     0.071498    -1.175855     
          
          
     sélection descendante globale sans NIDS     
     =======================================     
          
     Start:  AIC=-3.73     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + DENS + ORIEN + HAUTDOM +     
         STRAT + MELAN     
          
               Df Sum of Sq    RSS     AIC     
     - DENS     1    0.0357 15.166 -5.6569     
     - ORIEN    1    0.0383 15.168 -5.6511     
     - MELAN    1    0.2389 15.369 -5.2177     
     - PINS     1    0.3208 15.451 -5.0422     
     - HAUTDOM  1    0.4023 15.532 -4.8685     
     <none>                 15.130 -3.7346     
     - STRAT    1    1.5190 16.649 -2.5774     
     - DIAM     1    3.3869 18.517  0.9317     
     - HAUT     1    3.6205 18.750  1.3453     
     - PENTE    1    4.1551 19.285  2.2730     
     - ALT      1    5.5717 20.702  4.6122     
          
     Step:  AIC=-5.66     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + ORIEN + HAUTDOM +     
         STRAT + MELAN     
          
               Df Sum of Sq    RSS     AIC     
     - ORIEN    1    0.0871 15.253 -7.4679     
     - MELAN    1    0.2063 15.372 -7.2110     
     - HAUTDOM  1    0.3742 15.540 -6.8525     
     <none>                 15.166 -5.6569     
     - PINS     1    0.9909 16.156 -5.5683     
     - STRAT    1    1.8513 17.017 -3.8560     
     - DIAM     1    3.5100 18.676 -0.7866     
     - HAUT     1    4.0339 19.199  0.1263     
     - PENTE    1    4.1218 19.287  0.2771     
     - ALT      1    5.5373 20.703  2.6142     
          
     Step:  AIC=-7.47     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + HAUTDOM + STRAT +     
         MELAN     
          
               Df Sum of Sq    RSS     AIC     
     - MELAN    1    0.2984 15.551 -8.8285     
     - HAUTDOM  1    0.3782 15.631 -8.6597     
     - PINS     1    0.9455 16.198 -7.4832     
     <none>                 15.253 -7.4679     
     - STRAT    1    1.7752 17.028 -5.8346     
     - PENTE    1    4.0447 19.297 -1.7058     
     - DIAM     1    4.7047 19.957 -0.5960     
     - HAUT     1    4.9416 20.194 -0.2066     
     - ALT      1    6.3926 21.645  2.0831     
          
     Step:  AIC=-8.83     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + HAUTDOM + STRAT     
          
               Df Sum of Sq    RSS     AIC     
     - HAUTDOM  1    0.4992 16.050 -9.7858     
     - PINS     1    0.7699 16.321 -9.2338     
     <none>                 15.551 -8.8285     
     - STRAT    1    1.9116 17.463 -7.0027     
     - PENTE    1    4.1556 19.707 -3.0132     
     - DIAM     1    4.5302 20.081 -2.3918     
     - HAUT     1    4.8757 20.427 -1.8288     
     - ALT      1    6.1082 21.659  0.1045     
          
     Step:  AIC=-9.79     
     LN_NID ~ ALT + PENTE + PINS + HAUT + DIAM + STRAT     
          
             Df Sum of Sq    RSS     AIC     
     - PINS   1    0.9272 16.977 -9.9325     
     <none>               16.050 -9.7858     
     - STRAT  1    1.4451 17.495 -8.9408     
     - PENTE  1    4.2645 20.315 -4.0102     
     - DIAM   1    4.3079 20.358 -3.9398     
     - HAUT   1    4.7457 20.796 -3.2377     
     - ALT    1    6.6711 22.721 -0.3156     
          
     Step:  AIC=-9.93     
     LN_NID ~ ALT + PENTE + HAUT + DIAM + STRAT     
          
             Df Sum of Sq    RSS      AIC     
     - STRAT  1    0.5259 17.503 -10.9257     
     <none>               16.977  -9.9325     
     - PENTE  1    4.0243 21.002  -4.9127     
     - DIAM   1    5.6406 22.618  -2.4660     
     - ALT    1    5.8582 22.836  -2.1500     
     - HAUT   1    6.2248 23.202  -1.6246     
          
     Step:  AIC=-10.93     
     LN_NID ~ ALT + PENTE + HAUT + DIAM     
          
             Df Sum of Sq    RSS      AIC     
     <none>               17.503 -10.9257     
     - PENTE  1    5.4683 22.972  -3.9540     
     - ALT    1    7.3067 24.810  -1.4135     
     - DIAM   1    8.6123 26.116   0.2790     
     - HAUT   1   11.2389 28.742   3.4414     
          
     Call:     
     lm(formula = LN_NID ~ ALT + PENTE + HAUT + DIAM, data = cheSN)     
          
     Coefficients:     
     (Intercept)          ALT        PENTE         HAUT         DIAM     
        7.732147    -0.003924    -0.057343    -1.356138     0.283059     
          
          
     sélection descendante et ascendante sans NIDS     
     =============================================     
          
     Start:  AIC=15.44     
     LN_NID ~ 1     
          
               Df Sum of Sq    RSS     AIC     
     + STRAT    1   17.4987 32.097  3.0848     
     + HAUTDOM  1   14.5357 35.060  5.9986     
     + ALT      1   14.1222 35.474  6.3856     
     + DENS     1   13.8412 35.755  6.6459     
     + PINS     1   13.2969 36.299  7.1445     
     + PENTE    1    9.1464 40.450 10.7172     
     + HAUT     1    8.9707 40.625 10.8602     
     <none>                 49.596 15.4443     
     + ORIEN    1    2.6164 46.980 15.6559     
     + DIAM     1    2.0025 47.594 16.0843     
     + MELAN    1    0.1958 49.400 17.3138     
          
     Step:  AIC=3.08     
     LN_NID ~ STRAT     
          
               Df Sum of Sq    RSS     AIC     
     + ALT      1    5.7642 26.333 -1.4474     
     + PENTE    1    3.0899 29.008  1.7446     
     <none>                 32.097  3.0848     
     + ORIEN    1    1.8378 30.260  3.1390     
     + HAUT     1    1.4585 30.639  3.5502     
     + HAUTDOM  1    0.2153 31.882  4.8627     
     + DIAM     1    0.0945 32.003  4.9875     
     + MELAN    1    0.0853 32.012  4.9971     
     + DENS     1    0.0368 32.061  5.0470     
     + PINS     1    0.0020 32.095  5.0828     
     - STRAT    1   17.4987 49.596 15.4443     
          
     Step:  AIC=-1.45     
     LN_NID ~ STRAT + ALT     
          
               Df Sum of Sq    RSS     AIC     
     + PENTE    1    3.0701 23.263 -3.5382     
     + DENS     1    2.2734 24.060 -2.4269     
     + PINS     1    1.9416 24.392 -1.9749     
     <none>                 26.333 -1.4474     
     + HAUT     1    0.6000 25.733 -0.2080     
     + ORIEN    1    0.5599 25.773 -0.1565     
     + HAUTDOM  1    0.0482 26.285  0.4922     
     + DIAM     1    0.0383 26.295  0.5046     
     + MELAN    1    0.0190 26.314  0.5288     
     - ALT      1    5.7642 32.097  3.0848     
     - STRAT    1    9.1407 35.474  6.3856     
          
     Step:  AIC=-3.54     
     LN_NID ~ STRAT + ALT + PENTE     
          
               Df Sum of Sq    RSS     AIC     
     + PINS     1    2.4280 20.835 -5.1757     
     + DENS     1    2.3402 20.923 -5.0370     
     <none>                 23.263 -3.5382     
     + ORIEN    1    1.2150 22.048 -3.3084     
     + HAUT     1    0.6449 22.618 -2.4660     
     + HAUTDOM  1    0.0771 23.186 -1.6477     
     + DIAM     1    0.0608 23.202 -1.6246     
     + MELAN    1    0.0000 23.263 -1.5382     
     - PENTE    1    3.0701 26.333 -1.4474     
     - STRAT    1    5.5013 28.764  1.4668     
     - ALT      1    5.7445 29.008  1.7446     
          
     Step:  AIC=-5.18     
     LN_NID ~ STRAT + ALT + PENTE + PINS     
          
               Df Sum of Sq    RSS     AIC     
     + ORIEN    1    1.3806 19.454 -5.4381     
     <none>                 20.835 -5.1757     
     + HAUT     1    0.4769 20.358 -3.9398     
     - PINS     1    2.4280 23.263 -3.5382     
     + MELAN    1    0.1301 20.705 -3.3824     
     + DENS     1    0.1029 20.732 -3.3390     
     + HAUTDOM  1    0.0717 20.763 -3.2895     
     + DIAM     1    0.0391 20.796 -3.2377     
     - PENTE    1    3.5565 24.392 -1.9749     
     - STRAT    1    6.5118 27.347  1.7992     
     - ALT      1    8.1381 28.973  3.7054     
          
     Step:  AIC=-5.44     
     LN_NID ~ STRAT + ALT + PENTE + PINS + ORIEN     
          
               Df Sum of Sq    RSS     AIC     
     <none>                 19.454 -5.4381     
     - ORIEN    1    1.3806 20.835 -5.1757     
     + HAUT     1    0.5057 18.949 -4.3073     
     + DENS     1    0.2830 19.172 -3.9218     
     + HAUTDOM  1    0.0867 19.368 -3.5855     
     + MELAN    1    0.0163 19.438 -3.4658     
     + DIAM     1    0.0000 19.454 -3.4382     
     - PINS     1    2.5935 22.048 -3.3084     
     - PENTE    1    4.3177 23.772 -0.8237     
     - ALT      1    6.4077 25.862  1.9570     
     - STRAT    1    6.6742 26.129  2.2954     
          
     Call:     
     lm(formula = LN_NID ~ STRAT + ALT + PENTE + PINS + ORIEN, data = cheSN)     
          
     Coefficients:     
     (Intercept)        STRAT          ALT        PENTE         PINS        ORIEN     
       11.107720    -1.745234    -0.004373    -0.054446     0.071498    -1.175855     
          
          

 

4. Régression logistique simple

Réaliser la régression de la variable binaire de cardiopathie coronarienne (CHD69) à l'aide de la seule variable AGE dans le dossier WCGS.

Quel est le rapport de côte pour une augmentation de l'age d'un an ? Et pour 10 ans ?

Quelle est la probabilité d'avoir l'évènement pour une personne de 55 ans ?

Données :   wcgs.dar.

Solution :  

Commençons par regarder les variables CHD69 et AGE séparément puis conjointement.


          
     # lecture des données     
          
     WCGS       <- lit.dar("wcgs.dar")     
     age        <- WCGS$AGE     
          
     # recodage de chd     
          
     chd        <- factor(WCGS$CHD69,labels=c("Non","Oui"))     
     chd69      <- ifelse(chd=="Non",0,1)     
          
     # statistiques élémentaires     
          
     decritQT("AGE dans WCGS",age,"ans",TRUE,"wcgs_age.png")     
     decritQL("CHD dans WCGS",chd69,c("Non","Oui"),TRUE,"wcgs_chd.png")     
          
     decritQTparFacteur(     
       titreQT="AGE dans WCGS",     
       nomVarQT=age,     
       unite="ans",     
       titreQL="CHD",     
       nomVarQL=chd69,     
       labels=c("Non","Oui"),     
       graphique=TRUE,     
       fichier_image="wcgs1.png"     
     ) # fin de decritQTparFacteur     
          
          
          
          

          
     DESCRIPTION STATISTIQUE DE LA VARIABLE  AGE dans WCGS     
          
      Taille              3154        individus     
      Moyenne               46.2787         ans     
      Ecart-type             5.5240         ans     
      Coef. de variation    12               %     
      1er Quartile          42.0000         ans     
      Mediane               45.0000         ans     
      3eme Quartile         50.0000         ans     
      iqr absolu             8.0000         ans     
      iqr relatif           18.0000          %     
      Minimum               39              ans     
      Maximum               59              ans     
          
      Tracé tige et feuilles     
          
       The decimal point is at the |     
          
       39 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+136     
       40 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+147     
       41 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+103     
       42 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+92     
       43 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+85     
       44 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+105     
       45 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+56     
       46 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+40     
       47 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+17     
       48 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+34     
       49 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+4     
       50 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+5     
       51 | 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000     
       52 | 00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000     
       53 | 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000     
       54 | 00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000     
       55 | 00000000000000000000000000000000000000000000000000000000000000000000000000000000     
       56 | 0000000000000000000000000000000000000000000000000000000000000000000000000000     
       57 | 000000000000000000000000000000000000000000000000000000000000000     
       58 | 00000000000000000000000000000000000000000000000000000000     
       59 | 00000000000000000000000000000000000000000000000     
          
          
      vous pouvez utiliser  wcgs_age.png     
          
     TRI A PLAT DE :  CHD dans WCGS (ordre des modalités)     
          
                        Non  Oui  Total     
      Effectif         2897  257   3154     
      Cumul Effectif   2897 3154   3154     
      Frequence (en %)   92    8    100     
      Cumul fréquences   92  100    100     
          
     VARIABLE :  CHD dans WCGS (par fréquence décroissante)     
          
                         Non  Oui  Total     
       Effectif         2897  257   3154     
       Cumul Effectif   2897 3154   3154     
       Frequence (en %)   92    8    100     
       Cumul fréquences   92  100    100     
          
     [1] "Non" "Oui"     
          
      vous pouvez utiliser  wcgs_chd.png     
          
          
          
     VARIABLE QT  AGE dans WCGS ,unit=ans     
     VARIABLE QL  CHD  labels :  Non Oui     
          
                 N       Moy Unite       Ect Cdv        Q1       Med        Q3       EIQ       Min     
     Global   3154    46.279   ans     5.523  12    42.000    45.000    50.000     8.000    39.000     
     Non      2897    46.082   ans     5.456  12    41.000    45.000    50.000     9.000    39.000     
     Oui       257    48.490   ans     5.790  12    44.000    49.000    53.000     9.000    39.000     
                  Max     
     Global    59.000     
     Non       59.000     
     Oui       59.000     
          
     Analysis of Variance Table     
          
     Response: nomVarQT     
                 Df Sum Sq Mean Sq F value    Pr(>F)     
     ql           1   1369 1368.52   45.48 1.827e-11 ***     
     Residuals 3152  94846   30.09     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     Kruskal-Wallis rank sum test     
          
     data:  nomVarQT by ql     
     Kruskal-Wallis chi-squared = 40.8787, df = 1, p-value = 1.62e-10     
          
          
wcgs_chd.png

          

wcgs_age.png

          wcgs1

On peut en conclure que pour les données présentées, l'AGE est assez homogène (cdv=12 %), peu étendu et non normal, et que la CHD est assez peu présente (8 %), avec une différence significative selon l'age (p<<0,0001).

Réalisons maintenant la régression logistique de la variable CHD69 (binaire) à l'aide de la seule variable AGE (continue) à l'aide de la fonction glm et dont les résultats peuvent être consultés via la fonction summary.glm. Comme il est d'usage de rapporter la log-vraisemblance du modèle, les intervalles de confiance des coefficients, le rapport de vraisemblance (LR) et le R2, nous utilisons, en plus des fonctions confint et logLik du package élémentaire stats les fonctions confint.glm et lrm des packages MASS et Design qui fournissent directement ces valeurs. Remarque : la fonction classique anova fournit aussi le LR (qu'elle nomme Deviance).


     # lecture des données     
          
     WCGS       <- lit.dar("wcgs.dar")     
     age        <- WCGS$AGE     
     chd        <- factor(WCGS$CHD69,labels=c("Non","Oui"))     
          
     # régression logistique     
          
     rlb <- glm( formula = chd ~ age, family = binomial() )     
     summary.glm( rlb )     
     anova( rlb )     
          
     # il manque la log-vraisemblance et les intervalles de confiance     
          
     logLik( rlb )     
          
     library(MASS)     
          
     confint( rlb )     
          
     # l'intervalle calculé par Stata dans le livre de Vittinghoff correspond     
     # à l'intervalle sous hypothèse de normalité asymptotique     
          
     confint.default( rlb )     
          
     library(Design)     
     lrm( rlb )     
          
          

          
     > rlb <- glm( formula = chd ~ age, family = binomial() )     
          
     > summary.glm( rlb )     
          
     Call:     
     glm(formula = chd ~ age, family = binomial())     
          
     Deviance Residuals:     
         Min       1Q   Median       3Q      Max     
     -0.6208  -0.4545  -0.3668  -0.3292   2.4835     
          
     Coefficients:     
                 Estimate Std. Error z value Pr(>|z|)     
     (Intercept) -5.93952    0.54932 -10.813  < 2e-16 ***     
     age          0.07442    0.01130   6.585 4.56e-11 ***     
     ---     
     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
          
     (Dispersion parameter for binomial family taken to be 1)     
          
         Null deviance: 1781.2  on 3153  degrees of freedom     
     Residual deviance: 1738.4  on 3152  degrees of freedom     
     AIC: 1742.4     
          
     Number of Fisher Scoring iterations: 5     
          
     > anova( rlb )     
          
     Analysis of Deviance Table     
          
     Model: binomial, link: logit     
          
     Response: chd     
          
     Terms added sequentially (first to last)     
          
          
          Df Deviance Resid. Df Resid. Dev     
     NULL                  3153     1781.2     
     age   1   42.888      3152     1738.4     
          
     > logLik( rlb )     
          
     'log Lik.' -869.178 (df=2)     
          
          
     > library(MASS)     
     > confint( rlb )     
                       2.5 %      97.5 %     
     (Intercept) -7.02449218 -4.86960992     
     age          0.05227351  0.09661197     
          
     > confint.default( rlb )     
          
                       2.5 %      97.5 %     
     (Intercept) -7.01616019 -4.86287169     
     age          0.05227038  0.09657473     
          
     > library(Design)     
     > lrm( rlb )     
          
     Logistic Regression Model     
          
     lrm(formula = rlb)     
          
     Frequencies of Responses     
      Non  Oui     
     2897  257     
          
     Obs   Max Deriv Model L.R.      d.f.        P        C      Dxy    Gamma    Tau-a       R2      Brier     
     3154      9e-08      42.89         1        0     0.62     0.24    0.252    0.036    0.031      0.074     
          
               Coef     S.E.    Wald Z P     
     Intercept -5.93952 0.54932 -10.81 0     
     age        0.07442 0.01130   6.58 0     
          
          
          

Le rapport de côte p/(1-p) pour une augmentation de l'age d'un an s'obtient en prenant l'exponentielle du coefficient β1, soit ici exp(0.074*1) qui vaut à peu près 1,077 ; pour 10 ans, c'est exp(0.074*10) soit 2,105.

La probabilité d'avoir un trouble coronarien pour une personne de 55 ans se calcule avec la fonction logistique :

p(55) = exp(-5.940+0.074*55)/(1 + exp(-5.940+0.074*55)) = 0.1335417

 

5. Régression logistique multiple

Réalisez la régression de la variable de cardiopathie coronarienne (CHD69) à l'aide des 5 variables age, chol, sbp, bmi et smoke dans le dossier WCGS.

Quel est le meilleur modèle ?

Solution :  

Voici, sans explication superflue, le code et les résultats :


     # lecture des données     
          
     WCGS       <- lit.dar("wcgs.dar")     
     lesv       <- c("AGE","CHOL","SBP","BMI","SMOKE")     
     wcgs       <- WCGS[,lesv]     
          
     # définition de la matrice des résultats     
          
     nbv <- ncol(wcgs)     
     mdr <- data.frame(matrix(nrow=nbv,ncol=2))     
     names(mdr)     <- c("AUROC","AIC")     
     row.names(mdr) <- names(wcgs)     
          
     # on utilise une boucle pour stocker les AUROC     
          
     library(limma) # pour la fonction auROC     
          
     for (idv in (1:nbv)) {     
       rlb <- glm(WCGS$CHD69 ~ wcgs[,idv])     
       re  <- summary(ml)     
       ypred <- predict(rlb,newdata=wcgs,type="response")     
       mdr[idv,1] <- auROC(WCGS$CHD69,ypred)     
       mdr[idv,2] <- rlb$aic     
     } # fin pour idv     
          
     # tri par AUROC décroissante     
          
     idx <- order(mdr[,1],decreasing=TRUE)     
     mdr <- mdr[idx,]     
     cats("Meilleures régressions logistiques simples pour CHD69 :")     
     print( round(mdr,2) )     
          
     # analyses complémentaires     
          
     wcgs <- cbind(WCGS$CHD69,wcgs)     
     names(wcgs)[1] <- "CHD69"     
          
     mdc(wcgs)     
          
     gr("wcgs_pairs1.png")     
     pairsi(wcgs)     
     dev.off()     
          
          

          
     Meilleures régressions logistiques simples pour CHD69 :     
     =======================================================     
          
           AUROC    AIC     
     CHOL   0.66 703.84     
     SBP    0.64 724.07     
     AGE    0.62 735.22     
     SMOKE  0.62 757.56     
     BMI    0.57 768.15     
          
     Matrice des corrélations au sens de  Pearson pour  3154  lignes et  6  colonnes     
          
           CHD69   AGE  CHOL   SBP    BMI SMOKE     
     CHD69 1.000     
     AGE   0.119 1.000     
     CHOL  0.163 0.089 1.000     
     SBP   0.133 0.166 0.123 1.000     
     BMI   0.062 0.027 0.071 0.288  1.000     
     SMOKE 0.085 0.005 0.097 0.003 -0.142 1.000     
          

          wcgs_pairs 1


     # lecture des données     
          
     WCGS       <- lit.dar("wcgs.dar")     
     lesv1      <- c("AGE","CHOL","SBP","BMI","SMOKE")     
     lesv2      <- c("ARCUS","DIBPAT")     
     wcgs       <- WCGS[,c(lesv1,lesv2)]     
          
     # analyses préliminaires     
          
     wcgs <- cbind(WCGS$CHD69,wcgs)     
     names(wcgs)[1] <- "CHD69"     
          
     mdc(wcgs)     
          
     gr("wcgs_pairs2.png")     
     pairsi(wcgs)     
     dev.off()     
          
     # meilleur modèle     
          
     rchModeleLogistique(df=wcgs)     
          

          
     Matrice des corrélations au sens de  Pearson pour  3154  lignes et  8  colonnes     
          
            CHD69   AGE  CHOL   SBP    BMI SMOKE ARCUS DIBPAT     
     CHD69  1.000     
     AGE    0.119 1.000     
     CHOL   0.163 0.089 1.000     
     SBP    0.133 0.166 0.123 1.000     
     BMI    0.062 0.027 0.071 0.288  1.000     
     SMOKE  0.085 0.005 0.097 0.003 -0.142 1.000     
     ARCUS  0.066 0.187 0.121 0.034 -0.035 0.073 1.000     
     DIBPAT 0.112 0.088 0.057 0.077  0.028 0.061 0.044  1.000     
     RStudioGD     
             2     
                 AUROC (Intercept)        AGE       CHOL        SBP        BMI     SMOKE     ARCUS    DIBPAT      AIC     
     M07V00  0.7503843  -12.274385 0.05718452 0.01038885 0.01821424 0.05948104 0.5986052 0.2190833 0.7021576 1594.445     
     M07V01  0.6199228   -5.939516 0.07442256         NA         NA         NA        NA        NA        NA 1742.356     
     M07V02  0.6624719   -5.359022         NA 0.01244778         NA         NA        NA        NA        NA 1706.423     
     M07V03  0.6313891   -5.926461         NA         NA 0.02667129         NA        NA        NA        NA 1736.441     
     M07V04  0.5661082   -4.508689         NA         NA         NA 0.08427681        NA        NA        NA 1773.426     
     M07V05  0.5775470   -2.763620         NA         NA         NA         NA 0.6298631        NA        NA 1762.381     
     M07V06  0.5551950   -2.599052         NA         NA         NA         NA        NA 0.4918141        NA 1762.216     
     M07V07  0.6027757   -2.934395         NA         NA         NA         NA        NA        NA 0.8641250 1744.344     
     M07V-01 0.7368158   -9.906505         NA 0.01033884 0.02161026 0.05373681 0.5563023 0.3531192 0.7427474 1614.084     
     M07V-02 0.7175185  -10.335688 0.05680841         NA 0.02009451 0.06780513 0.6615986 0.3156082 0.7252050 1643.590     
     M07V-03 0.7411357  -11.165558 0.06534138 0.01082220         NA 0.09108547 0.6135754 0.1868962 0.7335877 1610.684     
     M07V-04 0.7455221  -11.070892 0.05591968 0.01057667 0.02058945         NA 0.5478250 0.2089027 0.7086461 1597.388     
     M07V-05 0.7396840  -11.526847 0.05360428 0.01093851 0.01857986 0.04111321        NA 0.2521781 0.7282620 1610.730     
     M07V-06 0.7490711  -12.291680 0.06037041 0.01072095 0.01807657 0.05512504 0.6024262        NA 0.6971798 1604.035     
     M07V-07 0.7345037  -12.313392 0.06118345 0.01052450 0.01937777 0.06209862 0.6284042 0.2202162        NA 1617.149     
     M07VBS  0.7503843  -12.274385 0.05718452 0.01038885 0.01821424 0.05948104 0.5986052 0.2190833 0.7021576 1594.445     
     M07VFS  0.7503843  -12.274385 0.05718452 0.01038885 0.01821424 0.05948104 0.5986052 0.2190833 0.7021576 1594.445     
     M07VMS  0.7503843  -12.274385 0.05718452 0.01038885 0.01821424 0.05948104 0.5986052 0.2190833 0.7021576 1594.445     
          

          wcgs_pairs 2

La fonction step() permet aussi de déterminer le meilleur modèle par sélection ascendante et descendante de variables pour les modèles linéiares généralisésvariables :


     # données sans NA     
          
     wcgsGLM <- na.omit(wcgs[,(1:6)])     
          
     # modèle minimal     
          
     glmMin <- glm(as.factor(CHD69)~1,data=wcgsGLM,family=binomial() )     
          
     # modèle maxnimal     
          
     glmMax <- glm(as.factor(CHD69)~.,data=wcgsGLM,family=binomial() )     
          
     # sélection forward du meilleur modèle     
          
     cats("sélection ascendante")     
     step(glmMin,scope=list(upper=glmMax),trace=1,direction="forward")     
          
     # sélection backward du meilleur modèle     
          
     cats("sélection descendante")     
     step(glmMax,scope=list(lower=glmMin),trace=1,direction="backward")     
          
          

          
     sélection ascendante     
     ====================     
          
     Start:  AIC=1781.2     
     as.factor(CHD69) ~ 1     
          
             Df Deviance    AIC     
     + CHOL   1   1702.4 1706.4     
     + SBP    1   1729.5 1733.5     
     + AGE    1   1736.0 1740.0     
     + SMOKE  1   1756.2 1760.2     
     + BMI    1   1767.2 1771.2     
     <none>       1779.2 1781.2     
          
     Step:  AIC=1706.42     
     as.factor(CHD69) ~ CHOL     
          
             Df Deviance    AIC     
     + SBP    1   1664.1 1670.1     
     + AGE    1   1665.8 1671.8     
     + SMOKE  1   1686.6 1692.6     
     + BMI    1   1693.7 1699.7     
     <none>       1702.4 1706.4     
          
     Step:  AIC=1670.08     
     as.factor(CHD69) ~ CHOL + SBP     
          
             Df Deviance    AIC     
     + AGE    1   1637.8 1645.8     
     + SMOKE  1   1647.2 1655.2     
     <none>       1664.1 1670.1     
     + BMI    1   1662.5 1670.5     
          
     Step:  AIC=1645.8     
     as.factor(CHD69) ~ CHOL + SBP + AGE     
          
             Df Deviance    AIC     
     + SMOKE  1   1619.3 1629.3     
     + BMI    1   1635.6 1645.6     
     <none>       1637.8 1645.8     
          
     Step:  AIC=1629.34     
     as.factor(CHD69) ~ CHOL + SBP + AGE + SMOKE     
          
            Df Deviance    AIC     
     + BMI   1   1614.6 1626.6     
     <none>      1619.3 1629.3     
          
     Step:  AIC=1626.61     
     as.factor(CHD69) ~ CHOL + SBP + AGE + SMOKE + BMI     
          
          
     Call:  glm(formula = as.factor(CHD69) ~ CHOL + SBP + AGE + SMOKE + BMI,     
         family = binomial(), data = wcgsGLM)     
          
     Coefficients:     
     (Intercept)         CHOL          SBP          AGE        SMOKE          BMI     
       -12.33630      0.01084      0.01930      0.06436      0.63269      0.05767     
          
     Degrees of Freedom: 3141 Total (i.e. Null);  3136 Residual     
     Null Deviance:	    1779     
     Residual Deviance: 1615 	AIC: 1627     
          
     sélection descendante     
     =====================     
          
     Start:  AIC=1626.61     
     as.factor(CHD69) ~ AGE + CHOL + SBP + BMI + SMOKE     
          
             Df Deviance    AIC     
     <none>       1614.6 1626.6     
     - BMI    1   1619.3 1629.3     
     - SMOKE  1   1635.6 1645.6     
     - SBP    1   1635.7 1645.7     
     - AGE    1   1643.6 1653.6     
     - CHOL   1   1668.8 1678.8     
          
     Call:  glm(formula = as.factor(CHD69) ~ AGE + CHOL + SBP + BMI + SMOKE,     
         family = binomial(), data = wcgsGLM)     
          
     Coefficients:     
     (Intercept)          AGE         CHOL          SBP          BMI        SMOKE     
       -12.33630      0.06436      0.01084      0.01930      0.05767      0.63269     
          
     Degrees of Freedom: 3141 Total (i.e. Null);  3136 Residual     
     Null Deviance:	    1779     
     Residual Deviance: 1615 	AIC: 1627     
          

 

 

Code-source php de cette page. Retour à la page principale du cours.

 

 

retour gH    Retour à la page principale de   (gH)