{"id":719,"date":"2016-10-24T07:27:23","date_gmt":"2016-10-24T07:27:23","guid":{"rendered":"http:\/\/blog.tiran.info\/?p=719"},"modified":"2016-10-24T07:27:23","modified_gmt":"2016-10-24T07:27:23","slug":"regression-logistique-avec-r","status":"publish","type":"post","link":"https:\/\/blog.tiran.stream\/?p=719","title":{"rendered":"R\u00e9gression logistique avec R"},"content":{"rendered":"

Dans la continuit\u00e9 du pr\u00e9c\u00e9dent billet<\/a>, la r\u00e9gression logistique est mise en \u0153uvre avec R cette fois. Les donn\u00e9es utilis\u00e9es sont identiques.<\/p>\n

Chargement et mise en forme des donn\u00e9es<\/h2>\n

A l’aide de ROracle<\/a>, on peuple les datasets jo_ck_prerio (JO pr\u00e9-Rio: \u00e9chantillon d’apprentissage) et jo_ck_rio (JO de Rio: \u00e9chantillon de test):<\/p>\n

> library(ROracle)\n> ora = Oracle()\n> cnx = dbConnect(ora, username="rafa", password="rafa", dbname="S1401037:1521\/STATPDB")\n> \n> jo_ck_prerio <- dbGetQuery(cnx, "select * from JO_CK_PRERIO")\n> jo_ck_rio <- dbGetQuery(cnx, "select * from JO_CK_RIO")\n> \n> dbDisconnect(cnx)\n[1] TRUE\n> \n<\/pre>\n
Les champs TYPE_EPREUVE, SEXE et PAYS sont ensuite convertis en tant que facteurs. On en profite aussi pour fixer la cat\u00e9gorie de r\u00e9f\u00e9rence (Slalom ici) du facteur TYPE_EPREUVE:<\/p>\n
> jo_ck_prerio$TYPE_EPREUVE <- as.factor(jo_ck_prerio$TYPE_EPREUVE)\n> jo_ck_rio$TYPE_EPREUVE <- as.factor(jo_ck_rio$TYPE_EPREUVE)\n> \n> jo_ck_prerio$TYPE_EPREUVE <- relevel(jo_ck_prerio$TYPE_EPREUVE,"Slalom")\n> jo_ck_rio$TYPE_EPREUVE <- relevel(jo_ck_rio$TYPE_EPREUVE,"Slalom")\n> \n> jo_ck_prerio$SEXE <- as.factor(jo_ck_prerio$SEXE)\n> jo_ck_rio$SEXE <- as.factor(jo_ck_rio$SEXE)\n> \n> jo_ck_prerio$PAYS <- as.factor(jo_ck_prerio$PAYS)\n> jo_ck_rio$PAYS <- as.factor(jo_ck_rio$PAYS)\n> \n> summary(jo_ck_prerio)\n     NOM            TYPE_EPREUVE      AGE        SEXE        TAILLE          POIDS               PAYS    \n Length:439         Slalom:163   Min.   :16.00   F:130   Min.   :154.0   Min.   : 50.00   Poland   : 25  \n Class :character   Sprint:276   1st Qu.:23.00   M:309   1st Qu.:172.0   1st Qu.: 66.00   Germany  : 23  \n Mode  :character                Median :26.00           Median :178.0   Median : 75.00   Australia: 20  \n                                 Mean   :25.78           Mean   :177.6   Mean   : 75.19   France   : 19  \n                                 3rd Qu.:28.00           3rd Qu.:184.5   3rd Qu.: 84.00   Slovakia : 18  \n                                 Max.   :42.00           Max.   :202.0   Max.   :109.00   Spain    : 18  \n                                                                                          (Other)  :316  \n>\n<\/pre>\n
Cr\u00e9ation du mod\u00e8le<\/h2>\n
La r\u00e9gression logistique binaire<\/a> est r\u00e9alis\u00e9e \u00e0 l’aide de la fonction glm pour\u00a0laquelle on pr\u00e9cise\u00a0l’argument family=binomial :<\/p>\n
> log_model <- glm(TYPE_EPREUVE ~ . - NOM,family=binomial,data=jo_ck_prerio)\n> summary(log_model)\n\nCall:\nglm(formula = TYPE_EPREUVE ~ . - NOM, family = binomial, data = jo_ck_prerio)\n\nDeviance Residuals: \n     Min        1Q    Median        3Q       Max  \n-2.38343  -0.43245   0.00008   0.42094   2.46605  \n\nCoefficients:\n                             Estimate Std. Error z value Pr(>|z|)    \n(Intercept)                -2.063e+01  1.075e+04  -0.002  0.99847    \nAGE                        -1.785e-03  3.923e-02  -0.045  0.96372    \nSEXEM                      -4.230e+00  6.332e-01  -6.680 2.38e-11 ***\nTAILLE                     -1.071e-01  3.840e-02  -2.790  0.00527 ** \nPOIDS                       3.649e-01  4.598e-02   7.935 2.11e-15 ***\nPAYSAngola                  3.770e+01  1.312e+04   0.003  0.99771    \nPAYSArgentina               1.807e+01  1.075e+04   0.002  0.99866    \nPAYSAustralia               1.549e+01  1.075e+04   0.001  0.99885    \nPAYSAustria                 1.508e+01  1.075e+04   0.001  0.99888    \nPAYSAzerbaijan              3.123e+01  1.521e+04   0.002  0.99836    \nPAYSBelarus                 3.333e+01  1.151e+04   0.003  0.99769    \nPAYSBelgium                 1.775e+01  1.075e+04   0.002  0.99868    \nPAYSBosnia and Herzegovina -2.293e+00  1.521e+04   0.000  0.99988    \nPAYSBrazil                  1.775e+01  1.075e+04   0.002  0.99868    \nPAYSBulgaria                3.459e+01  1.148e+04   0.003  0.99759    \nPAYSCanada                  1.635e+01  1.075e+04   0.002  0.99879    \nPAYSChile                   1.502e+01  1.075e+04   0.001  0.99889    \nPAYSChina                   1.585e+01  1.075e+04   0.001  0.99882    \nPAYSCook Islands            1.072e+01  1.075e+04   0.001  0.99920    \nPAYSCote d'Ivoire           3.460e+01  1.521e+04   0.002  0.99818    \nPAYSCroatia                 1.397e+01  1.075e+04   0.001  0.99896    \nPAYSCuba                    3.522e+01  1.171e+04   0.003  0.99760    \nPAYSCzech Republic          1.522e+01  1.075e+04   0.001  0.99887    \nPAYSDenmark                 3.350e+01  1.180e+04   0.003  0.99774    \nPAYSEcuador                 3.018e+01  1.521e+04   0.002  0.99842    \nPAYSEgypt                   3.867e+01  1.521e+04   0.003  0.99797    \nPAYSEstonia                 3.445e+01  1.521e+04   0.002  0.99819    \nPAYSFinland                 3.491e+01  1.233e+04   0.003  0.99774    \nPAYSFrance                  1.629e+01  1.075e+04   0.002  0.99879    \nPAYSGermany                 1.472e+01  1.075e+04   0.001  0.99891    \nPAYSGreat Britain           1.534e+01  1.075e+04   0.001  0.99886    \nPAYSGreece                  1.569e+01  1.075e+04   0.001  0.99884    \nPAYSGuam                    3.637e+01  1.521e+04   0.002  0.99809    \nPAYSHungary                 3.423e+01  1.110e+04   0.003  0.99754    \nPAYSIndonesia               3.509e+01  1.521e+04   0.002  0.99816    \nPAYSIran                    3.708e+01  1.199e+04   0.003  0.99753    \nPAYSIreland                 1.538e+01  1.075e+04   0.001  0.99886    \nPAYSIsrael                  3.390e+01  1.219e+04   0.003  0.99778    \nPAYSItaly                   1.675e+01  1.075e+04   0.002  0.99876    \nPAYSJapan                   1.685e+01  1.075e+04   0.002  0.99875    \nPAYSKazakhstan              1.840e+01  1.075e+04   0.002  0.99864    \nPAYSLatvia                  3.435e+01  1.225e+04   0.003  0.99776    \nPAYSLithuania               3.352e+01  1.130e+04   0.003  0.99763    \nPAYSMacedonia              -1.953e+00  1.282e+04   0.000  0.99988    \nPAYSMexico                  3.276e+01  1.296e+04   0.003  0.99798    \nPAYSMorocco                -3.109e+00  1.304e+04   0.000  0.99981    \nPAYSMyanmar                 3.630e+01  1.521e+04   0.002  0.99810    \nPAYSNetherlands            -4.471e+00  1.166e+04   0.000  0.99969    \nPAYSNew Zealand             1.684e+01  1.075e+04   0.002  0.99875    \nPAYSNigeria                -4.972e+00  1.521e+04   0.000  0.99974    \nPAYSNorway                  3.477e+01  1.216e+04   0.003  0.99772    \nPAYSPoland                  1.586e+01  1.075e+04   0.001  0.99882    \nPAYSPortugal                1.663e+01  1.075e+04   0.002  0.99877    \nPAYSRomania                 3.376e+01  1.120e+04   0.003  0.99760    \nPAYSRussia                  1.669e+01  1.075e+04   0.002  0.99876    \nPAYSSamoa                   3.304e+01  1.521e+04   0.002  0.99827    \nPAYSSao Tome and Principe   3.767e+01  1.521e+04   0.002  0.99802    \nPAYSSenegal                 3.653e+01  1.310e+04   0.003  0.99777    \nPAYSSerbia                  3.314e+01  1.162e+04   0.003  0.99772    \nPAYSSeychelles              3.297e+01  1.521e+04   0.002  0.99827    \nPAYSSingapore               3.412e+01  1.521e+04   0.002  0.99821    \nPAYSSlovakia                1.552e+01  1.075e+04   0.001  0.99885    \nPAYSSlovenia                1.509e+01  1.075e+04   0.001  0.99888    \nPAYSSouth Africa            1.869e+01  1.075e+04   0.002  0.99861    \nPAYSSouth Korea             3.658e+01  1.285e+04   0.003  0.99773    \nPAYSSpain                   1.590e+01  1.075e+04   0.001  0.99882    \nPAYSSweden                  3.487e+01  1.137e+04   0.003  0.99755    \nPAYSSwitzerland             1.512e+01  1.075e+04   0.001  0.99888    \nPAYSThailand               -1.510e+00  1.521e+04   0.000  0.99992    \nPAYSTogo                   -2.431e+00  1.521e+04   0.000  0.99987    \nPAYSTunisia                 3.535e+01  1.183e+04   0.003  0.99762    \nPAYSUkraine                 3.473e+01  1.200e+04   0.003  0.99769    \nPAYSUnited States           1.492e+01  1.075e+04   0.001  0.99889    \nPAYSUruguay                 3.554e+01  1.521e+04   0.002  0.99814    \nPAYSUzbekistan              3.422e+01  1.216e+04   0.003  0.99775    \nPAYSVenezuela               3.137e+01  1.521e+04   0.002  0.99835    \nPAYSVietnam                 3.778e+01  1.521e+04   0.002  0.99802    \n---\nSignif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1\n\n(Dispersion parameter for binomial family taken to be 1)\n\n    Null deviance: 579.17  on 438  degrees of freedom\nResidual deviance: 267.27  on 362  degrees of freedom\nAIC: 421.27\n\nNumber of Fisher Scoring iterations: 18\n\n>\n<\/pre>\n
A l\u2019instar de la fonctionnalit\u00e9 de \u00ab feature selection \u00bb d\u2019ODM, on peut utiliser la fonction step<\/a> pour r\u00e9aliser une d\u00e9termination pas \u00e0 pas de la meilleure combinaison de pr\u00e9dicteurs:<\/p>\n
> step(log_model, dir="backward")\nStart:  AIC=421.27\nTYPE_EPREUVE ~ (NOM + AGE + SEXE + TAILLE + POIDS + PAYS) - NOM\n\n         Df Deviance    AIC\n- PAYS   72   407.80 417.80\n- AGE     1   267.28 419.28\n<none>        267.27 421.27\n- TAILLE  1   275.58 427.58\n- SEXE    1   329.21 481.21\n- POIDS   1   376.09 528.09\n\nStep:  AIC=417.8\nTYPE_EPREUVE ~ AGE + SEXE + TAILLE + POIDS\n\n         Df Deviance    AIC\n- AGE     1   407.93 415.93\n<none>        407.80 417.80\n- TAILLE  1   422.42 430.42\n- SEXE    1   474.27 482.27\n- POIDS   1   541.79 549.79\n\nStep:  AIC=415.93\nTYPE_EPREUVE ~ SEXE + TAILLE + POIDS\n\n         Df Deviance    AIC\n<none>        407.93 415.93\n- TAILLE  1   422.77 428.77\n- SEXE    1   474.39 480.39\n- POIDS   1   541.81 547.81\n\nCall:  glm(formula = TYPE_EPREUVE ~ SEXE + TAILLE + POIDS, family = binomial, \n    data = jo_ck_prerio)\n\nCoefficients:\n(Intercept)        SEXEM       TAILLE        POIDS  \n    -0.6010      -3.1906      -0.1019       0.2901  \n\nDegrees of Freedom: 438 Total (i.e. Null);  435 Residual\nNull Deviance:\t    579.2 \nResidual Deviance: 407.9 \tAIC: 415.9\n> \n<\/pre>\n
La fonction step parvient \u00e0 la conclusion que seuls les pr\u00e9dicteurs SEXE, TAILLE et POIDS doivent \u00eatre maintenus dans le mod\u00e8le final:<\/p>\n
> log_model <- glm(TYPE_EPREUVE ~ TAILLE + SEXE + POIDS,family=binomial,data=jo_ck_prerio)\n> summary(log_model)\n\nCall:\nglm(formula = TYPE_EPREUVE ~ TAILLE + SEXE + POIDS, family = binomial, \n    data = jo_ck_prerio)\n\nDeviance Residuals: \n    Min       1Q   Median       3Q      Max  \n-2.8981  -0.7873   0.3187   0.7275   2.2077  \n\nCoefficients:\n            Estimate Std. Error z value Pr(>|z|)    \n(Intercept) -0.60105    3.59533  -0.167 0.867233    \nTAILLE      -0.10189    0.02772  -3.675 0.000238 ***\nSEXEM       -3.19061    0.44163  -7.225 5.03e-13 ***\nPOIDS        0.29005    0.03179   9.123  < 2e-16 ***\n---\nSignif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1\n\n(Dispersion parameter for binomial family taken to be 1)\n\n    Null deviance: 579.17  on 438  degrees of freedom\nResidual deviance: 407.93  on 435  degrees of freedom\nAIC: 415.93\n\nNumber of Fisher Scoring iterations: 5\n\n>\n<\/pre>\n
Scoring du mod\u00e8le<\/h2>\n
La fonction predict permet d’appliquer notre mod\u00e8le sur le dataframe de test jo_ck_rio:<\/p>\n
> rio_pred <- predict(log_model, jo_ck_rio)\n> head(rio_pred)\n         1          2          3          4          5          6 \n-1.9928707 -1.9928707 -1.5458741 -1.9534234 -0.9501602 -1.0520475 \n> \n<\/pre>\n
Les r\u00e9sultats en sortie correspondent aux valeurs rendues par la fonction logistique<\/a> utilisant les coefficients obtenus plus haut:<\/p>\n
\"P<\/p>\n
Ici, la variable d\u00e9pendante est un facteur et l’aide de la fonction glm (?glm) permet de trouver la convention utilis\u00e9e par R pour d\u00e9finir quel niveau est associ\u00e9 avec le succ\u00e8s 1 ou l’\u00e9chec 0:<\/p>\n
\u00ab\u00a0For binomial and quasibinomial families the response can also be specified as a factor (when the first level<\/strong> denotes failure and all others success)\u00a0\u00bb<\/em><\/p>\n
Ici, le premier niveau du facteur TYPE_EPREUVE est Slalom:<\/p>\n
> levels(jo_ck_prerio$TYPE_EPREUVE)[1]\n[1] "Slalom"\n> levels(jo_ck_rio$TYPE_EPREUVE)[1]\n[1] "Slalom"\n> \n<\/pre>\n
Le point d’inflexion de la fonction logistique \u00e9tant \u00e0 0.5, on consid\u00e9rera que les valeurs sup\u00e9rieures correspondent \u00e0 une pr\u00e9diction de la cat\u00e9gorie Sprint.<\/p>\n
> rio_pred.bin <- ifelse(rio_pred > 0.5,"Sprint","Slalom")\n>\n<\/pre>\n
On peut alors produire une matrice de confusion afin de mesurer la qualit\u00e9 de nos pr\u00e9dictions:<\/p>\n
> table(rio_pred.bin, jo_ck_rio$TYPE_EPREUVE)\n            \nrio_pred.bin Slalom Sprint\n      Slalom     46     15\n      Sprint     15    108\n> \n<\/pre>\n
Le taux de r\u00e9ussite des pr\u00e9dictions avec ce mod\u00e8le logistique binaire est de l’ordre de 84%<\/strong> (46+108)\/(46+15+15+108).<\/p>\n
Repr\u00e9sentation graphique<\/h2>\n
La repr\u00e9sentation graphique des donn\u00e9es permet de mieux appr\u00e9hender l’origine des d\u00e9fauts de classification.\u00a0On remarque bien\u00a0au passage que les comp\u00e9titeurs des deux disciplines ont globalement des corpulences diff\u00e9rentes:<\/p>\n
> prerio <- jo_ck_prerio[,c("TAILLE","POIDS","SEXE","TYPE_EPREUVE")]\n> prerio$COLOR <- ifelse(prerio$TYPE_EPREUVE == "Slalom", "blue","forestgreen")\n> prerio$PCH <- ifelse(prerio$TYPE_EPREUVE == "Slalom", 3,4)\n> prerio <- prerio[,c("TAILLE","POIDS","SEXE","COLOR","PCH")]\n> \n> rio <- jo_ck_rio[,c("TAILLE","POIDS","SEXE","TYPE_EPREUVE")]\n> rio$COLOR <- ifelse(rio$TYPE_EPREUVE == "Slalom", "blue","forestgreen")\n> rio$PCH <- ifelse(rio$TYPE_EPREUVE == "Slalom", 3,4)\n> rio$pred <- as.factor(rio_pred.bin)\n> \n> rio$COLOR[which(rio$pred!=rio$TYPE_EPREUVE)] <- "red"\n> rio <- rio[,c("TAILLE","POIDS","SEXE","COLOR","PCH")]\n> \n> jo_ck <- rbind(prerio,rio)\n> \n> library(scatterplot3d)\n> library(extrafont)\nRegistering fonts with R\n> loadfonts(device="win")\nAgency FB already registered with windowsFonts().\n[...]\nWingdings 3 already registered with windowsFonts().\n> \n>\n> par(family = "Verdana")\n>\n> scatterplot3d(jo_ck$POIDS, jo_ck$TAILLE, jo_ck$SEXE,\n+ main = "Comp\u00e9titeurs Individuels Canoe-Kayak au JO depuis 2000",\n+ xlab = "Poids (kg)",\n+ ylab = "Taille (cm)",\n+ zlab = "",\n+ color = jo_ck$COLOR,\n+ angle = 25,\n+ pch=jo_ck$PCH,\n+ cex.symbols=0.6,\n+ z.ticklabs=c("Femmes","Hommes"),\n+ lab.z=1,\n+ cex.lab=1.2,\n+ font.lab=2)\n>\n>\n> legend("topleft",\n+ bty="n",\n+ cex=0.8,\n+ c("Slalomeur", "Sprinteur", "Slalomeur mal class\u00e9", "Sprinteur mal class\u00e9"),\n+ col=c("blue", "forestgreen", "red","red"),\n+ pch=c(3,4,3,4))\n> \n<\/pre>\n
 <\/p>\n
\"jo_ck_logit\"<\/a><\/p>\n
 <\/p>\n","protected":false},"excerpt":{"rendered":"
Dans la continuit\u00e9 du pr\u00e9c\u00e9dent billet, la r\u00e9gression logistique est mise en \u0153uvre avec R cette fois. Les donn\u00e9es utilis\u00e9es<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"colormag_page_container_layout":"default_layout","colormag_page_sidebar_layout":"default_layout","footnotes":""},"categories":[3,12,13,17],"tags":[],"class_list":["post-719","post","type-post","status-publish","format-standard","hentry","category-classification","category-r","category-regression","category-visualisation"],"_links":{"self":[{"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=\/wp\/v2\/posts\/719","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=719"}],"version-history":[{"count":0,"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=\/wp\/v2\/posts\/719\/revisions"}],"wp:attachment":[{"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=719"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=719"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blog.tiran.stream\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=719"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}