Valid XHTML     Valid CSS2    

Module de Biostatistiques,

partie 2

Ecole Doctorale Biologie Santé

gilles.hunault "at" univ-angers.fr

      

Les dates de formation sont ici.

 

Solutions pour la séance numéro 4 (énoncés)

  1. La cote est liée à la probabilité par la la relation cote = proba / (1-proba) soit encore proba = cote / (1+cote). Ainsi une cote de 3 pour 1 correspond à une probabilité de 0,75 et une proba de 0,5 à une cote de 1 pour 1. Voici les fonctions associées :

    
         proba2cote <- function(proba) {     
              
          return( proba/(1-proba) )     
              
         } # fin de fonction proba2cote     
              
         ####################################     
              
         cote2proba <- function(cote) {     
              
          return( cote/(1+cote) )     
              
         } # fin de fonction cote2proba     
              
         ####################################     
         ####################################     
              
         cat(" à une cote de de 3/1 correspond une proba de ",cote2proba(3/1),"\n")     
              
         cat(" à une probabilité de 0,5 correspond une cote de ",proba2cote(0.5),"\n")     
              
              
    
    
          à une cote de de 3/1 correspond une proba de  0.75     
          à une probabilité de 0,5 correspond une cote de  1     
              
    

    Le rapport de chances ou rapport de cotes ou OR (odds ratio) est une valeur numérique qui quantifie la dépendance et l'importance d'évènements binaires.

    Voir à ce sujet et dans cet ordre le wiki français (court) puis le wiki anglais (long).

    Il ne faut pas confondre le rapport de cotes et le risque relatif RR, en anglais : relative risk. Pour faire court : RR est le rapport des risques et OR est le rapport des cotes.

  2. Si dans la régression linéaire classique on modélise la variable continue quantitative Y par une combinaison linéaire de variables xi (notée βX), en régression logistique, on modélise une variable Y qualitative ce qui inclut le cas d'une variable Y binaire ou («logique»). Le modèle linéaire de base s'écrit Y = βX alors que le modèle linéaire généralisé s'écrit g(Y) = βXg est nommée fonction de lien. Pour la régression logistique, la fonction de lien est la fonction logit définie par logit(p)=log(p/(1-p)) dont l'inverse est la fonction logistique qui est donc définie par logis(x)=exp(x)/(1+exp(x)).

    On lira avec profit le document MODGEN du CNAM pour comprendre ce que sont les modèles linéaires généralisés, puis le document semin-R pour voir un exemple d'application.

    Un long document en français (270 pages) sur le sujet la régression logistique est celui de R. Rakotomalala disponible à l'adresse

    http://eric.univ-lyon2.fr/~ricco/cours/cours/pratique_regression_logistique.pdf

    dont une copie locale est ici.

    Pour réaliser une régression logistique en R, on utilise la fonction glm() du package stats au lieu de la fonction lm() du package stats pour la régression linéaire et on précise que la famille de fonctions à utiliser est binomial.

    L'affichage du modèle via la fonction générique print() du package base n'est en général pas suffisant et on lui préfère la fonction summary.glm() du package stats mais comme pour la fonction summary.lm() de ce même package stats on peut se contenter d'utiliser la fonction générique summary() du package base. En d'autres termes, si rlb est un modèle issu d'une régression logistique, on peut utiliser summary(rlb) mais si on veut de l'aide sur la "vraie" fonction summary() utilisée pour ce modèle, il faut écrire help(summary.glm). Il est en général conseillé d'appliquer aussi la fonction générique anova() sur un modèle issu d'une régression logistique, ce qui a pour but d'exécuter la fonction anova.glm().

    Comme pour la fonction lm() et sa fonction predict.lm(), la fonction glm() dispose d'une fonction predict.glm() disponible via la fonction générique predict().

    De la même façon, comme pour la fonction lm() et sa fonction plot.lm(), la fonction glm() dispose d'une fonction plot.glm() disponible via la fonction générique plot() qui permet de réaliser une analyse graphique des résidus.

  3. Il est bien sûr très fortement conseillé d'étudier les deux variables TAILLE et GROUPE séparément puis conjointement avant d'effectuer la régression logistique binaire de GROUPE en fonction de TAILLE, soit le code R :

    
         ## avant la régression logistique binaire simple de GROUPE en fonction de TAILLE     
         ## dans les données du fichier pg.dar, on analyse les variables     
              
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données du fichier-texte     
              
         pg <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/pg.dar")     
              
         # description de TAILLE     
              
         gr("pgTAILLE.png")     
         decritQT("TAILLE",pg$TAILLE,"cm",TRUE)     
         dev.off()     
              
         # description de GROUPE     
              
         gr("pgGROUPE.png")     
         decritQL("GROUPE",pg$GROUPE,"petits grands",TRUE)     
         dev.off()     
              
         # description de TAILLE par GROUPE     
              
         gr("pgTAILLEparGROUPE.png")     
         decritQTparFacteur("TAILLE",pg$TAILLE,"cm","GROUPE",pg$GROUPE,"petits grands",TRUE)     
         dev.off()     
              
    
    
         DESCRIPTION STATISTIQUE DE LA VARIABLE  TAILLE     
              
          Taille                30        individus     
          Moyenne              159.9000          cm     
          Ecart-type            18.8869          cm     
          Coef. de variation    12               %     
          1er Quartile         140.0000          cm     
          Mediane              165.0000          cm     
          3eme Quartile        175.0000          cm     
          iqr absolu            35.0000          cm     
          iqr relatif           21.0000          %     
          Minimum               130.000          cm     
          Maximum               190.000          cm     
              
          Tracé tige et feuilles     
              
           The decimal point is 1 digit(s) to the right of the |     
              
           13 | 0025566     
           14 | 001     
           15 |     
           16 | 012455889     
           17 | 0355789     
           18 | 012     
           19 | 0     
              
              
         TRI A PLAT DE :  GROUPE (ordre des modalités)     
              
                           petits grands  Total     
          Effectif             10     20     30     
          Cumul Effectif       10     30     30     
          Frequence (en %)     33     67    100     
          Cumul fréquences     33    100    100     
              
         VARIABLE :  GROUPE (par fréquence décroissante)     
              
                            grands petits  Total     
           Effectif             20     10     30     
           Cumul Effectif       20     30     30     
           Frequence (en %)     67     33    100     
           Cumul fréquences     67    100    100     
              
         [1] "petits" "grands"     
              
              
              
         VARIABLE QT  TAILLE ,unit=cm     
         VARIABLE QL  GROUPE  labels :  petits grands     
              
                     N       Moy Unite       Ect Cdv        Q1       Med        Q3       EIQ       Min       Max     
         Global     30   159.900    cm    18.569  12   140.000   165.000   175.000    35.000   130.000   190.000     
         petits     10   137.600    cm     8.789   6   132.750   135.500   139.000     6.250   130.000   162.000     
         grands     20   171.050    cm    10.278   6   165.000   171.500   178.250    13.250   141.000   190.000     
              
         Analysis of Variance Table     
              
         Response: nomVarQT     
                   Df Sum Sq Mean Sq F value    Pr(>F)     
         ql         1 7459.4  7459.4  72.387 2.998e-09 ***     
         Residuals 28 2885.4   103.0     
         ---     
         Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
         ANOVA : there is a significant diffence at the level 0.05 for  TAILLE   GROUPE     
              
         	Kruskal-Wallis rank sum test     
              
         data:  nomVarQT by ql     
         Kruskal-Wallis chi-squared = 18.239, df = 1, p-value = 1.948e-05     
              
         KRUSKAL-WALLIS : there is a significant diffence at the level 0.05 for  TAILLE   GROUPE     
              
              
    

    Une lecture détaillée de ces résultats et des graphiques montre qu'il y a bien sans doute deux groupes de taille.

              non su

    Le code R suivant réalise la régression logistique de GROUPE en fonction de TAILLE

    
         ## exemple de régression logistique binaire simple :     
         ## on modélise GROUPE en fonction de TAILLE dans les données     
         ## du fichier pg.dar     
              
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données du fichier-texte     
              
         pg <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/pg.dar")     
              
         # calcul de la régression logistique     
              
         rlb <- glm(GROUPE ~ TAILLE,data=pg,family="binomial")     
              
         # affichage simplifié (coefficients)     
              
         cats("Régression logistique de GROUPE en fonction de TAILLE")     
              
         print(rlb)     
              
         # affichage détaillé     
              
         cats("Détails de la régression logistique")     
         print(summary(rlb))     
              
    
    
         Régression logistique de GROUPE en fonction de TAILLE     
         =====================================================     
              
              
         Call:  glm(formula = GROUPE ~ TAILLE, family = "binomial", data = pg)     
              
         Coefficients:     
         (Intercept)       TAILLE     
            -27.2103       0.1812     
              
         Degrees of Freedom: 29 Total (i.e. Null);  28 Residual     
         Null Deviance:	    38.19     
         Residual Deviance: 10.89 	AIC: 14.89     
              
         Détails de la régression logistique     
         ===================================     
              
              
         Call:     
         glm(formula = GROUPE ~ TAILLE, family = "binomial", data = pg)     
              
         Deviance Residuals:     
             Min       1Q   Median       3Q      Max     
         -2.1225  -0.2589   0.1087   0.2728   1.9168     
              
         Coefficients:     
                      Estimate Std. Error z value Pr(>|z|)     
         (Intercept) -27.21029    8.95762  -3.038  0.00238 **     
         TAILLE        0.18118    0.05893   3.074  0.00211 **     
         ---     
         Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
              
         (Dispersion parameter for binomial family taken to be 1)     
              
             Null deviance: 38.191  on 29  degrees of freedom     
         Residual deviance: 10.895  on 28  degrees of freedom     
         AIC: 14.895     
              
         Number of Fisher Scoring iterations: 6     
              
              
    

    Comme pour une régression linéaire, il est possible de connaitre les valeurs calculées par le modèle et d'en prédire de nouvelles :

    
         ## régression logistique binaire simple et prédictions     
         ## on prédit des données pour les tailles 120 et 150 cm     
              
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données du fichier-texte     
              
         pg <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/pg.dar")     
              
         # calcul de la régression logistique     
              
         rlb <- glm(GROUPE ~ TAILLE,data=pg,family="binomial")     
              
         # affichage des valeurs calculées par le modèle     
         # (données triées par groupe puis taille croissante)     
              
         npg <- cbind(pg,rlb$fitted.values)     
         ord <- order(pg$GROUPE,pg$TAILLE)     
         print(npg[ord,])     
              
         # prédictions     
              
         cats("Valeurs prédites")     
         nTailles <- data.frame( TAILLE=c(120,150) )     
         nGroupes <- predict(rlb,nTailles,"response")     
         print( cbind(nTailles,nGroupes) )     
              
         # tracé des valeurs initiales et prédites     
              
         coulGrp <- ifelse(pg$GROUPE==0,"blue","red")     
         plot(GROUPE ~ TAILLE,data=pg,xlim=c(110,200),col=coulGrp,type="p",pch=19)     
         title(main="Exemple de régression logistique")     
         points( as.numeric(nGroupes) ~ unlist(nTailles), col="green", cex=1.5,pch=19 )     
              
         # tracé de la fonction logistique du modèle     
              
         x <- data.frame( TAILLE=seq(from=min(pg$TAILLE),to=max(pg$TAILLE),length.out=1000) )     
         y <- predict(rlb,x,"response")     
         points( as.numeric(y) ~ unlist(x), col="black", type="l")     
              
         # valeur-seuil par défaut : 0.5     
              
         abline(h=0.5)     
              
    
    
             TAILLE GROUPE rlb$fitted.values     
         A01    130      0        0.02517286     
         A04    130      0        0.02517286     
         A09    132      0        0.03577320     
         A05    135      0        0.06005399     
         A07    135      0        0.06005399     
         A03    136      0        0.07113425     
         A08    136      0        0.07113425     
         A02    140      0        0.13650072     
         A06    140      0        0.13650072     
         C01    162      0        0.89485898     
         A10    141      1        0.15929538     
         C02    160      1        0.85557305     
         C04    161      1        0.87655250     
         C20    164      1        0.92440288     
         C03    165      1        0.93613049     
         C07    165      1        0.93613049     
         C08    168      1        0.96189141     
         C12    168      1        0.96189141     
         C16    169      1        0.96800461     
         C11    170      1        0.97316451     
         C15    173      1        0.98423973     
         C13    175      1        0.98897761     
         C19    175      1        0.98897761     
         C10    177      1        0.99230232     
         C17    178      1        0.99356976     
         C18    179      1        0.99462964     
         C05    180      1        0.99551561     
         C14    181      1        0.99625597     
         C09    182      1        0.99687448     
         C06    190      1        0.99926469     
              
         Valeurs prédites     
         ================     
              
           TAILLE    nGroupes     
         1    120 0.004200577     
         2    150 0.491792541     
              
    

              non su

  4. Un simple appel de la fonction ifelse() du package base permet de déterminer les groupes prédits une fois le seuil de séparation choisi. La matrice de confusion est alors simplement le tri croisé des variables groupesOriginaux et groupesPredits. Le taux de bien classés s'obtient comme le pourcentage associé aux valeurs de la diagonale :

    
         ## qualité de la régression sur les données pg.dar avec un seuil de séparation à 0,5
         
         # calcul de la régression
         
         rlb   <- glm(GROUPE ~ TAILLE,data=pg,family="binomial")
         
         # valeurs prédites
         
         ypred <- predict(rlb,newdata=data.frame(TAILLE=pg$TAILLE),type="response")
         
         # groupes prédits
         
         tdr   <- data.frame(
                   taille=pg$TAILLE,
                   groupeOriginal=pg$GROUPE,
                   groupePredit= ifelse( ypred<0.5,0,1)
         ) # fin de data.frame
         
         print(tdr)
         
         cats("Matrice de confusion")
         
         with(tdr,print(table(
          groupeOriginal,
          groupePredit
         
         )))
         
         cats("Nombre et taux de bien classés")
         
         nbTot <- nrow(tdr)
         nbBc  <- sum(tdr$groupeOriginal==tdr$groupePredit)
         pctBc <- 100*nbBc/nbTot
         
         cat(nbBc,"bien classés sur",nbTot,"individus en tout, soit",pctBc,"%\n")
         
    
    
            taille groupeOriginal groupePredit
         1     130              0            0
         2     140              0            0
         3     162              0            1
         4     160              1            1
         5     136              0            0
         6     165              1            1
         7     130              0            0
         8     135              0            0
         9     140              0            0
         10    135              0            0
         11    161              1            1
         12    136              0            0
         13    180              1            1
         14    190              1            1
         15    132              0            0
         16    141              1            0
         17    165              1            1
         18    168              1            1
         19    182              1            1
         20    177              1            1
         21    170              1            1
         22    168              1            1
         23    175              1            1
         24    181              1            1
         25    173              1            1
         26    169              1            1
         27    178              1            1
         28    179              1            1
         29    175              1            1
         30    164              1            1
         
         Matrice de confusion
         ====================
         
                       groupePredit
         groupeOriginal  0  1
                      0  9  1
                      1  1 19
         
         Nombre et taux de bien classés
         ==============================
         
         28 bien classés sur 30 individus en tout, soit 93.33333 %
         
    

    Avec un seuil de séparation à 0,15, tous les individus du groupe 1 sont bien prédits et le taux de bien classés est alors meilleur, mais bien sûr ce seuil est arbitraire et lié aux données...

    
         ## qualité de la régression sur les données pg.dar avec un seuil de séparation à 0,15
         
         # calcul de la régression
         
         rlb2   <- glm(GROUPE ~ TAILLE,data=pg,family="binomial")
         
         # valeurs prédites
         
         ypred <- predict(rlb2,newdata=data.frame(TAILLE=pg$TAILLE),type="response")
         
         # groupes prédits
         
         tdr2   <- data.frame(
                   taille=pg$TAILLE,
                   groupeOriginal=pg$GROUPE,
                   groupePredit= ifelse( ypred<0.15,0,1)
         ) # fin de data.frame
         
         cats("Matrice de confusion")
         
         with(tdr2,print(table(
          groupeOriginal,
          groupePredit
         
         )))
         
         cats("Nombre et taux de bien classés")
         
         nbTot <- nrow(tdr2)
         nbBc  <- sum(tdr2$groupeOriginal==tdr2$groupePredit)
         pctBc <- 100*nbBc/nbTot
         
         cat(nbBc,"bien classés sur",nbTot,"individus en tout, soit",pctBc,"%\n")
         
    
    
         
         Matrice de confusion
         ====================
         
                       groupePredit
         groupeOriginal  0  1
                      0  9  1
                      1  0 20
         
         Nombre et taux de bien classés
         ==============================
         
         29 bien classés sur 30 individus en tout, soit 96.66667 %
         
    
  5. Commençons par regarder les variables CHD69 et AGE séparément puis conjointement.

    
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données     
              
         WCGS       <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/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   Max     
         Global   3154    46.279   ans     5.523  12    42.000    45.000    50.000     8.000    39.000   59.000     
         Non      2897    46.082   ans     5.456  12    41.000    45.000    50.000     9.000    39.000   59.000     
         Oui       257    48.490   ans     5.790  12    44.000    49.000    53.000     9.000    39.000   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 concure 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 de CHD 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 [pseudo]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 β1fois 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

    On pourra comparer ces résultats à ceux fournis par Sas :

    
              
         OPTIONS nodate LINESIZE=256 PAGESIZE=50 NOCENTER FORMDLIM=' ' FORMCHAR='-----------' nonumber ;     
              
         PROC IMPORT datafile="wcgs.xls" dbms=excel out=wcgs replace ;     
              
         * Etude de l'age ;     
              
         PROC MEANS data=wcgs N MEAN STD CV MIN MAX MAXDEC=2 ;     
              var age ;     
              
         * Etude de la cardiopathie coronarienne ;     
              
         PROC FREQ data=wcgs ;     
              table chd69 ;     
              
         * Régression logistique de la cardiopathie coronarienne;     
         * en fonction de l'age ;     
              
         PROC LOGISTIC data=wcgs simple descending  ;     
              model chd69 = age / expb rsquare clparm=both ;     
              units age = 1 5 10 20 ;     
              
         RUN ;     
              
              
    
    
              
         NOTE: SAS (r) Proprietary Software 9.2 (TS2M3)     
         NOTE: Copyright (c) 2002-2008 by SAS Institute Inc., Cary, NC, USA.     
               Licensed to LICENCE GRATUITE ETUDIANTS PROFS MENR, Site 50101271.     
         NOTE: La session est exécutée sur la plate-forme XP_PRO .     
              
              
              
         NOTE: L'initialisation de SAS a utilisé :     
               temps réel          0.65 secondes     
               temps UC            0.60 secondes     
              
         1          /*  TimeStamp  (dos) : 30 Octobre 10 17:33  ;  z:\wcgs.sas */     
         2     
         3          OPTIONS nodate LINESIZE=256 PAGESIZE=50 NOCENTER FORMDLIM=' ' FORMCHAR='-----------' nonumber ;     
         4     
         5          PROC IMPORT datafile="wcgs.xls" dbms=excel out=wcgs replace ;     
         6     
         7          * Etude de l'age ;     
         8     
              
         NOTE: WORK.WCGS data set was successfully created.     
         NOTE: Procédure IMPORT a utilisé (Durée totale du traitement) :     
               temps réel          0.95 secondes     
               temps UC            0.67 secondes     
              
              
         9          PROC MEANS data=wcgs N MEAN STD CV MIN MAX MAXDEC=2 ;     
         10              var age ;     
         11     
         12         * Etude de la cardiopathie coronarienne ;     
         13     
              
         NOTE:  3154 observation(s) lue(s) dans la table WORK.WCGS.     
         NOTE: La procédure MEANS a été imprimée en page 1.     
         NOTE: Procédure MEANS a utilisé (Durée totale du traitement) :     
               temps réel          0.25 secondes     
               temps UC            0.13 secondes     
              
              
         14         PROC FREQ data=wcgs ;     
         15              table chd69 ;     
         16     
         17         * Régression logistique de la cardiopathie coronarienne;     
         18         * en fonction de l'age ;     
         19     
              
         NOTE:  3154 observation(s) lue(s) dans la table WORK.WCGS.     
         NOTE: La procédure FREQ a été imprimée en page 2.     
         Le Système SAS     
              
         NOTE: Procédure FREQ a utilisé (Durée totale du traitement) :     
               temps réel          0.03 secondes     
               temps UC            0.03 secondes     
              
              
         20         PROC LOGISTIC data=wcgs simple descending  ;     
         21              model chd69 = age / expb rsquare clparm=both ;     
         22     
              
         NOTE: PROC LOGISTIC is modeling the probability that chd69=1.     
         NOTE: Critère de convergence (GCONV=1E-8) respecté.     
         NOTE:  3154 observation(s) lue(s) dans la table WORK.WCGS.     
         NOTE: La procédure LOGISTIC a été imprimée en pages 3-5.     
         NOTE: Procédure LOGISTIC a utilisé (Durée totale du traitement) :     
               temps réel          0.59 secondes     
               temps UC            0.57 secondes     
              
              
         23         PROC LOGISTIC data=wcgs descending  ;     
         24              model chd69 = age ;     
         25              units age = 1 5 10 20 ;     
         26     
         27         RUN ;     
              
         NOTE: PROC LOGISTIC is modeling the probability that chd69=1.     
         NOTE: Critère de convergence (GCONV=1E-8) respecté.     
         NOTE:  3154 observation(s) lue(s) dans la table WORK.WCGS.     
         NOTE: La procédure LOGISTIC a été imprimée en pages 6-7.     
         NOTE: Procédure LOGISTIC a utilisé (Durée totale du traitement) :     
               temps réel          0.51 secondes     
               temps UC            0.50 secondes     
              
              
         NOTE: SAS Institute Inc., SAS Campus Drive, Cary, NC USA 27513-2414     
         NOTE: Le Système SAS a utilisé :     
               temps réel          3.12 secondes     
               temps UC            2.56 secondes     
              
              
    
    
              
         Le Système SAS     
              
         Procédure MEANS     
              
         Variable d'analyse : age     
                                                      Coef de     
            N         Moyenne      Ecart-type       variation         Minimum         Maximum     
         ------------------------------------------------------------------------------------     
         3154           46.28            5.52           11.94           39.00           59.00     
         ------------------------------------------------------------------------------------     
              
         Procédure FREQ     
              
                                                      Fréquence    Pctage.     
         chd69    Fréquence           Pourcentage     cumulée      cumulé     
         -----------------------------------------------------------------     
             0        2897            91.85               2897      91.85     
             1         257             8.15               3154     100.00     
              
         Procédure LOGISTIC     
              
         Informations sur le modèle     
              
         Table                            WORK.WCGS     
         Variable de réponse              chd69     
         Nombre de niveaux de réponse     2     
         Modèle                           logit binaire     
         Technique d'optimisation         Score de Fisher     
              
         Nombre d'observations lues         3154     
         Nombre d'observations utili        3154     
              
         Profil de réponse     
              
           Valeur                  Fréquence     
         ordonnée        chd69        totale     
              
                1            1           257     
                2            0          2897     
              
         La probabilité modélisée est chd69=1.     
              
         Statistiques descriptives pour variables continues     
                                                          Ecart-     
         Variable     chd69             Moyenne             type          Minimum          Maximum     
              
         age                 1        48.490272         5.801477        39.000000        59.000000     
                             0        46.082499         5.456675        39.000000        59.000000     
         ------------------------------------------------------------------------------------------     
                         Total        46.278694         5.524045        39.000000        59.000000     
              
              
         Etat de convergence du modèle     
         Critère de convergence (GCONV=1E-8) respecté.     
              
         Statistiques d'ajustement du modèle     
              
                                       Constante     
                       Constante              et     
         Critère      uniquement     covariables     
              
         AIC            1783.244        1742.356     
         SC             1789.300        1754.469     
         -2 Log L       1781.244        1738.356     
              
              
         R carré    0.0135    R carré remis à l'échelle max.    0.0313     
              
              
         Test de l'hypothèse nulle globale : BETA=0     
              
         Test                      Khi-2      DDL     Pr > Khi-2     
         Rapp. de vrais.         42.8876        1         <.0001     
         Score                   44.8616        1         <.0001     
         Wald                    43.3584        1         <.0001     
              
         Estimations par l'analyse du maximum de vraisemblance     
              
                               Valeur      Erreur         Khi-2     
         Paramètre    DDL     estimée        type       de Wald    Pr > Khi-2    Exp(Est)     
              
         Intercept      1     -5.9395      0.5493      116.9098        <.0001       0.003     
         age            1      0.0744      0.0113       43.3584        <.0001       1.077     
              
              
         Estimations des rapports de cotes     
              
                    Valeur     
                   estimée    Intervalle de confiance     
         Effet    du point         de Wald à 95 %     
              
         age         1.077       1.054          1.101     
              
         Association des probabilités prédites et des réponses observées     
              
         Pourcentage concordant      59.5    D de Somers    0.240     
         Pourcentage discordant      35.6    Gamma          0.252     
         Pourcentage lié              4.9    Tau-a          0.036     
         Paires                    744529    c              0.620     
              
              
         Intervalle de confiance de vraisemblance du profil pour les paramètres     
              
                         Valeur     
         Paramètre      estimée     Intervalle de confiance à 95 %     
              
         Intercept      -5.9395      -7.0245               -4.8696     
         age             0.0744       0.0523                0.0966     
              
         Intervalle de confiance de Wald pour les paramètres     
              
                         Valeur     
         Paramètre      estimée     Intervalle de confiance à 95 %     
              
         Intercept      -5.9395      -7.0162               -4.8629     
         age             0.0744       0.0523                0.0966     
              
         Rapports de cotes     
                                  Valeur     
         Effet        Unité      estimée     
              
         age         1.0000        1.077     
         age         5.0000        1.451     
         age        10.0000        2.105     
         age        20.0000        4.430     
              
              
    

    Pour savoir si c'est un bon modèle, on peut essayer un seuil de séparation à 0,5 et calculer le taux de bien-classés, comme précédemment  :

    
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données     
              
         WCGS       <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/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 )     
              
         # prédictions et taux total de bien classés     
              
         ypred <- predict(rlb,newdata=data.frame(CHD69=WCGS$CHD69),type="response")     
         tdr   <- data.frame(chdOrg=WCGS$CHD69,chdPred=ifelse(ypred<0.5,0,1))     
              
         cats("Matrice de confusion")     
              
         with(tdr,print(table(     
               chdOrg,     
               chdPred     
              
         )))     
              
         cats("Nombre et taux de bien classés")     
              
         nbTot <- nrow(tdr)     
         nbBc  <- sum(tdr$chdOrg==tdr$chdPred)     
         pctBc <- 100*nbBc/nbTot     
              
         cat(nbBc,"bien classés sur",nbTot,"individus en tout, soit",pctBc,"%\n")     
              
    
    
         Call:     
         glm(formula = chd ~ age, family = binomial())     
              
         Deviance Residuals:     
             Min       1Q   Median       3Q      Max     
         -0.6209  -0.4545  -0.3669  -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     
              
         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     
              
         Matrice de confusion     
         ====================     
              
               chdPred     
         chdOrg    0     
              0 2897     
              1  257     
              
         Nombre et taux de bien classés     
         ==============================     
              
         2897 bien classés sur 3154 individus en tout, soit 91.85162 %     
              
    
  6. Comme pour une régression linéaire, il n'est pas possible d'utiliser la fonction plot() pour afficher le modèle parce que cette fonction plot() affiche l'analyse graphique des résidus. Il faut donc construire des points où x (l'age) couvre une grande plage de variation et où y est la valeur prédite par ce modèle, soit le code R suivant :

    
         ## tracé d'une régression logistique binaire
         
         # chargement des fonctions gH :
         
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")
         
         # lecture des données
         
         WCGS       <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/wcgs.dar")
         age        <- WCGS$AGE
         chd        <- factor(WCGS$CHD69,labels=c("Non","Oui"))
         
         # calcul du modèle régression logistique
         
         rlb <- glm( formula = chd ~ age, family = binomial() )
         
         # tracé de la variable binaire en fonction de la variable explicative
         
         chd01 <- WCGS$CHD69
         couleurCHD <- ifelse(chd01==0,"green","red")
         plot( chd01 ~ age,type="p",main="CHD en fonction de AGE",pch=19,col=couleurCHD,xlim=c(0,200))
         
         # tracé de la fonction logistique du modèle
         
         x <- data.frame( age=seq(from=0,to=200,length.out=1000))
         y <- predict(rlb,x,"response")
         points( as.numeric(y) ~ unlist(x), col="black", type="l")
         
         # valeur-seuil par défaut : 0.5
         
         abline(h=0.5,col="blue")
         
         
    

    On se rend sans doute compte ici que la régression logistique n'est pas "bonne".

              non su

  7. La sensibilité et la spécificité sont des indicateurs diagnostiques qui permettent de décrire la performance d'une prédiction. Ce sont deux indicateurs parmi de nombreux autres, obtenus à partir des 4 valeurs entières VP, FP, FN et VP. Voir notre page étude de valeurs diagnostiques pour plus d'explications avec des exemples numériques.

    ROC signifie Receiver Operating Characteristic (mais la page wiki en est plus complète, comme souvent) alors que AUROC signifie Area Under the ROC Curve. Une courbe ROC modélise les performances de la prédiction en fonction des différents seuils possibles et l'AUROC est l'aire sous la courbe. Cette aire est utilisée comme indicateur global de la performance, moins sensible à la prévalence que la sensibilité et la spécificité.

    On lira avec intérêt les transparents de E. Rakotomalala et, pour une approche plus bioinformatique : évaluation proteome avant de regarder notre page de calcul d'AUROC en direct pour se familiariser avec ces notions.

    Pour calculer une AUROC avec R, de nombreux packages existent, par exemple le package limma avec sa fonction auROC() et le package AUC et sa fonction auc() :

    
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données     
              
         WCGS       <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/wcgs.dar")     
         age    <- WCGS$AGE     
         chdbin <- WCGS$CHD69 # en 0/1     
              
         # recodage de chd     
              
         chd    <- factor(WCGS$CHD69,labels=c("Non","Oui"))     
              
         # régression logistique     
              
         rlb <- glm( formula = chd ~ age, family = binomial() )     
              
         # prédiction     
              
         ypred <- predict(rlb,newdata=data.frame(age),type="response")     
              
         # auroc via le package limma     
              
         library(limma)     
              
         cat(" AUROC (pkg limma) = ", auROC(WCGS$CHD69,ypred),"\n")     
              
         detach("package:limma") ;     
              
         # auroc via le package AUC     
              
         library(AUC)     
              
         cat(" AUROC (pkg AUC)   = ", auc(roc(ypred,chd)),"\n")     
              
         detach("package:AUC") ;     
              
    
    
          AUROC (pkg limma) =  0.6238682     
          AUROC (pkg AUC)   =  0.6199228     
              
    

    Il est donc relativement clair que cette modélisation n'est pas très bonne, puisque l'AUROC ne vaut "que" environ 0,6. A titre de comparaison, l'AUROC de la modélisation précédente, pour les données pg.dar était de 0,985..

    Si on veut disposer de l'intervalle de confiance de l'AUROC, il faut se tourner vers le package pROC avec sa fonction ci.auc() ; ce package fournit également le tracé de la courbe ROC via la fonction plot.roc(). Pour calculer l'odds ratio, le package questionp fournit la fonction odds.ratio(). A l'aide de notre fonction sensSpec() il est possible de connaitre les différents seuils intéressants :

    
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données     
              
         WCGS   <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/wcgs.dar")     
         age    <- WCGS$AGE     
         chdbin <- WCGS$CHD69 # en 0/1     
              
         # recodage de chd     
              
         chd    <- factor(WCGS$CHD69,labels=c("Non","Oui"))     
              
         # régression logistique     
              
         rlb <- glm( formula = chd ~ age, family = binomial() )     
              
         # prédiction     
              
         ypred <- predict(rlb,newdata=data.frame(age),type="response")     
              
         ####################################################     
              
         # odds ratio     
              
         library(questionr)     
         odds.ratio(rlb)     
         detach("package:questionr") ;     
              
         # intervalle de confiance de l'AUROC (package pROC)     
              
         library(pROC)     
         print(ci.auc(chdbin,age))     
              
         # tracé de la courbe ROC (package pROC)     
              
         gr("auroc1.png") ;     
         plot.roc(x=chdbin,predictor=age,print.thres="best",print.auc=TRUE)     
         dev.off()     
              
         detach("package:pROC") ;     
              
         # seuils remarquables avec notre fonction sensSpec     
              
         sensSpec(chdbin,age,FALSE)     
              
         # tracés avec seuils (package PresenceAbsence)     
              
         library(PresenceAbsence)     
              
         dfauc <- data.frame(cbind(1:length(age),chdbin,ypred))     
         gr("auroc2.png") ;     
         auc.roc.plot(dfauc,opt.methods=1:11)     
         dev.off()     
              
         detach("package:PresenceAbsence") ;     
              
    
    
                        OR 2.5 % 97.5 %         p     
         (Intercept) 0.003 0.001  0.008 < 2.2e-16 ***     
         age         1.077 1.054  1.101 4.558e-11 ***     
         ---     
         Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
         95% CI: 0.5832-0.6566 (DeLong)     
              
         Call:     
         plot.roc.default(x = chdbin, predictor = age, print.thres = "best",     print.auc = TRUE)     
              
         Data: age in 2897 controls (chdbin 0) < 257 cases (chdbin 1).     
         Area under the curve: 0.6199     
              
         Seuils remarquables parmi les  24  seuils vus     
              
             CUTOFF   CRIT  VP   FP  FN   VN TOTAL        SPEC      1-SPEC        SENS         VPP         VPN      YOUDEN          OA     
             0.5000 IO     257 2897   0    0  3154    0.000000    1.000000    1.000000    0.081484    1.000000    0.000000    0.081484     
            42.0000 IN     209 1947  48  950  3154    0.327925    0.672075    0.813230    0.096939    0.951904    0.141155    0.367470     
            47.0000 IY     152 1051 105 1846  3154    0.637211    0.362789    0.591440    0.126351    0.946181    0.228651    0.633481     
            58.0000 IP      11   36 246 2861  3154    0.987573    0.012427    0.042802    0.234043    0.920824    0.030375    0.910590     
            59.0000 OA       0    0 257 2897  3154    1.000000    0.000000    0.000000    1.000000    0.918516    0.000000    0.918516     
              
    

              auroc 1          auroc 2

    Une fois le seuil choisi, on peut produire à partir de la RLB une nouvelle classification et la comparer à la classification originale. Notre fonction decritModeleLogistiqueBinaire() se charge d'automatiser tous les calculs. Elle utilise le package ROCR :

    
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données     
              
         WCGS   <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/wcgs.dar")     
         age    <- WCGS$AGE     
         chdbin <- WCGS$CHD69 # en 0/1     
              
         # recodage de chd     
              
         chd    <- factor(WCGS$CHD69,labels=c("Non","Oui"))     
              
         # régression logistique     
              
         rlb <- glm( formula = chd ~ age, family = binomial() )     
              
         # prédiction     
              
         ypred <- predict(rlb,newdata=data.frame(age),type="response")     
              
         ####################################################     
              
         # détails sur la RLB avec la fonction (gH)     
         # qui utilise le package ROCR     
              
         decritModeleLogistiqueBinaire(     
          dfRLB=data.frame(chdbin,age),     
          details=TRUE,     
          sigvar=TRUE,     
          seuil=0.0801     
         ) # fin de decritModeleLogistiqueBinaire     
              
         detach("package:ROCR") ;     
              
    
    
         TRI A PLAT DE :  Cible 0/1 (ordre des modalités)     
              
                              0    1  Total     
          Effectif         2897  257   3154     
          Cumul Effectif   2897 3154   3154     
          Frequence (en %)   92    8    100     
          Cumul fréquences   92  100    100     
              
         VARIABLE :  Cible 0/1 (par fréquence décroissante)     
              
                               0    1  Total     
           Effectif         2897  257   3154     
           Cumul Effectif   2897 3154   3154     
           Frequence (en %)   92    8    100     
           Cumul fréquences   92  100    100     
              
         [1] "0" "1"     
              
              
              
         Call:     
         glm(formula = qlbin ~ ., family = "binomial", data = dfRLB[,     
             -1, drop = FALSE])     
              
         Deviance Residuals:     
             Min       1Q   Median       3Q      Max     
         -0.6209  -0.4545  -0.3669  -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     
              
              
         Auroc (ROCR) =  0.6199228     
                        Estimate Std. Error    z value     Pr(>|z|)       2.5 %      97.5 %    etendue     
         (Intercept) -5.93951594 0.54931839 -10.812520 3.003058e-27 -7.02449218 -4.86960992 2.15488226     
         age          0.07442256 0.01130234   6.584705 4.557900e-11  0.05227351  0.09661197 0.04433846     
              
         Détails par seuil (gH)     
              
         Liste des seuils et de leurs caractéristiques     
              
             CUTOFF   CRIT  VP   FP  FN   VN TOTAL        SPEC      1-SPEC        SENS         VPP         VPN      YOUDEN          OA     
             0.0000        257 2897   0    0  3154    0.000000    1.000000    1.000000    0.081484    1.000000    0.000000    0.081484     
             0.0458        238 2650  19  247  3154    0.085261    0.914739    0.926070    0.082410    0.928571    0.011331    0.153773     
             0.0491        238 2650  19  247  3154    0.085261    0.914739    0.926070    0.082410    0.928571    0.011331    0.153773     
             0.0527        226 2385  31  512  3154    0.176735    0.823265    0.879377    0.086557    0.942910    0.056112    0.233989     
             0.0566 IN     209 1947  48  950  3154    0.327925    0.672075    0.813230    0.096939    0.951904    0.141155    0.367470     
             0.0607        196 1745  61 1152  3154    0.397653    0.602347    0.762646    0.100979    0.949711    0.160299    0.427394     
             0.0651        183 1523  74 1374  3154    0.474284    0.525716    0.712062    0.107268    0.948895    0.186346    0.493659     
             0.0698        171 1349  86 1548  3154    0.534346    0.465654    0.665370    0.112500    0.947368    0.199716    0.545022     
             0.0747        171 1349  86 1548  3154    0.534346    0.465654    0.665370    0.112500    0.947368    0.199716    0.545022     
             0.0801 IY     152 1051 105 1846  3154    0.637211    0.362789    0.591440    0.126351    0.946181    0.228651    0.633481     
             0.0857        152 1051 105 1846  3154    0.637211    0.362789    0.591440    0.126351    0.946181    0.228651    0.633481     
             0.0917        133  906 124 1991  3154    0.687263    0.312737    0.517510    0.128008    0.941371    0.204772    0.673431     
             0.0981        112  793 145 2104  3154    0.726269    0.273731    0.435798    0.123757    0.935527    0.162066    0.702600     
             0.1049         84  563 173 2334  3154    0.805661    0.194339    0.326848    0.129830    0.930993    0.132509    0.766646     
             0.1121         69  465 188 2432  3154    0.839489    0.160511    0.268482    0.129213    0.928244    0.107972    0.792961     
             0.1197         69  465 188 2432  3154    0.839489    0.160511    0.268482    0.129213    0.928244    0.107972    0.792961     
             0.1278         46  276 211 2621  3154    0.904729    0.095271    0.178988    0.142857    0.925494    0.083717    0.845593     
             0.1363         46  276 211 2621  3154    0.904729    0.095271    0.178988    0.142857    0.925494    0.083717    0.845593     
             0.1453         36  206 221 2691  3154    0.928892    0.071108    0.140078    0.148760    0.924107    0.068970    0.864616     
             0.1548         24  142 233 2755  3154    0.950984    0.049016    0.093385    0.144578    0.922021    0.044369    0.881103     
             0.1648 IP      11   36 246 2861  3154    0.987573    0.012427    0.042802    0.234043    0.920824    0.030375    0.910590     
             0.1753 OA       0    0 257 2897  3154    1.000000    0.000000    0.000000    1.000000    0.918516    0.000000    0.918516     
             0.5000 IO       0    0 257 2897  3154    1.000000    0.000000    0.000000    1.000000    0.918516    0.000000    0.918516     
             1.0000          0    0 257 2897  3154    1.000000    0.000000    0.000000    1.000000    0.918516    0.000000    0.918516     
              
         Seuils remarquables parmi les  24  seuils vus     
              
             CUTOFF   CRIT  VP   FP  FN   VN TOTAL        SPEC      1-SPEC        SENS         VPP         VPN      YOUDEN          OA     
             0.0566 IN     209 1947  48  950  3154    0.327925    0.672075    0.813230    0.096939    0.951904    0.141155    0.367470     
             0.0801 IY     152 1051 105 1846  3154    0.637211    0.362789    0.591440    0.126351    0.946181    0.228651    0.633481     
             0.1648 IP      11   36 246 2861  3154    0.987573    0.012427    0.042802    0.234043    0.920824    0.030375    0.910590     
             0.1753 OA       0    0 257 2897  3154    1.000000    0.000000    0.000000    1.000000    0.918516    0.000000    0.918516     
             0.5000 IO       0    0 257 2897  3154    1.000000    0.000000    0.000000    1.000000    0.918516    0.000000    0.918516     
         voici les  2  variables significatives : ","(Intercept)","age     
              
          Analyse des valeurs binaires prédites pour le seuil  0.0801     
         ============================================================     
              
          pct bien classés =  63.34813     
              
         Comparaison binaire org/prédite     
              
         Matrice de pondération     
             C01 C02     
         L01   0   1     
         L02   1   0     
              
         Fréquences initiales de référence =  score1     
         freqref     
            0    1     
         2897  257     
         Fréquences correspondantes calculées = score2     
            0    1     
         1846  152     
              
         Matrice de discordance     
                0    1  Sum     
         0   1846 1051 2897     
         1    105  152  257     
         Sum 1951 1203 3154     
              
         score1 vs score2  score moyen de discordance  0.367   (0.3665187) ; 1998  bien classés sur 3154  soit % BC =  63.3     
          taux des  1998 bien classés sur 3154 :  63.34813     
          taux des  1846 bien classés sur 2897 pour la modalité 0 :  63.72109     
          taux des  152 bien classés sur 257 pour la modalité 1 :  59.14397     
         [1] 0.6199228     
              
    
  8. Il n'y a a priori aucun problème pour modéliser la dépendance de la variable binaire CHD à partir de la variable ARCUS par une régression logistique, mais il est bon de commencer par réaliser le tri croisé de ces variables.

    
         # données     
              
         arcus      <- factor(WCGS$ARCUS,labels=c("Non","Oui"))     
         chd        <- factor(WCGS$CHD69,labels=c("Non","Oui"))     
              
         # tri croisé     
              
         tc <- table(chd,arcus)     
         print(tc)     
         print(chisq.test(tc))     
              
         triCroise("ARCUS",WCGS$ARCUS,"Non Oui","CHD",WCGS$CHD69,"Non oui")     
              
         # régression logistique     
              
         rlb <- glm( formula = chd ~ arcus, family = binomial() )     
         print( summary.glm( rlb ) )     
         print( anova( rlb ) )     
              
         ####################################################     
              
         # auroc     
              
         library(pROC)     
         print(auc(chdbin,as.ordered(arcus)))     
         print(ci.auc(chdbin,as.ordered(arcus)))     
         detach("package:pROC") ;     
              
         # odds ratio     
              
         library(questionr)     
         odds.ratio(rlb)     
         detach("package:questionr") ;     
              
         # risque relatif etc.     
              
         library(epiR)     
         epi.2by2(tc)     
         detach("package:epiR") ;     
              
              
    
    
              arcus     
         chd    Non  Oui     
           Non 2058  839     
           Oui  153  102     
              
         	Pearson's Chi-squared test with Yates' continuity correction     
              
         data:  tc     
         X-squared = 13.1161, df = 1, p-value = 0.0002928     
              
              
         TRI CROISE DES QUESTIONS :     
               ARCUS  (en ligne)     
               CHD  (en colonne)     
         Effectifs     
              Non oui     
         Non 2058 153     
         Oui  839 102     
         Effectifs avec totaux     
                Non oui Total     
         Non   2058 153  2211     
         Oui    839 102   941     
         Total 2897 255  3152     
              
         Valeurs en % du total     
               Non oui TOTAL     
         Non    65   5    70     
         Oui    27   3    30     
         TOTAL  92   8   100     
              
         Call:     
         glm(formula = chd ~ arcus, family = binomial())     
              
         Deviance Residuals:     
             Min       1Q   Median       3Q      Max     
         -0.4790  -0.4790  -0.3787  -0.3787   2.3112     
              
         Coefficients:     
                     Estimate Std. Error z value Pr(>|z|)     
         (Intercept)  -2.5991     0.0838 -31.016  < 2e-16 ***     
         arcusOui      0.4918     0.1342   3.664 0.000248 ***     
         ---     
         Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
              
         (Dispersion parameter for binomial family taken to be 1)     
              
             Null deviance: 1771.2  on 3151  degrees of freedom     
         Residual deviance: 1758.2  on 3150  degrees of freedom     
           (2 observations deleted due to missingness)     
         AIC: 1762.2     
              
         Number of Fisher Scoring iterations: 5     
              
         Analysis of Deviance Table     
              
         Model: binomial, link: logit     
              
         Response: chd     
              
         Terms added sequentially (first to last)     
              
              
               Df Deviance Resid. Df Resid. Dev     
         NULL                   3151     1771.2     
         arcus  1   12.984      3150     1758.2     
         Area under the curve: 0.5552     
         95% CI: 0.524-0.5864 (DeLong)     
                        OR 2.5 % 97.5 %         p     
         (Intercept) 0.074 0.063  0.087 < 2.2e-16 ***     
         arcusOui    1.635 1.254  2.124 0.0002483 ***     
         ---     
         Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
                      Disease +    Disease -      Total        Inc risk *        Odds     
         Exposed +         2058          839       2897              71.0        2.45     
         Exposed -          153          102        255              60.0        1.50     
         Total             2211          941       3152              70.1        2.35     
              
         Point estimates and 95 % CIs:     
         ---------------------------------------------------------     
         Inc risk ratio                         1.18 (1.07, 1.31)     
         Odds ratio                             1.64 (1.24, 2.14)     
         Attrib risk *                          11.04 (4.8, 17.27)     
         Attrib risk in population *            10.15 (3.92, 16.37)     
         Attrib fraction in exposed (%)         15.54 (6.39, 23.8)     
         Attrib fraction in population (%)      14.46 (5.87, 22.28)     
         ---------------------------------------------------------     
          * Cases per 100 population units     
              
    

    Par contre, si on veut utiliser la variable AGEC, il faut commencer par la convertir en facteur afin de pouvoir réaliser la RLB sur chaque modalité.

    
         # données     
              
         age        <- WCGS$AGE     
         agec       <- factor(WCGS$AGEC,labels=paste("Classe",0:4),ordered=TRUE)     
         chd        <- factor(WCGS$CHD69,labels=c("Non","Oui"))     
              
         # tri croisé     
              
         tc <- table(chd,agec)     
         print(tc)     
         print(chisq.test(tc))     
              
         triCroise("AGE (en 5 classes)",WCGS$AGEC,"Classe0 Classe1 Classe2 Classe3 Classe4","CHD",WCGS$CHD69,"Non oui")     
              
         ####################################################     
              
         # régression logistique FAUSSE (age est traitée en quantitative)     
              
         rlb <- glm( formula = chd ~ age, family = binomial() )     
         print( summary.glm( rlb ) )     
         print( anova( rlb ) )     
              
         ####################################################     
              
         # régression logistique CORRECTE (agec est traitée en qualitative)     
              
         rlb <- glm( formula = chd ~ agec, family = binomial() )     
         print( summary.glm( rlb ) )     
         print( anova( rlb ) )     
              
         # auroc     
              
         library(pROC)     
         catn()     
         print(auc(chdbin,agec))     
         print(ci.auc(chdbin,agec))     
         catn()     
              
         # odds ratio     
              
         library(questionr)     
         odds.ratio(rlb)     
         detach("package:questionr") ;     
              
              
    
    
              agec     
         chd   Classe 0 Classe 1 Classe 2 Classe 3 Classe 4     
           Non      512     1036      680      463      206     
           Oui       31       55       70       65       36     
              
         	Pearson's Chi-squared test     
              
         data:  tc     
         X-squared = 46.6534, df = 4, p-value = 1.801e-09     
              
              
         TRI CROISE DES QUESTIONS :     
               AGE (en 5 classes)  (en ligne)     
               CHD  (en colonne)     
         Effectifs     
                  Non oui     
         Classe0  512  31     
         Classe1 1036  55     
         Classe2  680  70     
         Classe3  463  65     
         Classe4  206  36     
         Effectifs avec totaux     
                  Non oui Total     
         Classe0  512  31   543     
         Classe1 1036  55  1091     
         Classe2  680  70   750     
         Classe3  463  65   528     
         Classe4  206  36   242     
         Total   2897 257  3154     
              
         Valeurs en % du total     
                 Non oui TOTAL     
         Classe0  16   1    17     
         Classe1  33   2    35     
         Classe2  22   2    24     
         Classe3  15   2    17     
         Classe4   7   1     8     
         TOTAL    92   8   100     
              
         Call:     
         glm(formula = chd ~ age, family = binomial())     
              
         Deviance Residuals:     
             Min       1Q   Median       3Q      Max     
         -0.6209  -0.4545  -0.3669  -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     
              
         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     
              
         Call:     
         glm(formula = chd ~ agec, family = binomial())     
              
         Deviance Residuals:     
             Min       1Q   Median       3Q      Max     
         -0.5676  -0.4427  -0.3429  -0.3216   2.4444     
              
         Coefficients:     
                     Estimate Std. Error z value Pr(>|z|)     
         (Intercept) -2.34428    0.06908 -33.938  < 2e-16 ***     
         agec.L       0.97791    0.17437   5.608 2.05e-08 ***     
         agec.Q       0.09326    0.16193   0.576   0.5647     
         agec.C      -0.27984    0.14615  -1.915   0.0555 .     
         agec^4       0.16808    0.13208   1.273   0.2032     
         ---     
         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: 1736.3  on 3149  degrees of freedom     
         AIC: 1746.3     
              
         Number of Fisher Scoring iterations: 5     
              
         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     
         agec  4   44.946      3149     1736.3     
              
         Area under the curve: 0.6141     
         95% CI: 0.5783-0.6499 (DeLong)     
              
                        OR 2.5 % 97.5 %         p     
         (Intercept) 0.096 0.084  0.110 < 2.2e-16 ***     
         agec.L      2.659 1.890  3.752 2.045e-08 ***     
         agec.Q      1.098 0.794  1.500   0.56467     
         agec.C      0.756 0.567  1.007   0.05553 .     
         agec^4      1.183 0.912  1.532   0.20317     
         ---     
         Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     
              
    
  9. Exercice volontairement non détaillé sur le site. On fournit juste le code R solution.

    
         # chargement des fonctions gH :     
              
         source("http://forge.info.univ-angers.fr/~gh/wstat/statgh.r",encoding="latin1")     
              
         # lecture des données     
              
         WCGS       <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/wcgs.dar")     
              
         # sélection des données (on en profite pour enlever les NA)     
              
         lesv       <- c("CHD69","AGE","CHOL","SBP","BMI","SMOKE")     
         wcgs       <- na.omit( WCGS[,lesv] )     
              
         # définition de la matrice des résultats     
              
         nbv <- ncol(wcgs) -1     
         mdr <- data.frame(matrix(nrow=nbv,ncol=2))     
         names(mdr)     <- c("AUROC","AIC")     
         row.names(mdr) <- names(wcgs)[-1]     
              
         # on utilise une boucle pour stocker les AUROC     
              
         library(limma) # pour la fonction auROC     
              
         for (idv in (1:nbv)) {     
           rlb <- glm( wcgs$CHD69 ~ wcgs[,1+idv], family="binomial" )     
           # re  <- summary(rlb)     
           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     
              
         mdc(wcgs) # matrice des corrélations linéaires     
              
         #gr("wcgs_pairs1.png")     
         #pairsi(wcgs)     
         #dev.off()     
              
         # meilleur sous-ensemble de régression     
              
         library(MASS) # pour la fonction stepAIC     
              
         mrlb1 <- stepAIC( glm( wcgs$CHD69 ~ . , data=wcgs[,-1], family="binomial" )  , direction="forward" , trace=5)     
         mrlb2 <- stepAIC( glm( wcgs$CHD69 ~ . , data=wcgs[,-1], family="binomial" )  , direction="backward", trace=5)     
         mrlb3 <- stepAIC( glm( wcgs$CHD69 ~ . , data=wcgs[,-1], family="binomial" )  , direction="both"    , trace=5)     
              
         #############################################################################################     
         #     
         # on utilise maintenant 7 variables     
         #     
         #############################################################################################     
              
         catss("Analyse avec 7 variables")     
              
         lesv2       <- c("CHD69","AGE","CHOL","SBP","BMI","SMOKE","ARCUS","DIBPAT")     
         wcgs        <- na.omit( WCGS[,lesv2] )     
              
         mdc(wcgs) # matrice des corrélations linéaires     
              
         nbv <- ncol(wcgs) -1     
         mdr <- data.frame(matrix(nrow=nbv,ncol=2))     
         names(mdr)     <- c("AUROC","AIC")     
         row.names(mdr) <- names(wcgs)[-1]     
              
         # on utilise une boucle pour stocker les AUROC     
              
         library(limma) # pour la fonction auROC     
              
         for (idv in (1:nbv)) {     
           rlb <- glm( wcgs$CHD69 ~ wcgs[,1+idv], family="binomial" )     
           # re  <- summary(rlb)     
           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) )     
              
         mrlb4 <- stepAIC( glm( wcgs$CHD69 ~ . , data=wcgs[,-1], family="binomial" )  , direction="forward" , trace=5)     
         mrlb5 <- stepAIC( glm( wcgs$CHD69 ~ . , data=wcgs[,-1], family="binomial" )  , direction="backward", trace=5)     
         mrlb6 <- stepAIC( glm( wcgs$CHD69 ~ . , data=wcgs[,-1], family="binomial" )  , direction="both"    , trace=5)     
              
         # auroc et AIC pour les 7 variables     
              
         rlb   <- glm( wcgs$CHD69 ~ .,data=wcgs[,-1], family="binomial" )     
         ypred <- predict(rlb,newdata=wcgs,type="response")     
         cat("Auroc globale \n")     
         print( c(auROC( wcgs$CHD69 ,ypred), rlb$aic ))     
              
    
    
              
         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

    
              
         + ======================== +     
         +                          +     
         + Analyse avec 7 variables +     
         +                          +     
         + ======================== +     
              
              
         Matrice des corrélations au sens de  Pearson pour  3140  lignes et  8  colonnes     
              
                CHD69   AGE  CHOL   SBP    BMI SMOKE ARCUS DIBPAT     
         CHD69  1.000     
         AGE    0.120 1.000     
         CHOL   0.161 0.089 1.000     
         SBP    0.134 0.170 0.123 1.000     
         BMI    0.065 0.026 0.071 0.288  1.000     
         SMOKE  0.086 0.003 0.098 0.004 -0.143 1.000     
         ARCUS  0.066 0.188 0.121 0.032 -0.034 0.076 1.000     
         DIBPAT 0.113 0.090 0.057 0.078  0.027 0.061 0.046  1.000     
              
         Meilleures régressions logistiques simples pour CHD69 :     
         =======================================================     
              
                AUROC     AIC     
         CHOL    0.66 1697.84     
         SBP     0.63 1724.10     
         AGE     0.62 1729.77     
         DIBPAT  0.60 1732.08     
         SMOKE   0.58 1750.06     
         BMI     0.57 1760.53     
         ARCUS   0.56 1760.23     
         Start:  AIC=1594.44     
         wcgs$CHD69 ~ AGE + CHOL + SBP + BMI + SMOKE + ARCUS + DIBPAT     
              
         Start:  AIC=1594.44     
         wcgs$CHD69 ~ AGE + CHOL + SBP + BMI + SMOKE + ARCUS + DIBPAT     
              
                  Df Deviance    AIC     
         <none>        1578.4 1594.4     
         - ARCUS   1   1580.8 1594.8     
         - BMI     1   1583.4 1597.4     
         - SBP     1   1596.7 1610.7     
         - SMOKE   1   1596.7 1610.7     
         - AGE     1   1600.1 1614.1     
         - DIBPAT  1   1603.2 1617.2     
         - CHOL    1   1627.0 1641.0     
         Start:  AIC=1594.44     
         wcgs$CHD69 ~ AGE + CHOL + SBP + BMI + SMOKE + ARCUS + DIBPAT     
              
                  Df Deviance    AIC     
         <none>        1578.4 1594.4     
         - ARCUS   1   1580.8 1594.8     
         - BMI     1   1583.4 1597.4     
         - SBP     1   1596.7 1610.7     
         - SMOKE   1   1596.7 1610.7     
         - AGE     1   1600.1 1614.1     
         - DIBPAT  1   1603.2 1617.2     
         - CHOL    1   1627.0 1641.0     
              
         Auroc globale  0.7503843 1594.4446497     
              
    

 

 

retour gH    Retour à la page principale de   (gH)