Les paramètres suivants doivent être définis dans R ou RStudio afin de reproduire les analyses présentées dans ce document. Les packages ggplot2, hrbrthemes, directlabels, cowplot, texreg, Hmisc, rms, mfp, multcomp et reshape2 ne font pas partie de la distribution R et doivent être installés si nécessaire :

Les dépendances de ces packages seront installées automatiquement. On supposera également que les instructions R sont exécutées dans un répertoire de travail avec les fichiers de données accessibles dans un sous-répertoire data/. Si ce n’est pas le cas, il suffit de créer un sous-répertoire data/ et d’y enregistrer les fichiers de données, ou de redéfinir les chemins d’accès dans les instructions de lecture des fichiers de données ci-après.

Régression linéaire simple

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 () et \(\beta_1\) la pente () 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, \(\hat\beta_0\) et \(\hat\beta_1\) : \[\begin{array}{l} \hat\beta_0 = \bar y - \hat\beta_1\bar x\\ \hat\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 (\(\hat\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é.

Chargement des données

Les données utilisées proviennent d’une étude sur le volume expiratoire maximum (en litre par seconde) et sont disponibles dans le fichier fev.rda :

## 'data.frame':    654 obs. of  6 variables:
##  $ id    : int  301 451 501 642 901 1701 1752 1753 1901 1951 ...
##  $ age   : atomic  9 8 7 9 9 8 6 6 8 9 ...
##   ..- attr(*, "units")= chr "years"
##  $ fev   : atomic  1.71 1.72 1.72 1.56 1.9 ...
##   ..- attr(*, "units")= chr "liters"
##  $ height: atomic  57 67.5 54.5 53 57 61 58 56 58.5 60 ...
##   ..- attr(*, "units")= chr "inches"
##  $ sex   : Factor w/ 2 levels "female","male": 1 1 1 2 2 1 1 1 1 1 ...
##  $ smoke : Factor w/ 2 levels "non-current smoker",..: 1 1 1 1 1 1 1 1 1 1 ...

Description des variables

Un résumé numérique pour l’ensemble des variables peut être obtenu avec summary (on pourrait utiliser la fonction describe du package Hmisc pour avoir des statistiques plus détaillées) :

##        id             age             fev            height    
##  Min.   :  201   Min.   : 3.00   Min.   :0.791   Min.   :46.0  
##  1st Qu.:15811   1st Qu.: 8.00   1st Qu.:1.981   1st Qu.:57.0  
##  Median :36071   Median :10.00   Median :2.547   Median :61.5  
##  Mean   :37170   Mean   : 9.93   Mean   :2.637   Mean   :61.1  
##  3rd Qu.:53638   3rd Qu.:12.00   3rd Qu.:3.119   3rd Qu.:65.5  
##  Max.   :90001   Max.   :19.00   Max.   :5.793   Max.   :74.0  
##      sex                     smoke    
##  female:318   non-current smoker:589  
##  male  :336   current smoker    : 65  
##                                       
##                                       
##                                       
##

On propose quelques pré-traitements pour faciliter les analyses subséquentes :

## FEV$age [years] 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      654        0       17    0.988    9.931      3.3        5        6 
##      .25      .50      .75      .90      .95 
##        8       10       12       14       15 
##                                                                       
## Value          3     4     5     6     7     8     9    10    11    12
## Frequency      2     9    28    37    54    85    94    81    90    57
## Proportion 0.003 0.014 0.043 0.057 0.083 0.130 0.144 0.124 0.138 0.087
##                                                     
## Value         13    14    15    16    17    18    19
## Frequency     43    25    19    13     8     6     3
## Proportion 0.066 0.038 0.029 0.020 0.012 0.009 0.005

On va s’intéresser à la relation entre la variable fev (ou son log, lfev) et les autres variables numériques, age et height. Voici quelques représentations graphiques préliminaires avec ggplot2, dans un premier temps univariées, puis bivariées :

Modélisation

Un premier modèle de régression que l’on peut considérer est simplement la relation entre FEV et taille sur l’ensemble de l’échantillon. On obtient les paramètres estimés d’un tel modèle à l’aide de lm et en utilisant une “formule” décrivant la relation entre les variables, comme dans le cas de l’analyse de variance :

## 
## Call:
## lm(formula = fev ~ height, data = FEV)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.752 -0.266 -0.004  0.245  2.119 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.43268    0.18146   -29.9   <2e-16
## height       0.05196    0.00116    44.7   <2e-16
## 
## Residual standard error: 0.431 on 652 degrees of freedom
## Multiple R-squared:  0.754,  Adjusted R-squared:  0.753 
## F-statistic: 1.99e+03 on 1 and 652 DF,  p-value: <2e-16

Les graphiques précédents suggérant toutefois une tendance non linéaire lorsque l’on stratifie sur le sexe, on peut regarder si les conditions de validité du modèle (linéarité de la relation entre la variable réponse et le prédicteur et constance de la variance) sont bien vérifiées, en particulier en analysant les résidus de ce modèle.

On pourrait utiliser une transformation de Box-Cox (\(\tfrac{y^{\lambda}-1}{\lambda}\)) ou de Tukey (\(y^{\lambda}\)) pour sélectionner une transformation optimale de la variable réponse pour linéariser la relation, mais on va simplement considérer que la FEV varie linéairement avec le cube de la taille, ce qui revient à modéliser la racine cubique de la FEV, \(\sqrt[3]{\text{fev}}\). En pratique, cela reste assez proche de ce que nous donnerait une transformation de Tukey (on trouve \(\lambda = 0.263 \approx 1/3\) en utilisant les instructions R v <- boxcox(m); v$x[which.max(v$y)]) :

Il reste quelques observations avec des résidus élevés aux deux extrêmes de la distribution des valeurs prédites (les observations avec des résidus supérieurs à 3 en valeur absolue sont annotés avec l’identifiant du sujet). On peut vérifier si ces observations influencent les paramètres estimés à l’aide des mesures d’influence par “jacknife” fournies par R :

## Influence measures of
##   lm(formula = I(fev^(1/3)) ~ height, data = FEV) :
## 
##        dfb.1_  dfb.hght     dffit cov.r   cook.d     hat inf
## 1   -3.38e-02  3.01e-02 -5.12e-02 1.002 1.31e-03 0.00234    
## 2    1.86e-01 -2.03e-01 -2.73e-01 0.943 3.61e-02 0.00343   *
## 3   -8.11e-03  7.54e-03 -9.94e-03 1.007 4.94e-05 0.00361    
## 4   -1.30e-02  1.23e-02 -1.50e-02 1.008 1.13e-04 0.00465    
## 5   -1.42e-02  1.27e-02 -2.15e-02 1.005 2.32e-04 0.00234    
## 6   -2.34e-03  5.00e-04 -1.99e-02 1.004 1.97e-04 0.00153    
## 7   -1.78e-02  1.53e-02 -3.17e-02 1.004 5.03e-04 0.00199    
## 8   -6.90e-02  6.28e-02 -9.37e-02 0.996 4.38e-03 0.00277    
## 9   -1.45e-02  1.21e-02 -2.88e-02 1.004 4.17e-04 0.00186    
## 10  -1.51e-02  1.04e-02 -5.27e-02 0.999 1.38e-03 0.00159    
## 11  -3.69e-03  3.48e-03 -4.25e-03 1.008 9.03e-06 0.00465    
## 12   2.92e-03 -2.73e-03  3.49e-03 1.007 6.11e-06 0.00393    
## 13  -1.26e-03  1.06e-03 -2.51e-03 1.005 3.15e-06 0.00186    
## 14  -7.68e-03  4.23e-03 -3.77e-02 1.002 7.10e-04 0.00155    
##  [ getOption("max.print") est atteint -- 640 lignes omises ]

Autres approches pour la non linéarité

Voici trois approches alternatives permettant de s’affranchir de la stricte linéarité telle qu’assumée dans le modèle de régression linéaire : dans un premier modèle (m1), on peut inclure un terme quadratique pour rendre compte du changement de pente apparent pour des tailles supérieures à 160 cm ; sur la même idée, on peut utiliser des polynômes de degré variable, ici d’ordre 3 (m2) ; enfin, une approche plus souple par rapport aux polynômes consisterait à utiliser des splines cubiques restreints (m3), ici à l’ordre 3 également. Dans le cas des modèles m2 et m3, les fonctions pol et rcs sont fournies par le package rms :

Une dernière approche consiste à utiliser des polynömes fractionnaires, l’avantage étant comme dans le cas des splines une plus grande flexibilité et une validation croisée interne permettant de sélectionner le degré optimal :

## Call:
## mfp(formula = fev ~ fp(height, df = 4, select = 0.05), data = FEV)
## 
## 
## Deviance table:
##           Resid. Dev
## Null model    490.92
## Linear model  120.934
## Final model   111.86
## 
## Fractional polynomials:
##        df.initial select alpha df.final power1 power2
## height          4   0.05  0.05        2      3      .
## 
## 
## Transformations of covariates:
##                  formula
## height I((height/100)^3)
## 
## Rescaled coefficients:
## Intercept   height.1  
## -1.74e-01   7.31e-07  
## 
## Degrees of Freedom: 653 Total (i.e. Null);  652 Residual
## Null Deviance:       491 
## Residual Deviance: 112   AIC: 707

Les prédictions de ces différents modèles sont affichées dans la figure suivante.

Régression, ANOVA et test t

Régression sur variable catégorielle

Avec R, il est important de s’assurer que les variables catégorielles sont bien traitées comme des facteurs, à l’exception des variables binaires codées en 0/1 qui ne posent pas de problème particulier (ni en tant que variable réponse, ni en tant que prédicteur). R utilise des contrastes de traitement par défaut, ce qui signifie que la catégorie de référence à laquelle sont comparés tous les autres niveaux est le premier niveau du facteur.

Voici une application dans laquelle on modélise la FEV en fonction du statut (fumeur/non-fumeur) :

## [1] "non-current smoker" "current smoker"
## 
## Call:
## lm(formula = fev ~ smoke, data = FEV)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.775 -0.634 -0.102  0.480  3.227 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)           2.5661     0.0347   74.04   <2e-16
## smokecurrent smoker   0.7107     0.1099    6.46    2e-10
## 
## Residual standard error: 0.841 on 652 degrees of freedom
## Multiple R-squared:  0.0602, Adjusted R-squared:  0.0588 
## F-statistic: 41.8 on 1 and 652 DF,  p-value: 1.99e-10
##                        2.5 %   97.5 %
## (Intercept)         2.498083 2.634202
## smokecurrent smoker 0.494835 0.926603

On obtiendra le même résultat, au signe près, avec un test t supposant l’égalité des variances :

## 
##  Two Sample t-test
## 
## data:  fev by smoke
## t = -6.464, df = 652, p-value = 1.99e-10
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.926603 -0.494835
## sample estimates:
## mean in group non-current smoker     mean in group current smoker 
##                          2.56614                          3.27686

Considérons maintenant l’âge, discrétiser en 4 classes d’effectifs équilibrés, et comparons les résultats d’une régression linéaire sur variable numérique (age), d’une régression linéaire sur variable catégorielle (age4) et d’une ANOVA :

##      age4     fev       SD
## 1 [ 4, 9) 1.85823 0.420073
## 2 [ 9,11) 2.55157 0.519614
## 3 [11,13) 3.11073 0.639349
## 4 [13,18] 3.59943 0.795803
## fev 
## [1] 0.693344 1.252500 1.741199

Le modèle m2 fournit 3 coefficients de régression pour la variable âge à 4 modalités ; chacun de ces coefficients représente la différence de moyenne de la catégorie en question par rapport à la première catégorie (4–9 ans) :

## 
## Call:
## lm(formula = fev ~ age4, data = FEV)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4167 -0.3882 -0.0632  0.3588  2.1936 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   1.8582     0.0395    47.1   <2e-16
## age4[ 9,11)   0.6933     0.0589    11.8   <2e-16
## age4[11,13)   1.2525     0.0620    20.2   <2e-16
## age4[13,18]   1.7412     0.0665    26.2   <2e-16
## 
## Residual standard error: 0.579 on 650 degrees of freedom
## Multiple R-squared:  0.556,  Adjusted R-squared:  0.554 
## F-statistic:  272 on 3 and 650 DF,  p-value: <2e-16

On retrouvera facilement les différentiels de moyennes via Hmisc :

##      age4     fev       SD
## 1 [ 4, 9) 1.85823 0.420073
## 2 [ 9,11) 2.55157 0.519614
## 3 [11,13) 3.11073 0.639349
## 4 [13,18] 3.59943 0.795803

Voici l’idée schématisée sous forme graphique (les boîtes à moustaches sont alignées sur les centres de classe de la variable age4) :

Régression multiple

Dans le cas où la variable réponse est continue, l’espérance mathématique de la distribution de la variable réponse s’écrit simplement comme une combinaison linéaire des \(p\) termes représentant les effets de chaque régresseur \(x_j\) (\(j=1,\dots,p\)) : \[\mathbb{E}(y\mid x_1,x_2,\dots) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots.\] Sa variance, \(\sigma^2\) peut être estimée par \(s^2 = \tfrac{1}{n-p-1}\sum_{i=1}^n(y_i - \hat y_i)^2\). On peut tester :

  1. l’égalité de l’ensemble des coefficients à l’aide du rapport entre le carré moyen de la régression et celui de la résiduelle (test F à \(p\) et \(n-p-1\) degrés de liberté), et l’hypothèse nulle se lit \(H_0:\, \beta_1 = \beta_2 = \dots = \beta_p = 0\) ;
  2. la nullité d’un coefficient à l’aide du rapport \(\tfrac{\hat\beta_j}{\text{SE}(\hat\beta_j)}\) (test t) ;
  3. la nullité d’un ensemble de coefficients (test F partiel par comparaison de modèles emboîtés à \(p\) et \(p-m\) degrés de libertés).

Les coefficients de régression représentent l’effet partiel de chaque prédicteur, c’est-à-dire l’effet de l’augmentation d’une unité de \(x_j\) sur \(y\) lorsque tous les autres prédicteurs sont maintenus constants.

Un autre modèle possible consiste à considérer le statut fumeur et un potentiel facteur de confusion, à savoir l’âge puisqu’en dessous de 9 ans il n’y a pas de fumeurs et qu’au-dessus on a bien deux sous-populations, comme on peut le vérifier dans la figure suivante.

Voici le modèle en question :

## 
## ===========================================
##                      Model 1     Model 2   
## -------------------------------------------
## (Intercept)            2.57 ***    0.35 ***
##                       (0.03)      (0.08)   
## smokecurrent smoker    0.71 ***   -0.21 ** 
##                       (0.11)      (0.08)   
## age                                0.23 ***
##                                   (0.01)   
## -------------------------------------------
## R^2                    0.06        0.58    
## Adj. R^2               0.06        0.58    
## Num. obs.            654         654       
## RMSE                   0.84        0.56    
## ===========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

Comme les modèles sont emboîtés, le test du coefficient pour la variable age peut se retrouver par simple comparaison de modèles (par test F ou par test du rapport de vraisemblance) :

## Analysis of Variance Table
## 
## Model 1: fev ~ smoke
## Model 2: fev ~ age + smoke
##   Res.Df   RSS Df Sum of Sq     F Pr(>F)
## 1    652 461.4                          
## 2    651 207.1  1     254.3 799.3 <2e-16

Et voici le modèle incluant un terme d’interaction entre les variables age et smoke :

## 
## Call:
## lm(formula = fev ~ age + smoke + age:smoke, data = FEV)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7757 -0.3471 -0.0327  0.3347  2.0575 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)              0.24206    0.08302    2.92   0.0037
## age                      0.24370    0.00837   29.11  < 2e-16
## smokecurrent smoker      1.90489    0.42493    4.48  8.7e-06
## age:smokecurrent smoker -0.15996    0.03159   -5.06  5.4e-07
## 
## Residual standard error: 0.554 on 650 degrees of freedom
## Multiple R-squared:  0.594,  Adjusted R-squared:  0.592 
## F-statistic:  317 on 3 and 650 DF,  p-value: <2e-16

Le graphique suivant résume la situation, avec des intervalles de confiance à 95 % (non simultanés !) :

Notons que l’on peut toujours utiliser la fonction ols au lieu de la fonction lm de base. Les sorties sont plus riches, il est possible d’introduire une régularisation de type “ridge” et rms fournit des procédure de validation et de calibration du modèle. les commandes de post-estimation (en langage Stata) sont plus intéressantes : on pourrait par exemple calculer les valeurs prédites de fev en fonction de smoke conditionnellement à un ensemble fixe et pré-spécifié de valeurs de age.

En dernier lieu, si l’on rajoute la variable height et height2, on remarquera que l’interaction age:smoke n’est plus significative :

## Linear Regression Model
##  
##  ols(formula = fev ~ age + smoke + age:smoke + pol(height, 2), 
##      data = FEV, x = TRUE)
##  
##                  Model Likelihood     Discrimination    
##                     Ratio Test           Indexes        
##  Obs     654    LR chi2    1025.54    R2       0.792    
##  sigma0.3974    d.f.             5    R2 adj   0.790    
##  d.f.    648    Pr(> chi2)  0.0000    g        0.878    
##  
##  Residuals
##  
##       Min       1Q   Median       3Q      Max 
##  -1.70600 -0.24141  0.00623  0.23525  1.81959 
##  
##  
##                             Coef    S.E.   t     Pr(>|t|)
##  Intercept                   7.9858 1.4768  5.41 <0.0001 
##  age                         0.0707 0.0099  7.11 <0.0001 
##  smoke=current smoker        0.1479 0.3166  0.47 0.6407  
##  height                     -0.1224 0.0192 -6.36 <0.0001 
##  height^2                    0.0005 0.0001  8.62 <0.0001 
##  age * smoke=current smoker -0.0226 0.0237 -0.96 0.3395  
## 
##   age              smoke height linear.predictors    se.fit
## 1   9 non-current smoker    120           1.61428 0.0865819
## 2  10 non-current smoker    120           1.68496 0.0927533
## 3  11 non-current smoker    120           1.75564 0.0995365
## 4  12 non-current smoker    120           1.82632 0.1068151
## 5  13 non-current smoker    120           1.89700 0.1144946
## 6  14 non-current smoker    120           1.96768 0.1224997

Approche alternative (spécifique au package rms) :

##    age height          smoke    yhat   lower   upper
## 1   18    160 current smoker 3.06612 2.82602 3.30621
## 2   18    162 current smoker 3.16466 2.92621 3.40310
## 3   18    164 current smoker 3.26746 3.03050 3.50442
## 4   18    166 current smoker 3.37453 3.13881 3.61024
## 5   18    168 current smoker 3.48586 3.25108 3.72064
## 6   18    170 current smoker 3.60145 3.36721 3.83570
## 7   18    172 current smoker 3.72131 3.48711 3.95552
## 8   18    174 current smoker 3.84544 3.61065 4.08023
## 9   18    176 current smoker 3.97383 3.73772 4.20993
## 10  18    178 current smoker 4.10648 3.86819 4.34477
## 11  18    180 current smoker 4.24339 4.00192 4.48486
## 
## Response variable (y): fev 
## 
## Limits are 0.95 confidence limits

Analyse de covariance

Soit \(y_{ij}\) la \(j\)ème observation dans le groupe \(i\). À l’image du modèle d’ANOVA à un facteur, le modèle d’ANCOVA s’écrit \[ y_{ij} = \mu+\alpha_i+\beta(x_{ij}-\bar x)+\varepsilon_{ij},\] où \(\beta\) est le coefficient de régression liant la réponse \(y\) et le cofacteur \(x\) (continu), avec \(\bar x\) la moyenne générale des \(x_{ij}\), et toujours un terme d’erreur \(\varepsilon_{ij}\sim \mathcal{N}(0,\sigma^2)\).

Notons que l’on fait l’hypothèse que \(\beta\) est le même dans chaque groupe. Cette hypothèse de parallélisme peut se vérifier en testant la significativité du terme d’interaction \(\alpha\beta\). La réponse moyenne ajustée pour l’effet du co-facteur numérique s’obtient simplement comme \(\bar\alpha_i+\hat\beta(\bar x_i-\bar x)\), où \(\bar x_i\) est la moyenne des \(x\) dans le \(i\)ème groupe.

Chargement des données

Les données anorexia portent sur une étude clinique dans laquelle on s’intéresse à l’amélioration de l’état de santé (mesurée par le gain de poids) de patientes anorexiques prises en charge par thérapie familiale (FT) ou thérapie comportementale (CBT). Un groupe contrôle (Cont) est également disponible. Les données sont disponibles dans le package MASS.

## 'data.frame':    72 obs. of  3 variables:
##  $ Treat : Factor w/ 3 levels "CBT","Cont","FT": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Prewt : num  80.7 89.4 91.8 74 78.1 88.3 87.3 75.1 80.6 78.4 ...
##  $ Postwt: num  80.2 80.1 86.4 86.3 76.1 78.1 75.1 86.7 73.5 84.6 ...

Description des variables

Voici un bref résumé numérique, suivi d’une description graphique des données :

## anorexia 
## 
##  3  Variables      72  Observations
## ---------------------------------------------------------------------------
## Treat 
##        n  missing distinct 
##       72        0        3 
##                             
## Value       Cont   CBT    FT
## Frequency     26    29    17
## Proportion 0.361 0.403 0.236
## ---------------------------------------------------------------------------
## Prewt 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##       72        0       58        1    37.08    2.645    33.18    34.43 
##      .25      .50      .75      .90      .95 
##    35.82    37.03    38.70    39.91    40.33 
## 
## lowest : 31.500 31.725 32.535 33.030 33.300, highest: 40.230 40.455 41.310 42.390 42.705
## ---------------------------------------------------------------------------
## Postwt 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##       72        0       67        1    38.33    4.148    33.05    33.94 
##      .25      .50      .75      .90      .95 
##    35.70    37.82    41.20    43.26    44.66 
## 
## lowest : 32.085 32.625 32.850 33.030 33.075, highest: 44.280 45.135 45.180 45.720 46.620
## ---------------------------------------------------------------------------

Modélisation

Dans un premier temps, après avoir exclus le groupe contrôle qui pourrait faire l’objet d’analyses séparées avec un ensemble de contrastes appropriés, on souhaite vérifier s’il existe ou non une interaction entre le type de thérapie et le poids de sortie, après avoir contrôlé pour le poids en entrée :

## Analysis of Variance Table
## 
## Model 1: Postwt ~ Prewt + Treat
## Model 2: Postwt ~ Prewt + Treat + Prewt:Treat
##   Res.Df   RSS Df Sum of Sq     F Pr(>F)
## 1     43 465.3                          
## 2     42 465.1  1    0.1897 0.017  0.896

En l’absence d’interaction, on pourra donc estimer le gain de poids moyen entre les deux groupes, quel que soit le poids en entrée :

## 
## Call:
## lm(formula = Postwt ~ Prewt + Treat, data = anorexia, subset = Treat != 
##     "Cont")
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.812 -1.622 -0.503  1.779  7.317 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    6.144      8.377    0.73  0.46724
## Prewt          0.871      0.225    3.88  0.00035
## TreatFT        1.947      1.006    1.94  0.05957
## 
## Residual standard error: 3.29 on 43 degrees of freedom
## Multiple R-squared:  0.314,  Adjusted R-squared:  0.282 
## F-statistic: 9.84 on 2 and 43 DF,  p-value: 0.000303
##                   2.5 %   97.5 %
## (Intercept) -10.7487231 23.03692
## Prewt         0.4184723  1.32402
## TreatFT      -0.0820634  3.97664

Exercices d’application

Dans cette application, on travaillera avec le data frame birthwt qui contient des données issues d’une étude rétrospective dans lequel on s’intéresse aux facteurs de risque d’accoucher d’un bébé avec un poids inférieur à la norme (< 2.5 kg selon les normes américaines des années 80). Les variables d’intérêt sont le poids du bébé (bwt), l’âge (age) et le poids (lwt, en pounds) de la mère ainsi que son origine ethnique (race), les antécédents d’hypertension (ht) et le status fumeur (smoke).

Les données sont disponibles dans le package MASS. Il est nécessaire de charger au préalable le package ou de le spécifier lors de l’importation :

##       low             age            lwt           race     
##  Min.   :0.000   Min.   :14.0   Min.   : 80   Min.   :1.00  
##  1st Qu.:0.000   1st Qu.:19.0   1st Qu.:110   1st Qu.:1.00  
##  Median :0.000   Median :23.0   Median :121   Median :1.00  
##  Mean   :0.312   Mean   :23.2   Mean   :130   Mean   :1.85  
##  3rd Qu.:1.000   3rd Qu.:26.0   3rd Qu.:140   3rd Qu.:3.00  
##  Max.   :1.000   Max.   :45.0   Max.   :250   Max.   :3.00  
##      smoke            ptl              ht               ui       
##  Min.   :0.000   Min.   :0.000   Min.   :0.0000   Min.   :0.000  
##  1st Qu.:0.000   1st Qu.:0.000   1st Qu.:0.0000   1st Qu.:0.000  
##  Median :0.000   Median :0.000   Median :0.0000   Median :0.000  
##  Mean   :0.392   Mean   :0.196   Mean   :0.0635   Mean   :0.148  
##  3rd Qu.:1.000   3rd Qu.:0.000   3rd Qu.:0.0000   3rd Qu.:0.000  
##  Max.   :1.000   Max.   :3.000   Max.   :1.0000   Max.   :1.000  
##       ftv             bwt      
##  Min.   :0.000   Min.   : 709  
##  1st Qu.:0.000   1st Qu.:2414  
##  Median :0.000   Median :2977  
##  Mean   :0.794   Mean   :2945  
##  3rd Qu.:1.000   3rd Qu.:3487  
##  Max.   :6.000   Max.   :4990

Voici une vue partielle des données individuelles et agrégées :

Sélection de modèle

On considérera le modèle de base suivant : bwt ~ lwt + race.

  1. Quelles sont les variables que l’on pourrait ajouter pour améliorer la portée ou la pertinence d’un tel modèle sur le plan de l’interprétation clinique ? Ces variables sont-elles significatives ?
  2. En utilisant une procédure de sélection automatique (step), indiquer parmi l’ensemble des variables lesquelles pourraient être considérées comme significatives au seuil de 5 %, sans tenir compte des problèmes inhérents à la sélection pas à pas.
  3. Comparer avec les résultats obtenus par une approche reposant sur une pénalisation de type “elasticnet” (package elasticnet ou glmnet).

Tests spécifiques

  1. Considérons le modèle bwt ~ lwt + race + smoke. Estimer les paramètres du modèle avec leurs intervalles de confiance à 95 %.
  2. On souhiate comparer le groupe white aux deux groupes réunis black+other tout en ajustant sur les deux autres variables du modèle. Proposer une approche reposant sur l’utilisation de contrastes (package multcomp ou rms).
  3. Existe t-il une différence signifciative entre les mères qui fument dans la catégorie white et celles qui fument dans la catégorie black ?

Estimation par bootstrap

  1. Estimer les paramètres du modèle bwt ~ age + ht et construire un intervalle de confiance à 95 % pour le coefficient de régression de ht.
  2. Comparer avec une approche par bootstrap (package boot ou procédure manuelle) pour la construction d’un intervalle de confiance.