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. [1] et Dupont [2] couvrent l'essentiel des notions présentées. Pour aller plus loin, le livre de [3] 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 [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
Ensuite, on utilise la commande scatter pour construire le diagramme de dispersion :
graph twoway scatter bwt lwt graph export "fig-03-scatter-bwt-lwt.eps", replace
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 [5] 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 [5], \((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", replace
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", replace
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", replace
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", replace
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", replace
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", replace
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 [bernard-1997-effec-ibupr]. Les données, disponibles dans le fichier .sepsis.dta, sont largement exploitées dans l'ouvrage de William Dupont [2]. 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", replace
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", replace
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. [6]. 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", replace
Figure 10 : Utilisation de polynômes fractionnaires pour la relation entre poids des bébés et âge de la mère
Approche robuste
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", replace
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) [7]. 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 [1] (
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", replace
Figure 12 : 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] | E. Vittinghoff, D. V. Glidden, S. C. Shiboski, and C. E. McCulloch, Regression Methods in Biostatistics. Linear, Logistic, Survival, and Repeated Measures Models. New-York: Springer, 2005. |
| [2] | W. D. Dupont, Statistical Modeling for Biomedical Researchers. Cambridge University Press, 2nd ed., 2009. |
| [3] | F. E. Harrell, Regression Modeling Strategies. Springer Series in Statistics, Cham: Springer International Publishing, 2 ed., 2015. |
| [4] | D. W. Hosmer and S. Lemeshow, Applied Logistic Regression. Wiley, 2000. |
| [5] | W. S. Cleveland, “Robust locally weighted regression and smoothing scatterplots,” Journal of the American Statistical Association, vol. 74, p. 829–836, Dec. 1979. |
| [6] | P. Royston and D. G. Altman, “Regression using fractional polynomials of continuous covariates: Parsimonious parametric modelling,” Applied Statistics, vol. 43, no. 3, pp. 429--467, 1994. |
| [7] | B. Jann, “ROBREG: Stata module providing robust regression estimators.” Statistical Software Components, Boston College Department of Economics, Jan. 2010. [ .html ] |