Stata : modélisation statistique (1)

Table des matières

Dans ce chapitre, 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 linéaire 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. [8] et Dupont [3] couvrent l’essentiel des notions présentées. Pour aller plus loin, le livre de Harrell [4] demeure une référence en ce qui concerne les techniques biostatistiques modernes appliquées dans un contexte médical. Bien que reposant sur le logiciel R, la plupart des techniques discutées dans cet ouvrage sont disponibles sous Stata.

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 du « meilleur » 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

Un exemple 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.

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 [5]. 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

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.eps", 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 [2] sur le diagramme de dispersion précédent afin d’évaluer visuellement les écarts à la linéarité concernant la relation entre ces deux variables. Une courbe « loess », ou lowess (locally weighted scatterplot smoothing), fait partie de la famille des régressions polynomiales « locales », c’est-à-dire des modèles non-paramétriques de régression tenant compte du voisinage des observations, un peu à l’image des techniques telles que les \(k\) plus proches voisins dans une tâche de classification (car l’estimation est pondérée) ou les moyennes mobiles en séries chronologiques (car on définit une fenêtre de lissage que l’on déplace sur l’axe des \(x\)). Le poids de chaque observation est déterminé par sa distance par rapport au point sur lequel se centre la régression, et généralement la fonction de pondération utilisée est de type tri-cubique comme proposée par [2], \((1 - (\text{dist}/\text{maxdist})^3)^3)\), ce qui revient à donner un poids maximal à l’observation centrale, les observations les plus distantes plus de cette dernière contribuant le moins à la régression polynomiale. Stata autorise n’importe quelle combinaison des options mean (utiliser la moyenne des observations, comme dans une moyenne mobile, au lieu des valeurs prédites par la régression) et noweight (l’utilisation d’une fonction de pondération tri-cubique ou non). À noter qu’il s’agit d’une approche quelque peu intensive sur le plan des calculs : pour \(n\) observations, \(n\) régressions sont réalisées.

La commande lowess peut être combinée à scatter, en utilisant la commande twoway et en regroupant les différentes instructions graphiques à l’aide de l’opérateur || ou de paires de parenthèses :

graph twoway scatter bwt lwt || lowess bwt lwt, legend(off)
graph export "fig-03-lowess-bwt-lwt.eps", 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 (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 entre eux ainsi que de \(x\)). 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 estimé à l’aide de \(\widehat\beta_0\) et que l’on notera \(b_0\). Ici, \(b_0=2373.5\). Cette valeur reste peu interprétable puisqu’elle représente le poids attendu pour un bébé lorsque le poids de la mère est de 0 kg. La pente (\(b1=9.8\)) indique de combien varie bwt lorsque lwt varie d’une unité, c’est-à-dire d’un kilogramme. 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

L’instruction _b[lwt] est une variable dite variable « système » stockées en mémoire et mise à jour après cahque commande d’estimation par Stata. Les variables _n et _rc sont d’autres exemples de telles variables système ((U) 13.4). Il est toutefois possible de sauvegarder ces résultats d’estimations à l’aide de la commande estimates ou en stockant la matrice virtuelle e(b) dans une macro locale, mais dans ce cas on ne plus indexer les valeurs par le nom des variables :

matrix b = e(b)
display b[1,1]
matrix b = e(b)
display b[1,1]
9.7837929

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.eps", fontface(DroidSans) replace

fig-03-qnorm-r.png

Figure 3 : Distribution des résidus simples

Un histogramme ou une courbe de densité permet d’examiner rapidement la forme de la distribution des résidus. Voici un exemple de courbe de densité construite avec kdensity, à laquelle on ajoute une courbe de densité normale à l’aide de l’option normal :

kdensity rs, normal normopts(lpat(--))
graph export "fig-03-kdensity-rs.eps", 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, qui fournira la même information dans le cas d’une régression avec un seul prédicteur) :

graph twoway scatter rs yhat, yline(0)
graph export "fig-03-scatter-rs-yhat.eps", 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.eps", 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 à l’aide du préfixe i., on obtiendra strictement le même résultat du fait du codage initial en 0/1 où une variation d’une unité correspond au passage de la classe « non fumeur » à la classe « fumeur » :

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 à présent 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. L’indicatruice exclue servira de catégorie de référence. 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 de référence. En utilisant le préfixe ib3, par exemple, 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.eps", 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 illustré 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 que 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.

Voici une autre illustration, cette fois-ci avec les données d’un essai clinique randomisé visant à évaluer l’effet de l’administration d’ibuprofène par voie intraveineuse sur la mortalité de patients en état septique sévère [1]. Les données, disponibles dans le fichier sepsis.dta, sont largement exploitées dans l’ouvrage de William Dupont [3]. Au total, le tableau de données est composé 22 variables dont 16 variables représentant une mesure de la température entre \(T_0\) et \(T_0 + 15 \times 2\) h, deux groupes de patients (« Placebo », n = 231 et « Ibuprofène », n = 224) et une mesure de morbidité (score APACHE).

use "data/sepsis.dta", replace
describe, simple
table treat, content(mean temp0 mean temp1 mean temp6) format(%5.1f)
use "data/sepsis.dta", replace
describe, simple
id        apache    followup  temp2     temp5     temp8     temp11    temp14
treat     o2del     temp0     temp3     temp6     temp9     temp12    temp15
race      fate      temp1     temp4     temp7     temp10    temp13
table treat, content(mean temp0 mean temp1 mean temp6) format(%5.1f)

-------------------------------------------------
Treatment | mean(temp0)  mean(temp1)  mean(temp6)
----------+--------------------------------------
  Placebo |       100.5        100.2         99.7
Ibuprofen |       100.4         99.5         98.4
-------------------------------------------------

Voici comment générer un aperçu des données individuelles longitudinale, en se limitant à la période 0-6 heures :

keep id treat temp0-temp6
reshape long temp, i(id) j(hour)
replace temp = (temp-32) / 1.8
graph twoway (scatter temp hour, ms(none) lcol(gs15) connect(l)) (scatter temp hour if hour < 2, ms(none) connect(l)), by(treat, legend(off)) xtitle(Time unit (x2 hours)) ytitle (Temperature (°C))
graph export "fig-03-scatter-temp-hour.eps", fontface(DroidSans) replace

fig-03-scatter-temp-hour.png

Figure 8 : Évolution de la température après la prise en charge dans les deux groupes de patients

Bien que la technique appropriée pour modéliser l’évolution de la température entre \(T_0\) et \(T_1\) entre les deux groupes soit une analyse de covariance, voici en attendant les questions auxquelles il est possible de répondre à l’aide de simples tests de Student. Premièrement, les deux groupes sont-ils comparables à \(T_0\) (\(H_0\) : temp0(ibuprofène) = temp0(placebo)) ? Voici l’instruction Stata correspondante :

ttest temp if hour == 0, by(treat)
ttest temp if hour == 0, by(treat)

Two-sample t test with equal variances
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
 Placebo |     231    38.05012    .0699302    1.062847    37.91233    38.18791
Ibuprofe |     224    37.97862    .0793881    1.188173    37.82217    38.13507
---------+--------------------------------------------------------------------
combined |     455    38.01492    .0527696    1.125614    37.91122    38.11862
---------+--------------------------------------------------------------------
    diff |            .0714991    .1056147               -.1360564    .2790546
------------------------------------------------------------------------------
    diff = mean(Placebo) - mean(Ibuprofe)                         t =   0.6770
Ho: diff = 0                                     degrees of freedom =      453

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.7506         Pr(|T| > |t|) = 0.4988          Pr(T > t) = 0.2494

Deuxièmement, les deux groupes sont-ils comparables à \(T_1\) en terme d’évolution \((T_0-T_1)\) (\(H_0\) : temp0−temp1(ibuprofène) = temp0−temp1(placebo)) ?

quietly reshape wide
gen difftemp = temp0 - temp1
ttest difftemp, by(treat)
quietly reshape wide
gen difftemp = temp0 - temp1
(35 missing values generated)
ttest difftemp, by(treat)

Two-sample t test with equal variances
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
 Placebo |     212    .1733228    .0408455     .594719    .0928054    .2538403
Ibuprofe |     208    .4913999    .0478614    .6902662    .3970417    .5857581
---------+--------------------------------------------------------------------
combined |     420    .3308467    .0323248    .6624603    .2673078    .3943856
---------+--------------------------------------------------------------------
    diff |            -.318077    .0628323               -.4415837   -.1945704
------------------------------------------------------------------------------
    diff = mean(Placebo) - mean(Ibuprofe)                         t =  -5.0623
Ho: diff = 0                                     degrees of freedom =      418

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.0000         Pr(|T| > |t|) = 0.0000          Pr(T > t) = 1.0000

Enfin, troisièment, on pourrait se demander s’il y a une évolution significative entre \(T_0\) et \(T_1\) pour le groupe traité : il s’agit cette fois d’un test t pour données appariées. Voici le code correspondant :

ttest temp0 == temp1 if treat == 1
ttest temp0 == temp1 if treat == 1

Paired t test
------------------------------------------------------------------------------
Variable |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
   temp0 |     208     38.0031    .0829791    1.196742    37.83951    38.16669
   temp1 |     208     37.5117    .0714197    1.030029     37.3709     37.6525
---------+--------------------------------------------------------------------
    diff |     208    .4913999    .0478614    .6902662    .3970417    .5857581
------------------------------------------------------------------------------
     mean(diff) = mean(temp0 - temp1)                             t =  10.2672
 Ho: mean(diff) = 0                              degrees of freedom =      207

 Ha: mean(diff) < 0           Ha: mean(diff) != 0           Ha: mean(diff) > 0
 Pr(T < t) = 1.0000         Pr(|T| > |t|) = 0.0000          Pr(T > t) = 0.0000

Par une approche de régression simple, on obtiendrait essentiellement des réponses similaires. Voici déjà une commande permettant d’estimer les paramètres du modèle dans les deux groupes :

quietly reshape long
bysort treat: regress temp hour, noheader
quietly reshape long
bysort treat: regress temp hour, noheader

-------------------------------------------------------------------------------
-> treat = Placebo
------------------------------------------------------------------------------
        temp |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        hour |   -.072687   .0138938    -5.23   0.000    -.0999401   -.0454339
       _cons |   38.00192    .050157   757.66   0.000     37.90353     38.1003
------------------------------------------------------------------------------

-------------------------------------------------------------------------------
-> treat = Ibuprofen
------------------------------------------------------------------------------
        temp |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        hour |  -.1675952   .0132796   -12.62   0.000    -.1936444   -.1415459
       _cons |   37.72096    .047841   788.46   0.000     37.62712    37.81481
------------------------------------------------------------------------------

À partir de là, on souhaite comparer les coefficients de régression entre les deux groupes. Pour cela, il y a deux approches possibles. D’un côté il est possible de reconnaître qu’il s’agit essentiellement d’un test de l’interaction entre les deux variables hour et treat, et c’est sans doute l’approche la plus simple de la question. Dans ce cas, il suffit de générer le terme d’interaction et de tester ses composantes directement :

quietly tabulate treat, gen(treat)
generate treat1hour = treat1*hour
generate treat2hour = treat2*hour
quietly regress temp treat1 treat2 treat1hour treat2hour
test treat1hour treat2hour
quietly tabulate treat, gen(treat)
generate treat1hour = treat1*hour
generate treat2hour = treat2*hour
quietly regress temp treat1 treat2 treat1hour treat2hour
test treat1hour treat2hour

 ( 1)  treat1hour = 0
 ( 2)  treat2hour = 0

       F(  2,  2966) =   89.02
            Prob > F =    0.0000

Le test ci-dessus est un test simultané (2 degrés de liberté) pour la nullité des termes d’interaction, tandis que le test ci-dessous permet d’évaluer l’égalité de ces deux termes :

test treat1hour = treat2hour
test treat1hour = treat2hour

 ( 1)  treat1hour - treat2hour = 0

       F(  1,  2966) =   24.33
            Prob > F =    0.0000

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. Revenons aux données sur les poids de naissance pour illustrer avec l’âge de la mère quelques-unes des approches possibles :

clear all
webuse lbw
scatter bwt age || qfitci bwt age, legend(off)
graph export "fig-03-scatter-bwt-age.eps", fontface(DroidSans) replace

fig-03-scatter-bwt-age.png

Figure 9 : 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 incluant l’âge et le carré de l’âge :

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
------------------------------------------------------------------------------

Le terme quadratique améliore t-il la qualité d’ajustement d’un tel modèle ? Ici, on voit que le \(R^2\) ajusté est de 3 %, ce qui ne change pas vraiment des résultats observés dans le cas d’une régression simple. On peut le vérifier également au niveau des indices AIC ou BIC :

quietly regress bwt age agesq
estimates store m1
quietly regress bwt age
estimates store m0
estimates stats m*
quietly regress bwt age agesq
estimates store m1
quietly regress bwt age
estimates store m0
estimates stats m*

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |    Obs    ll(null)   ll(model)     df          AIC         BIC
-------------+---------------------------------------------------------------
          m1 |    189   -1513.509    -1509.64      3      3025.28    3035.005
          m0 |    189   -1513.509   -1512.748      2     3029.497     3035.98
-----------------------------------------------------------------------------
               Note:  N=Obs used in calculating BIC; see [R] BIC note

Une autre approche repose sur l’utilisation de polynômes fractionnaires, qui ont été largement développés et popularisés par Royston et coll. [7]. L’idée générale est de considérer des polynômes dont les exposants sont pris dans un ensemble prédéfini de valeurs \(P = {-2, -1, -0.5, 0, 0.5, 1, 2, 3}\), où par convention \(x^{(0)} = \ln(x)\). Un polynôme fractionnaire de degré \(m\) se construit comme \(\text{FPm} = \beta_0 + \sum_{j=1}^m \beta_jx^{(p_j)}\), où \(p_j \in P\). On notera que pour un polynôme de degré \(m\), une même puissance peut être répétée \(m\) fois.

Stata 13 dispose de la commande fracpoly mais il est recommendé d’utiliser les commandes fp (cas univarié) et mfp (cas multivarié) qui permettent de construire automatiquement les termes d’un ou plusieurs polynômes fractionnaires pour une variable numérique donnée. Voici un exemple d’application sur la variable age :

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.eps", fontface(DroidSans) replace

fig-03-fpplot-bwt-age.png

Figure 10 : 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.eps", fontface(DroidSans) replace

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

Figure 11 : 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.11
       Model |  1566727.54     1  1566727.54           Prob > F      =  0.0826
    Residual |  32771445.4    65  504176.084           R-squared     =  0.0456
-------------+------------------------------           Adj R-squared =  0.0309
       Total |    34338173    66  520275.348           Root MSE      =  710.05

------------------------------------------------------------------------------
         bwt |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         lwt |   6.133087   3.479153     1.76   0.083     -.815261    13.08143
       _cons |   2067.861   426.5168     4.85   0.000     1216.048    2919.674
------------------------------------------------------------------------------
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) [6]. Dans le cas d’une approche par M-estimation, les estimés des paramètres du modèle de régression sont obtenus en minimisant une fonction de coût, \(\rho\), reposant sur la valeur des résidus sur l’ensemble des valeurs de \(X\). Spécifiquement, on recherche une fonction \(\rho(e) \ge 0\), symétrique et monotone,

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 [8] (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.eps", fontface(DroidSans) replace

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

Figure 12 : Estimation MCO versus M-estimateur

La régression linéaire multiple

#TODO

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

#TODO

Références

[1] G. R. Bernard, A. P. Wheeler, J. A. Russell, R. Schein, W. R. Summer, K. P. Steinberg, W. J. Fulkerson, P. E. Wright, B. W. Christman, W. D. Dupont, S. B. Higgins, and B. B. Swindell. The effects of ibuprofen on the physiology and survival of patients with sepsis. the ibuprofen in sepsis study group. New England Journal of Medicine, 336(13):912--918, 1997.
[2] William S. Cleveland. Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association, 74(368):829–836, Dec 1979. [ DOI ]
Keywords: dataviz
[3] W. D. Dupont. Statistical Modeling for Biomedical Researchers. Cambridge University Press, 2nd edition, 2009. [ .html ]
Keywords: Stata
[4] Frank E Harrell. Regression Modeling Strategies. Springer Series in Statistics. Springer International Publishing, Cham, 2 edition, 2015.
[5] David W. Hosmer and Stanley Lemeshow. Applied Logistic Regression. Wiley, 2000.
[6] Ben Jann. ROBREG: Stata module providing robust regression estimators. Statistical Software Components, Boston College Department of Economics, January 2010. [ .html ]
[7] P. Royston and D. G. Altman. Regression using fractional polynomials of continuous covariates: Parsimonious parametric modelling. Applied Statistics, 43(3):429--467, 1994.
[8] 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.
Keywords: Stata