Stata : modélisation statistique (1)

Table des matières

Dans cet exposé, on s’intéressera à la construction de modèles de régression explicatifs ou prédictifs. La première partie se concentre sur le modèle linéaire et ses applications. La seconde partie porte sur le modèle liénaire généralisé. Il existe un très bon ouvrage traitant du modèle linéaire généralisé, à présent dans sa quatrième édition et publié chez Stata Press : Generalized Linear Models and Extensions de Hardin & Hilbe. Concernant la modélisation statistique en général, les ouvrages de Vittinghoff et coll. [6] et Dupont [2] couvrent l’essentiel des notions présentées. Pour aller plus loin, le livre de Harrell [3] demeure une référence en ce qui concerne les techniques biostatistiques modernes appliquées dans un contexte médical.

La mise en oeuvre d’un modèle de régression a déjà été discutée brièvement dans le tutoriel d’introduction à Stata. Dans ce chapitre, on va s’intéresser à l’estimation des paramètres d’un modèle de régression linéaire, à la sélection de modèle dans un cadre explicatif, au diagnostic du modèle, et à la prédiction ponctuelle ou par intervalles. On prendra pour base des données observationnelles issues d’enquêtes ou d’études cliniques transversales. Les séries chronologiques et les données longitudinales seront traitées dans des chapitres séparés.

Le modèle de régression linéaire simple

Dans un premier temps, procédons à quelques rappels concernant la régression linéaire simple, la corrélation linéaire et le test de Student. Les notions connexes telles que les associations non linéaires ou les approches non paramétriques seront traitées ultérieurement.

Un exemple de régression linéaire simple

Les données d’illustration peuvent être chargées directement depuis internet à l’aide de la commande webuse. Il s’agit d’une enquête épidémiologique rétrospective dans laquelle on s’intéresse aux facteurs de risque d’un bébé ayant un poids inférieur à la norme, selon les normes américaines des années 90. Ces données sont extensivement analysées dans l’ouvrage de Hosmer & Lemeshow [4]. Au total, le tableau de données contient 11 variables, la variable réponse étant le poids des bébés (bwt) :

webuse lbw
describe, simple
set more off
webuse lbw
(Hosmer & Lemeshow data)
describe, simple
id     low    age    lwt    race   smoke  ptl    ht     ui     ftv    bwt

La relation entre le poids des bébés (en grammes) et le poids des mères (initialement en livres, converti en kilogrammes) est représentée dans la figure suivante sous la forme d’un simple diagramme de dispersion. Pour faciliter la lecture, le poids des mères est converti en kilogrammes, sans modifier le mode de représentation de la variable :

replace lwt = lwt / 2.204623, nopromote
replace lwt = lwt / 2.204623, nopromote
(189 real changes made)
(189 values truncated because of storage type)

Ensuite, on utilise la commande scatter pour construire le diagramme de dispersion :

set scheme plotplain
graph twoway scatter bwt lwt
graph export "fig-03-scatter-bwt-lwt.pdf", fontface(DroidSans) replace

fig-03-scatter-bwt-lwt.png

Figure 1 : Relation entre le poids des bébés et le poids des mères

Il est tout à fait possible et largement recommendé de superposer une courbe loess [1] sur le diagramme de dispersion précédent afin d’évaluer visuellement les écarts à la linéarité concernant la relation entre ces deux variables. La commande lowess peut être combinée à scatter :

graph twoway scatter bwt lwt || lowess bwt lwt, legend(off)
graph export "fig-03-lowess-bwt-lwt.pdf", fontface(DroidSans) replace

fig-03-lowess-bwt-lwt.png

Figure 2 : Relation entre le poids des bébés et le poids des mères (courbe loess)

La corrélation entre ces deux variables s’obtient grâce à correlate. Notons que cette commande fonctionne avec deux variables ou une liste de variables de sorte qu’elle pourra également être utilisée pour construire une matrice de corrélation :

summarize bwt lwt
correlate bwt lwt
summarize bwt lwt

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
         bwt |        189    2944.286     729.016        709       4990
         lwt |        189    58.34392    13.88108         36        113
correlate bwt lwt
(obs=189)

             |      bwt      lwt
-------------+------------------
         bwt |   1.0000
         lwt |   0.1863   1.0000

Voici une formulation simplifiée du modèle de régression linéaire. Soit \(y_i\) la réponse observée sur l’individu \(i\), et \(x_i\) sa valeur observée pour le prédicteur \(x\). Le modèle de régression linéaire s’écrit : \[y_i = \beta_0+\beta_1x_i+\varepsilon_i,\] où \(\beta_0\) représente l’ordonnée à l’origine (intercept) et \(\beta_1\) la pente (\emph{slope}) de la droite de régression, et \(\varepsilon_i\sim\mathcal{N}(0,\sigma^2)\) est un terme d’erreur (résidus, supposés indépendants). En minimisant les différences quadratiques entre les valeurs observées et les valeurs prédites (principe des MCO), on peut estimer les coefficients de régression, \(\widehat\beta_0\) et \(\widehat\beta_1\) : \[\begin{array}{l} \widehat\beta_0 = \bar y - \widehat\beta_1\bar x\\ \widehat\beta_1 = \sum(y_i-\bar y)(x_i-\bar x)/\sum(x_i-\bar x)^2\\ \end{array}\]

Sous \(H_0\), le rapport entre l’estimé de la pente (\(\widehat\beta_1\), de variance \(\frac{\text{SSR}/(n-2)}{(n-1)s_x^2}\)) et son erreur standard suit une loi de Student à \((n-2)\) degrés de liberté.

Les paramètres d’un tel modèle de régression, \(\widehat\beta_0\) et \(\widehat\beta_1\), peuvent être estimés grâce à la commande regress, en indiquant la variable à prédire et la ou les variables explicatives. Pour un modèle de régression linéaire simple, on se retrouve donc avec l’expression la plus simple qui soit :

regress bwt lwt
regress bwt lwt

      Source |       SS           df       MS      Number of obs   =       189
-------------+----------------------------------   F(1, 187)       =      6.72
       Model |  3467517.36         1  3467517.36   Prob > F        =    0.0103
    Residual |  96447781.2       187  515763.536   R-squared       =    0.0347
-------------+----------------------------------   Adj R-squared   =    0.0295
       Total |  99915298.6       188  531464.354   Root MSE        =    718.17

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         lwt |   9.783793   3.773317     2.59   0.010     2.340054    17.22753
       _cons |   2373.461    226.263    10.49   0.000     1927.105    2819.817
------------------------------------------------------------------------------

Les résultats fournis par regress se composent de deux tableaux : le tableau d’analyse de variance du modèle de régression, qui peut être supprimé via l’option noheader, et le tableau des coefficients de régression. La ligne _cons désigne le terme d’ordonnée à l’origine, \(\widehat\beta_0\). Le résultat du test de Student associé à lwt (\(\widehat\beta_1\)) peut se retrouver manuellement une fois que l’on a extrait les valeurs d’intérêt :

local tstat = _b[lwt] / _se[lwt]
display "t = " %4.2f `tstat' " p = " %4.3f 2*ttail(187, `tstat')
local tstat = _b[lwt] / _se[lwt]
display "t = " %4.2f `tstat' " p = " %4.3f 2*ttail(187, `tstat')
t = 2.59 p = 0.010

Diagnostic du modèle

La commande predict permet non seulement de calculer les valeurs ajustées du modèle mais également les résidus du modèle (\(e_i = \tilde y_i - y_i\)) ainsi que d’autres statistiques utiles pour diagnostiquer la qualité d’ajustement du modèle de régression.

predict double yhat
predict double rs, rstudent
summarize rs
predict double yhat
(option xb assumed; fitted values)
predict double rs, rstudent
summarize rs

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
          rs |        189   -.0010581    1.008052  -3.133601   2.961918

Dans le cas ci-dessus, ce sont les résidus studentisés, \(r_i = e_i / (s_{(i)}\sqrt{1-h_i})\), qui ont été calculés. D’autres options sont également disponibles mais ce type de résidus facilite à la fois l’interprétation et la détection de valeurs extrêmes (voir également l’aide en ligne, help regress postestimation). Par exemple, voici un diagramme de quantiles pour les résidus simples :

predict double r, resid
qnorm r
graph export "fig-03-qnorm-r.pdf", fontface(DroidSans) replace

fig-03-qnorm-r.png

Figure 3 : Distribution des résidus studentisés

Un histogramme ou une courbe de densité permet également d’examiner rapidement la forme de la distribution des résidus. Voici un exemple avec kdensity, pour lequel une courbe de densité normale a été ajoutée à l’aide de l’option normal :

kdensity rs, normal normopts(lpat(--))
graph export "fig-03-kdensity-rs.pdf", fontface(DroidSans) replace

fig-03-kdensity-rs.png

Figure 4 : Distribution des résidus studentisés

Le graphique suivant est plus informatif car il renseigne à la fois sur la distribution des résidus et la corrélation entre les valeurs prédites par le modèle et ces derniers, qui selon l’hypothèse du modèle doit être nulle. Ici, on utilise les valeurs de post-estimation calculées plus haut, mais il serait tout à fait possible d’utiliser directement la commande de post-estimation rvfplot (ou rvpplot, mais dans le cas d’une régression avec un seul prédicteur cela ne change rien) :

graph twoway scatter rs yhat, yline(0)
graph export "fig-03-scatter-rs-yhat.pdf", fontface(DroidSans) replace

fig-03-scatter-rs-yhat.png

Figure 5 : Relation entre valeurs ajustées et résidus

Cas de la régression sur une variable catégorielle

On a vu dans le chapitre sur la gestion des données comment représenter les variables catégorielles sous Stata : dans le cas des variables binaires, un codage sous forme de 0 et de 1 est parfaitement adéquat, tandis que dans le cas des variables à plus de deux modalités, on assigne à chaque niveau un code numérique en débutant à 1. Ainsi, pour une variable à trois modalités, le premier niveau sera représenté par la valeur 1 tandis que le troisième et dernier niveau prendra la valeur 3. On associera éventuellement des étiquettes à chacun des niveaux afin de mieux identifier les différentes classes.

Considérons la variable smoke qui indique si la mère fumait pendant le premier trimestre de sa grossesse :

tabulate smoke, nolabel
tabstat bwt, by(smoke) stat(mean sd n)
tabulate smoke, nolabel

     smoked |
     during |
  pregnancy |      Freq.     Percent        Cum.
------------+-----------------------------------
          0 |        115       60.85       60.85
          1 |         74       39.15      100.00
------------+-----------------------------------
      Total |        189      100.00
tabstat bwt, by(smoke) stat(mean sd n)

Summary for variables: bwt
     by categories of: smoke (smoked during pregnancy)

    smoke |      mean        sd         N
----------+------------------------------
nonsmoker |  3054.957   752.409       115
   smoker |  2772.297  659.8075        74
----------+------------------------------
    Total |  2944.286   729.016       189
-----------------------------------------
graph box bwt, over(smoke)
graph export "fig-03-box-bwt-smoke.pdf", fontface(DroidSans) replace

fig-03-box-bwt-smoke.png

Figure 6 : Relation entre poids des bébés et statut fumeur

Le modèle de régression suivant considère la variable smoke comme une variable numérique et le coefficient de régression pour cette variable représente la variation de poids lorsque smoke varie d’une unité (de 0 à 1) :

regress bwt smoke
regress bwt smoke

      Source |       SS           df       MS      Number of obs   =       189
-------------+----------------------------------   F(1, 187)       =      6.98
       Model |  3597444.33         1  3597444.33   Prob > F        =    0.0089
    Residual |  96317854.2       187  515068.739   R-squared       =    0.0360
-------------+----------------------------------   Adj R-squared   =    0.0308
       Total |  99915298.6       188  531464.354   Root MSE        =    717.68

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       smoke |  -282.6592   106.9544    -2.64   0.009    -493.6515   -71.66693
       _cons |   3054.957   66.92428    45.65   0.000     2922.933     3186.98
------------------------------------------------------------------------------

En indiquant à Stata que la variable smoke doit être traitée comme une variable catégorielle et de générer l’ensemble de variables indicatrices correspondant, on obtiendra strictement le même résultat du fait du codage initial en 0/1 :

regress bwt i.smoke
regress bwt i.smoke

      Source |       SS           df       MS      Number of obs   =       189
-------------+----------------------------------   F(1, 187)       =      6.98
       Model |  3597444.33         1  3597444.33   Prob > F        =    0.0089
    Residual |  96317854.2       187  515068.739   R-squared       =    0.0360
-------------+----------------------------------   Adj R-squared   =    0.0308
       Total |  99915298.6       188  531464.354   Root MSE        =    717.68

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       smoke |
     smoker  |  -282.6592   106.9544    -2.64   0.009    -493.6515   -71.66693
       _cons |   3054.957   66.92428    45.65   0.000     2922.933     3186.98
------------------------------------------------------------------------------

Considérons la variable race qui a trois niveaux. Il est tout à fait possible de générer l’ensemble des indicatrices associées à cette variable à l’aide de tabulate :

quietly tabulate race, gen(irace)
list race irace* in 1/5
quietly tabulate race, gen(irace)
list race irace* in 1/5

     +----------------------------------+
     |  race   irace1   irace2   irace3 |
     |----------------------------------|
  1. | black        0        1        0 |
  2. | other        0        0        1 |
  3. | white        1        0        0 |
  4. | white        1        0        0 |
  5. | white        1        0        0 |
     +----------------------------------+

Ensuite, il suffira d’inclure deux indicatrices parmi les trois dans le modèle de régression, par exemple regress bwt irace2 irace3. Mais comme on l’a vu plus haut, l’opérateur i. permet de générer automatiquement un ensemble d’indicatrices pour n’importe quelle variable catégorielle :

regress bwt i.race
regress bwt i.race

      Source |       SS           df       MS      Number of obs   =       189
-------------+----------------------------------   F(2, 186)       =      4.95
       Model |  5048361.06         2  2524180.53   Prob > F        =    0.0081
    Residual |  94866937.5       186  510037.298   R-squared       =    0.0505
-------------+----------------------------------   Adj R-squared   =    0.0403
       Total |  99915298.6       188  531464.354   Root MSE        =    714.17

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        race |
      black  |  -383.3181   157.8914    -2.43   0.016    -694.8064   -71.82985
      other  |  -298.9955   113.6899    -2.63   0.009    -523.2829   -74.70811
             |
       _cons |    3103.01   72.88956    42.57   0.000     2959.214    3246.807
------------------------------------------------------------------------------

Par défaut, le premier niveau de la variable catégorielle (ici, white) sert de niveau de référence, mais il est tout à fait possible de modifier ce comportement en indiquant la catégorie à utiliser pour la « baseline ». En utilisant le préfixe ib3, on indique à Stata que le troisième niveau de race servira de catégorie de référence :

regress bwt ib3.race
regress bwt ib3.race

      Source |       SS           df       MS      Number of obs   =       189
-------------+----------------------------------   F(2, 186)       =      4.95
       Model |  5048361.06         2  2524180.53   Prob > F        =    0.0081
    Residual |  94866937.5       186  510037.298   R-squared       =    0.0505
-------------+----------------------------------   Adj R-squared   =    0.0403
       Total |  99915298.6       188  531464.354   Root MSE        =    714.17

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        race |
      white  |   298.9955   113.6899     2.63   0.009     74.70811    523.2829
      black  |  -84.32262   165.0131    -0.51   0.610    -409.8604    241.2152
             |
       _cons |   2804.015   87.24962    32.14   0.000     2631.889    2976.141
------------------------------------------------------------------------------

On retrouvera bien les différences de moyennes par simple estimation de contrastes grâce à contrast ou margins :

contrast r.race, nowald effects
contrast r.race, nowald effects

Contrasts of marginal linear predictions

Margins      : asbalanced

------------------------------------------------------------------------------
             |   Contrast   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        race |
     (white  |
         vs  |
     other)  |   298.9955   113.6899     2.63   0.009     74.70811    523.2829
     (black  |
         vs  |
     other)  |  -84.32262   165.0131    -0.51   0.610    -409.8604    241.2152
------------------------------------------------------------------------------

Lien avec le test de Student

La différence de moyennes utilisée pour former la statistique de test de Student et qui est rappelée dans la sortie de ttest ci-dessous correspond strictement à la pente de la droite de régression estimée dans la section précédente :

ttest bwt, by(smoke)
ttest bwt, by(smoke)

Two-sample t test with equal variances
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
nonsmoke |     115    3054.957     70.1625     752.409    2915.965    3193.948
  smoker |      74    2772.297    76.70106    659.8075    2619.432    2925.162
---------+--------------------------------------------------------------------
combined |     189    2944.286    53.02811     729.016    2839.679    3048.892
---------+--------------------------------------------------------------------
    diff |            282.6592    106.9544                71.66693    493.6515
------------------------------------------------------------------------------
    diff = mean(nonsmoke) - mean(smoker)                          t =   2.6428
Ho: diff = 0                                     degrees of freedom =      187

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.9955         Pr(|T| > |t|) = 0.0089          Pr(T > t) = 0.0045

On peut d’ailleurs visualiser très facilement ce différentiel de moyennes à l’aide d’un simple diagramme de dispersion en considérant la variable binaire sur l’axe des abscisses. Plutôt que d’utiliser scatter et de redéfinir l’axe des x, il est plus simple d’utiliser un diagramme un point tel que proposé par la commande externe stripplot (à installer au préalable, ssc install stripplot) :

stripplot bwt, over(smoke) vertical jitter(1 0) addplot(lfit bwt smoke)
graph export "fig-03-stripplot-bwt-smoke.pdf", fontface(DroidSans) replace

fig-03-stripplot-bwt-smoke.png

Figure 7 : Relation entre poids des bébés et statut fumeur

Une manière de vérifier graphiquement l’hypothèse d’égalité des variances, nécessaire dans le test ci-dessus afin de recouvrer les résultats du test du coefficient de régression, consisterait à comparer les fonctions de répartition empirique des deux groupes comme suggéré sur le forum Stata.

Dans le cas d’une variable catégorielle à plus de deux niveaux, telle que race, il est toujours possible de former l’ensemble des tests de Student pour la comparaison des différentes paires de moyennes à l’aide de pwmean comme illsutré ci-dessous :

pwmean bwt, over(race) effects
pwmean bwt, over(race) effects

Pairwise comparisons of means with equal variances

over         : race

------------------------------------------------------------------------------
             |                            Unadjusted           Unadjusted
         bwt |   Contrast   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        race |
      black  |
         vs  |
      white  |  -383.3181   157.8914    -2.43   0.016    -694.8064   -71.82985
      other  |
         vs  |
      white  |  -298.9955   113.6899    -2.63   0.009    -523.2829   -74.70811
      other  |
         vs  |
      black  |   84.32262   165.0131     0.51   0.610    -241.2152    409.8604
------------------------------------------------------------------------------

L’option mcompare() permet d’adapter le type de statistique de test (Tukey, Dunnett), mais dans le cas du modèle de régression précédent il n’y a pas lieu d’appliquer de correction pour les tests multiples ou de modifier la statistique de test. La commande pwmean fournit les mêmes résultats et accepte les mêmes options que pwcompare. La seule différence est qye cette dernière s’utilise en tant que commande de post-estimation et sa syntaxe est plus souple dans le cas des modèles à plusieurs prédicteurs, incluant d’éventuels termes d’interaction.

Traitement de la non linéarité

Il existe plusieurs approches pour traiter le cas d’une relation non linéaire entre la variable réponse et un prédicteur continu.

Voici une illustration avec l’âge de la mère :

scatter bwt age || qfitci bwt age, legend(off)
graph export "fig-03-scatter-bwt-age.pdf", fontface(DroidSans) replace

fig-03-scatter-bwt-age.png

Figure 8 : Relation entre poids des bébés et âge de la mère

L’estimation des paramètres du modèle de régression ne pose pas de difficulté lorsque l’on suppose une simple relation linéaire :

gen agesq = age^2
regress bwt age agesq
gen agesq = age^2
regress bwt age agesq

      Source |       SS           df       MS      Number of obs   =       189
-------------+----------------------------------   F(2, 186)       =      3.89
       Model |   4007561.9         2  2003780.95   Prob > F        =    0.0222
    Residual |  95907736.7       186  515632.993   R-squared       =    0.0401
-------------+----------------------------------   Adj R-squared   =    0.0298
       Total |  99915298.6       188  531464.354   Root MSE        =    718.08

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |  -151.2227    66.3142    -2.28   0.024    -282.0474   -20.39809
       agesq |   3.253674   1.304625     2.49   0.014     .6799088    5.827439
       _cons |   4610.534   817.5263     5.64   0.000     2997.718     6223.35
------------------------------------------------------------------------------

Voici une approche reposant sur des polynômes fractionnaires :

fp <age> : regress bwt <age>
fp <age> : regress bwt <age>
(fitting 44 models)
(....10%....20%....30%....40%....50%....60%....70%....80%....90%....100%)

Fractional polynomial comparisons:
-------------------------------------------------------------------------------
         age |   df    Deviance  Res. s.d.   Dev. dif.   P(*)   Powers
-------------+-----------------------------------------------------------------
     omitted |    0   3027.017    729.016      9.161    0.062               
      linear |    1   3025.497    728.029      7.641    0.059   1           
       m = 1 |    2   3023.143    723.510      5.288    0.076   3           
       m = 2 |    4   3017.856    715.375      0.000       --   3 3         
-------------------------------------------------------------------------------
(*) P = sig. level of model with m = 2 based on F with 184
    denominator dof.

      Source |       SS           df       MS      Number of obs   =       189
-------------+----------------------------------   F(2, 186)       =      4.62
       Model |  4727609.47         2  2363804.73   Prob > F        =    0.0110
    Residual |  95187689.1       186  511761.769   R-squared       =    0.0473
-------------+----------------------------------   Adj R-squared   =    0.0371
       Total |  99915298.6       188  531464.354   Root MSE        =    715.38

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       age_1 |  -.2336151   .1059006    -2.21   0.029    -.4425358   -.0246944
       age_2 |   .0660855   .0287679     2.30   0.023     .0093323    .1228388
       _cons |    3198.04   191.0215    16.74   0.000     2821.192    3574.887
------------------------------------------------------------------------------
fp plot, residuals(none)
graph export "fig-03-fpplot-bwt-age.pdf", fontface(DroidSans) replace

fig-03-fpplot-bwt-age.png

Figure 9 : Utilisation de polynômes fractionnaires pour la relation entre poids des bébés et âge de la mère

Approche robuste

FIXME Find a better illustration + provide more background

Plutôt que de minimiser les écarts quadratiques entre les valeurs prédites et les valeurs observées, il est tout à fait possible d’utiliser un autre type d’estimateur.

Considérons la relation entre le poids des bébés et le poids des mères dont l’ethnicité est black. La commande suivante permet d’afficher un simple diagramme de dispersion ainsi que la droite de régression associée :

twoway (scatter bwt lwt) (lfit bwt lwt) if race == 3
graph export "fig-03-scatter-bwt-lwt-race3.pdf", fontface(DroidSans) replace

fig-03-scatter-bwt-lwt-race3.png

Figure 10 : Relation entre poids des bébés et taille de la mère

Les valeurs ajustées du modèle de régression peuvent être obtenues à l’aide de predict :

regress bwt lwt if race == 3
predict yhols
regress bwt lwt if race == 3

      Source |       SS           df       MS      Number of obs   =        67
-------------+----------------------------------   F(1, 65)        =      3.03
       Model |  1528195.75         1  1528195.75   Prob > F        =    0.0866
    Residual |  32809977.2        65  504768.881   R-squared       =    0.0445
-------------+----------------------------------   Adj R-squared   =    0.0298
       Total |    34338173        66  520275.348   Root MSE        =    710.47

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         lwt |   13.33073   7.661444     1.74   0.087    -1.970233    28.63169
       _cons |   2085.349   422.0535     4.94   0.000      1242.45    2928.249
------------------------------------------------------------------------------
predict yhols
(option xb assumed; fitted values)

La commande robreg disponible dans le package du même nom (ssc install moremata robreg) permet d’estimer les paramètres d’un modèle linéaire en utilisant des M-estimateurs (Huber ou bisquare) [5]. La syntaxe est identique à celle de regress mais il faut faut préciser le type d’estimateur après le nom de la commande : robreg m signifie par exemple une régression avec un estimateur de Huber tandis que robreg s indique à Stata d’utiliser un S-estimateur. Un exemple d’application est disponible dans [6] (FIXME check the reference carefully). Dans le cas présent, on utilisera l’instruction suivante :

quietly robreg m bwt lwt if race == 3
predict yhm

On peut superposer les prédictions de ces deux modèles sur le diagramme de dispersion précédent comme illustré ci-dessous :

twoway (scatter bwt lwt if race == 3) (line yhols yhm lwt, lwidth(*2 *2)), legend(order(2 "OLS" 3 "Huber"))
graph export "fig-03-scatter-bwt-lwt-race3-2.pdf", fontface(DroidSans) replace

fig-03-scatter-bwt-lwt-race3-2.png

Figure 11 : Estimation MCO versus M-estimateur

La régression linéaire multiple

Exemple de base

Diagnostic du modèle

Tests joints et intervalles de confiance simultanés

Spécification de contrastes

Comparaison de modèles emboîtés

Ces des données en cluster

Modèle linéaire et applications

Références

[1] William S. Cleveland. Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association, 74(368):829–836, Dec 1979. [ DOI ]
[2] W. D. Dupont. Statistical Modeling for Biomedical Researchers. Cambridge University Press, 2nd edition, 2009. [ .html ]
[3] Frank E Harrell. Regression Modeling Strategies. Springer Series in Statistics. Springer International Publishing, Cham, 2 edition, 2015.
[4] David W. Hosmer and Stanley Lemeshow. Applied Logistic Regression. Wiley, 2000.
[5] Ben Jann. ROBREG: Stata module providing robust regression estimators. Statistical Software Components, Boston College Department of Economics, January 2010. [ .html ]
[6] E. Vittinghoff, D. V. Glidden, S. C. Shiboski, and C. E. McCulloch. Regression Methods in Biostatistics. Linear, Logistic, Survival, and Repeated Measures Models. Springer, New-York, 2005.