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 et rms ne font pas partie de la distribution R et doivent être téléchargés au préalable si nécessaire :

install.packages(c("ggplot2", "hrbrthemes", "directlabels", "cowplot", "texreg", "rms")

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.

library(hrbrthemes)
library(directlabels)
library(ggfortify)
library(cowplot)
library(texreg)
library(survival)
library(rms)            ## Imports: Hmisc, ggplot2
options(digits = 6, show.signif.stars = FALSE)
theme_set(theme_ipsum(base_size = 11))

Régression logistique

Si l’on note $\pi$ la probabilité d’observer l’événement $y=1$ (vs. 0), alors le log odds (transformation ) peut s’exprimer comme une fonction linéaire des paramètres du modèle à $p$ prédicteurs : \[g(x)=\log\left(\frac{\pi}{1-\pi}\right)=\beta_0+\beta_1x_1+\dots+\beta_px_p,\] et la probabilité prédite s’écrit alors \[P(y=1\mid x_1,x_2,\dots,x_p)=\hat y_i=\frac{\exp(\hat\beta_0+\hat\beta_1x_1+\dots+\hat\beta_px_p)}{1+\exp(\hat\beta_0+\hat\beta_1x_1+\dots+\hat\beta_px_p)}.\]

Dans ce type de modèle, on fait l’hypothèse que $y_i$ suit une distribution binomiale et que les observations sont indépendantes (aucune hypothèse sur la variance qui n’est pas un paramètre dans ce cas de figure). Notons également l’absence de terme d’erreur. L’estimation d’un tel type de modèle se fait par la méthode du maximum de vraisemblance. D’autres fonctions de lien existent (probit, log-log).

Chargement des données

Les données birthwt ont déjà été examinées sous l’angle de la régression linéaire en considérant la variable bwt (poids du bébé) comme variable à expliquer. Dans les analyses qui suivent, on va s’intéresser à la variable low qui représente une indicatrice pour le cas bwt<2500.

data(birthwt, package = "MASS")
str(birthwt)

## 'data.frame':    189 obs. of  10 variables:
##  $ low  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : int  2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: int  0 0 1 1 1 0 0 0 1 1 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ui   : int  1 0 0 1 1 0 0 0 0 0 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...

Voici quelques pré-traitements sur les variables :

birthwt$lwt <- birthwt$lwt * 0.45
birthwt$race <- factor(birthwt$race, levels = 1:3, labels = c("white", "black", "other"))
birthwt$smoke <- factor(birthwt$smoke, levels = 0:1, labels = c("non-smoker", "smoker"))
birthwt$ftv[birthwt$ftv > 2] <- 2

Voici un résumé de la structure de données en considérant la variable bwt et en surlignant les observations pour lesquelles low=1 :

Description des variables

Voici un tableau récapitualitif des variables d’intérêt construit avec Hmisc :

s <- summary(low ~ age + lwt + race + smoke + ht + ui + ftv,
             data = birthwt, method = "reverse")
print(s, exclude1 = FALSE, npct = "both")

## 
## 
## Descriptive Statistics by low
## 
## +------------------+--------------------+--------------------+
## |                  |0                   |1                   |
## |                  |(N=130)             |(N=59)              |
## +------------------+--------------------+--------------------+
## |age               |   19.0/23.0/28.0   |   19.5/22.0/25.0   |
## +------------------+--------------------+--------------------+
## |lwt               |50.850/55.575/66.150|46.800/54.000/58.500|
## +------------------+--------------------+--------------------+
## |race : white      |     56% 73/130     |     39% 23/ 59     |
## +------------------+--------------------+--------------------+
## |    black         |     12% 15/130     |     19% 11/ 59     |
## +------------------+--------------------+--------------------+
## |    other         |     32% 42/130     |     42% 25/ 59     |
## +------------------+--------------------+--------------------+
## |smoke : non-smoker|     66% 86/130     |     49% 29/ 59     |
## +------------------+--------------------+--------------------+
## |    smoker        |     34% 44/130     |     51% 30/ 59     |
## +------------------+--------------------+--------------------+
## |ht : 0            |     96% 125/130    |     88%  52/ 59    |
## +------------------+--------------------+--------------------+
## |    1             |      4%   5/130    |     12%   7/ 59    |
## +------------------+--------------------+--------------------+
## |ui : 0            |     89% 116/130    |     76%  45/ 59    |
## +------------------+--------------------+--------------------+
## |    1             |     11%  14/130    |     24%  14/ 59    |
## +------------------+--------------------+--------------------+
## |ftv : 0           |     49% 64/130     |     61% 36/ 59     |
## +------------------+--------------------+--------------------+
## |    1             |     28% 36/130     |     19% 11/ 59     |
## +------------------+--------------------+--------------------+
## |    2             |     23% 30/130     |     20% 12/ 59     |
## +------------------+--------------------+--------------------+

Une synthèse graphique de ces mêmes données est proposée ci-dessous :

d <- aggregate(low ~ race + cut2(lwt, g = 3) + smoke, data = birthwt, sum)
d$pct <- d$low / nrow(birthwt)
names(d)[2] <- "age"
p <- ggplot(data = d, aes(x = pct, y = age)) +
       geom_line(aes(group = age), color = grey(.7)) +
       geom_point(aes(color = smoke)) +
       scale_color_manual("", values = c("steelblue", "orange")) +
       facet_wrap(~ race, ncol = 3) +
       labs(x = "Proportion weight < 2.5 kg", y = "Mother age")
p

Modélisation

Le modèle de régression logistique s’écrit :

fm <- low ~ lwt + race + smoke + ui
m <- glm(fm, data = birthwt, family = binomial)
summary(m)

## 
## Call:
## glm(formula = fm, family = binomial, data = birthwt)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.647  -0.863  -0.559   1.100   2.079  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.3866     0.8979   -0.43   0.6668
## lwt          -0.0266     0.0141   -1.89   0.0587
## raceblack     1.3104     0.5097    2.57   0.0101
## raceother     0.9557     0.4173    2.29   0.0220
## smokesmoker   1.0292     0.3820    2.69   0.0071
## ui            0.7785     0.4392    1.77   0.0763
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 211.91  on 183  degrees of freedom
## AIC: 223.9
## 
## Number of Fisher Scoring iterations: 4

Une fois les paramètres du modèle estimés, il est possible de calculer les intervalles de confiance ou transformer les coefficients de régression sous forme d’odds-ratio :

confint(m)

## Waiting for profiling to be done...

##                  2.5 %       97.5 %
## (Intercept) -2.1103878  1.433976583
## lwt         -0.0560017 -0.000365032
## raceblack    0.3131717  2.327859864
## raceother    0.1518227  1.797076721
## smokesmoker  0.2938670  1.800202392
## ui          -0.0894648  1.645071327

exp(coef(m)["smokesmoker"])

## smokesmoker 
##     2.79874

Voici les prédictions de ce modèle dans le cas “le plus défavorable” où smoke=smoker et ui=1 :

d <- expand.grid(lwt = seq(40, 100), race = factor(levels(birthwt$race)),
                 smoke = "smoker", ui = 1)
d$yhat <- predict(m, d, type = "response")
p <- ggplot(data = d, aes(x = lwt, y = yhat, color = race)) +
       geom_line(aes(group = race), size = 1) +
       scale_color_ipsum() +
       guides(color = FALSE) +
       labs(x = "Mother weight (kg)", y = "Pr(low = 1)", caption = "Predicted response curves")
direct.label(p + aes(label = race), method = "smart.grid")

Le package rms offre la commande lrm pour construire des modèles de régression logistique (cas binomial ou variables ordinales). Cette commande fournit des indices supplémentaires concernant la qualité d’ajustement et la qualité discriminante du modèle, et des procédures de post-estimation identiques à celles décrites dans le cas de la régression linéaire :

ddist <- with(birthwt, datadist(lwt, race, smoke, ui))  ## mandatory for Predict with ggplot, see below
options(datadist = "ddist")
m <- lrm(fm, data = birthwt, x = TRUE, y = TRUE)
m

## Logistic Regression Model
##  
##  lrm(formula = fm, data = birthwt, x = TRUE, y = TRUE)
##  
##                       Model Likelihood     Discrimination    Rank Discrim.    
##                          Ratio Test           Indexes           Indexes       
##  Obs           189    LR chi2     22.76    R2       0.160    C       0.702    
##   0            130    d.f.            5    g        0.927    Dxy     0.404    
##   1             59    Pr(> chi2) 0.0004    gr       2.527    gamma   0.406    
##  max |deriv| 2e-06                         gp       0.177    tau-a   0.174    
##                                            Brier    0.192                     
##  
##               Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept    -0.3866 0.8979 -0.43  0.6668  
##  lwt          -0.0266 0.0141 -1.89  0.0587  
##  race=black    1.3104 0.5097  2.57  0.0101  
##  race=other    0.9557 0.4173  2.29  0.0220  
##  smoke=smoker  1.0292 0.3820  2.69  0.0071  
##  ui            0.7785 0.4392  1.77  0.0763  
##

r <- anova(m)

Outre la possibilité d’utiliser Predict (P majuscule) pour former des prédictions marginales avec ou sans ajustement, il existe une interface avec le package ggplot2. Voici un exemple d’application :

ggplot(Predict(m, lwt, race="white", smoke, ui=1), anova = r, pval = TRUE)

Enfin, en plus des options disponibles pour la calibration d’un modèle pronostique (calibrate'), rms fournit la commande validate qui permet d’évaluer la stabilité des paramètres par validation croisée (par défaut, bootstrap) :

validate(m)

##           index.orig training    test optimism index.corrected  n
## Dxy           0.4039   0.4253  0.3776   0.0477          0.3562 40
## R2            0.1595   0.1783  0.1376   0.0407          0.1188 40
## Intercept     0.0000   0.0000 -0.0716   0.0716         -0.0716 40
## Slope         1.0000   1.0000  0.8750   0.1250          0.8750 40
## Emax          0.0000   0.0000  0.0421   0.0421          0.0421 40
## D             0.1151   0.1310  0.0977   0.0333          0.0818 40
## U            -0.0106  -0.0106  0.0028  -0.0133          0.0028 40
## Q             0.1257   0.1416  0.0949   0.0467          0.0790 40
## B             0.1919   0.1871  0.1978  -0.0107          0.2026 40
## g             0.9271   0.9983  0.8426   0.1557          0.7714 40
## gp            0.1767   0.1839  0.1641   0.0198          0.1569 40

Cas des données groupées

À titre d’illustration, voici un exemple d’analyse de données médicales pour lesquelles on dispose de la pression systolique (bpress) et de la concentration sanguine en cholestérol (chol) chez des patients ayant eu ou non un infarctus du myocarde (hdis). Le chargement des données ne pose pas de difficulté car les données sont déjà au format R :

load("data/hdis.rda")
str(hdis)

## 'data.frame':    56 obs. of  4 variables:
##  $ bpress: Factor w/ 8 levels "<117","117-126",..: 1 1 1 1 1 1 1 2 2 2 ...
##  $ chol  : Factor w/ 7 levels "<200","200-209",..: 1 2 3 4 5 6 7 1 2 3 ...
##  $ hdis  : int  2 0 0 0 0 1 0 0 2 1 ...
##  $ total : int  53 21 15 20 14 22 11 66 27 25 ...

Voici un résumé de la structure de données :

round(xtabs(hdis/total ~ bpress + chol, data = hdis), 2)

##          chol
## bpress    <200 200-209 210-219 220-244 245-259 260-284 >284
##   <117    0.04    0.00    0.00    0.00    0.00    0.05 0.00
##   117-126 0.00    0.07    0.04    0.12    0.00    0.23 0.05
##   127-136 0.03    0.00    0.10    0.02    0.00    0.08 0.14
##   137-146 0.02    0.00    0.00    0.07    0.13    0.06 0.17
##   147-156 0.05    0.00    0.00    0.10    0.11    0.25 0.06
##   157-166 0.08    0.00    0.00    0.07    0.00    0.15 0.33
##   167-186 0.14    0.00    0.00    0.07    0.40    0.38 0.21
##   >186    0.20    0.00    0.50    0.10    0.14    0.14 0.14

Pour simplifier l’analyse, on va remplacer les catégories de pression systolique par leur centre de classe et analyser cette variable en mode numérique. La deuxième étape consiste donc à agréger les données sur les différents niveaux de cholestérol :

labs <- c(111.5,121.5,131.5,141.5,151.5,161.5,176.5,191.5)
hdis$bpress <- rep(labs, each = 7)
d <- aggregate(cbind(hdis, total) ~ bpress, data = hdis, sum)

Voici comment calculer les proportions de cas :

d$prop <- d$hdis/d$total

Enfin, voici le modèle de régression logistique où, dans le cas des données groupées, il est nécessaire d’indiquer la variable qui code l’événement positif (celui qui nous intéresse) et celle qui code l’événement négatif (son complémentaire) :

m <- glm(cbind(hdis, total-hdis) ~ bpress, data = d, family = binomial)
d$yhat <- predict(m, type = "response")
d

##   bpress hdis total      prop      yhat
## 1  111.5    3   156 0.0192308 0.0333004
## 2  121.5   17   252 0.0674603 0.0420903
## 3  131.5   12   284 0.0422535 0.0530730
## 4  141.5   16   271 0.0590406 0.0667218
## 5  151.5   12   139 0.0863309 0.0835709
## 6  161.5    8    85 0.0941176 0.1041998
## 7  176.5   16    99 0.1616162 0.1435229
## 8  191.5    8    43 0.1860465 0.1944642

Les proportions observées et les valeurs prédites sont représentées dans la figure ci-dessous, la fonction f() permettant simplement de passer de l’échelle du log odds aux proportions pour le modèle m :

f <- function(x) 1/(1 + exp(-(coef(m)[1] + coef(m)[2]*x)))
p <- ggplot(data = d, aes(x = bpress, y = hdis/total)) +
       geom_point() +
       stat_function(fun = f, col = "lightcoral", size = 1) +
       labs(x = "Blodd pressure (mm Hg)", y = "Proportion CHD", caption = "Observed vs. fitted proportions")
p

Régression de Poisson

La régression de poisson permet de modéliser le log d’un taux ou une variation sous forme d’une combinaison linéaire de $p$ variables explicatives : \[ \log\big(\lambda\mid x_1, x_2, \dots\big)=\lambda_0+\beta_1x_1+\beta_2x_2+\dots, \] où $\lambda_0$ représente le taux attendu pour un individu pour lequel tous les cofacteurs valent 0. Sur l’échelle naturelle, le modèle prend, comme dans le cas de la régression logistique, une forme multiplicative où $\lambda = \lambda_0 \times \text{RR}_1^{x_1} \times \text{RR}_2^{x_2} \times \dots$.

Chargement des données

La régression de Poisson peut être utilisée dans plusieurs contexte différents : en lien avec la régression logistique (conevrgence de la loi de Poisson vers la loi binomiale), dans le cas du calcul de risque relatif dans les modèles de survie, dans le cadre des procédures de standardisation indirecte, ou plus simplement dans le cas des données de comptage. C’est ce dernier cas qui nous intéresse avec les données schoolinf.rda qui indiquent le nombre de jours après lequel on a diagnostiqué chez des étudiants la survenue d’une maladie infectieuse.

load("data/schoolinf.rda")
str(schoolinf)

## 'data.frame':    109 obs. of  2 variables:
##  $ days : int  1 2 3 3 4 4 4 6 7 8 ...
##  $ count: int  6 8 12 9 3 3 11 5 7 3 ...

Description des variables

p <- ggplot(data = schoolinf, aes(x = days, y = count)) +
       geom_point() +
       geom_smooth(method = "loess", color = "lightcoral") +
       labs(x = "No. days", y = "No. students diagnosed")
p

Modélisation

Le modèle de Poisson s’écrit comme dans le cas du modèle logistique, en adaptant simplement le type de famille dans la fonciton glm:

m1 <- glm(count ~ days, data = schoolinf, family = poisson)
summary(m1)

## 
## Call:
## glm(formula = count ~ days, family = poisson, data = schoolinf)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0048  -0.8572  -0.0933   0.6397   1.7370  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  1.99023    0.08394    23.7   <2e-16
## days        -0.01746    0.00173   -10.1   <2e-16
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 215.36  on 108  degrees of freedom
## Residual deviance: 101.17  on 107  degrees of freedom
## AIC: 393.1
## 
## Number of Fisher Scoring iterations: 5

On voit que malgré la significativité du coefficient de régression, la variance résiduelle est très élevée (trop proche des degrés de liberté), ce qui suggère un problème de sur-dispersion. Voici une solution pour prendre en compte cette variance non expliquée :

m2 <- glm(count ~ days, data = schoolinf, family = quasipoisson)
summary(m2)

## 
## Call:
## glm(formula = count ~ days, family = quasipoisson, data = schoolinf)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0048  -0.8572  -0.0933   0.6397   1.7370  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.99023    0.07479    26.6   <2e-16
## days        -0.01746    0.00154   -11.3   <2e-16
## 
## (Dispersion parameter for quasipoisson family taken to be 0.793944)
## 
##     Null deviance: 215.36  on 108  degrees of freedom
## Residual deviance: 101.17  on 107  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

On voit qu’il s’agit en fait d’un problème de sous-dispersion (le paramètre de dispersion, une fois estimé et non fixé à 1, est < 1).

Voici les prédictions de ce modèle :

d <- data.frame(days = seq(0, 150, 2))
d$yhat <- predict(m2, d, type = "response")
p <- ggplot(data = schoolinf, aes(x = days, y = count)) +
       geom_point() +
       geom_line(data = d, aes(x = days, y = yhat), color = "lightcoral", size = 1) +
       labs(x = "No. days", y = "No. students diagnosed")
p

Modèle de Cox

Le modèle de régression de Cox est un modèle semi-paramétrique simple dans lequel on modélise le risque ou la fréquence instantanée de l’événement d’intérêt, appelé “hazard”, en fonction de variables explicatives sous la forme : \[ \log\big(h(t\mid x_1,x_2,\dots)\big) = \log\big(h_0(t)\big) + \beta_1x_1 + \beta_2x_2 + \dots. \] À la différence du modèle de Poisson, le risque de base, $h_0(t)$, peut varier librement d’un individu à l’autre et il n’est pas estimé : on s’intéresse juste au rapport de risques, $\frac{h(t\mid X=x+1)}{h(t\mid X=x)}=\frac{h_0(t)e^{\beta_1(x+1)}}{h_0(t)e^{\beta_1x}}=e^{\beta_1}$.

Chargement des données

Le package survival permet de traiter le cas de ce type de données, appelées données de survie, pour lesquelles sont associées une variable codant une durée ou une date chronologique et une variable codant un événement binaire (0 = l’événement n’a pas été observé à la dernière date d’observation, 1 = l’événement a été observé ; ici, l’événement est le décès). Voici comment procéder sour R pour définir cette variable composite :

data(lung, package = "survival")
lung$sex <- factor(lung$sex, levels = 1:2, labels = c("Male", "Female"))
st <- with(lung, Surv(time=time, event=status))
head(st)

## [1]  306   455  1010+  210   883  1022+

Description des variables

Toutes les commandes R doivent utiliser la variable composite définie par Surv, raison pour laquelle il est plus commode de la stocker une fois pour toutes dans une variable, ici st. La médiane de survie et son intervalle de confiance pour l’ensemble des individus s’obtient de la manière suivante :

s <- survfit(st ~ 1, data=lung)
s

## Call: survfit(formula = st ~ 1, data = lung)
## 
##       n  events  median 0.95LCL 0.95UCL 
##     228     165     310     285     363

Pour obtenir la médiane de survie dans deux groupes, on remplacera le terme d’intercept (~ 1) par le facteur correspondant. En réalité, survfit permet de construire la table de mortalité, que l’on peut afficher entièrement ou partiellement en appliquant la méthode summary.survfit :

summary(s, times=seq(1, 200, by=20))

## Call: survfit(formula = st ~ 1, data = lung)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    228       0    1.000  0.0000        1.000        1.000
##    21    220       8    0.965  0.0122        0.941        0.989
##    41    217       3    0.952  0.0142        0.924        0.980
##    61    211       7    0.921  0.0179        0.887        0.957
##    81    205       7    0.890  0.0207        0.851        0.932
##   101    196       6    0.864  0.0227        0.821        0.910
##   121    189       6    0.837  0.0245        0.791        0.887
##   141    184       5    0.815  0.0257        0.766        0.867
##   161    176       8    0.780  0.0275        0.728        0.836
##   181    159      15    0.713  0.0301        0.656        0.774

Enfin, on peut résumer graphiquement cette table de mortalité à l’aide d’une courbe de Kaplan-Meier (KM) :

s <- survfit(st ~ sex, data=lung)
p <- autoplot(s, censor = TRUE, censor.colour = grey(.5)) +
       scale_color_manual("", values = c("steelblue", "orange")) +
       scale_fill_manual("", values = c("steelblue", "orange")) +
       guides(fill = FALSE) +
       theme(legend.position = c(.9, .9)) +
       labs(x = "Followup time (hr.)", y = "Probability survival")
p

Modélisation

Le test du log-rank permettant de tester l’égalité de deux courbes de KM s’obtient avec la commande survdiff :

survdiff(st ~ sex, data=lung)

## Call:
## survdiff(formula = st ~ sex, data = lung)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=Male   138      112     91.6      4.55      10.3
## sex=Female  90       53     73.4      5.68      10.3
## 
##  Chisq= 10.3  on 1 degrees of freedom, p= 0.00131

Quant au modèle de Cox, la syntaxe est la même que dans le cas des autres modèles, à ceci près que la variable réponse est la variable composite définie plus haut, st :

m1 <- coxph(st ~ sex, data=lung)
m1

## Call:
## coxph(formula = st ~ sex, data = lung)
## 
##             coef exp(coef) se(coef)     z      p
## sexFemale -0.531     0.588    0.167 -3.18 0.0015
## 
## Likelihood ratio test=10.6  on 1 df, p=0.00111
## n= 228, number of events= 165

Il est possible de traiter un cofacteur comme une variable de stratification, de sorte que ses paramètres sont estimés mais ne contribuent pas au degré de liberté du modèle :

m2 <- coxph(st ~ sex + strata(age), data=lung)
m2

## Call:
## coxph(formula = st ~ sex + strata(age), data = lung)
## 
##             coef exp(coef) se(coef)     z      p
## sexFemale -0.645     0.525    0.212 -3.04 0.0024
## 
## Likelihood ratio test=9.79  on 1 df, p=0.00176
## n= 228, number of events= 165

Pour vérifier l’hypothèse de risques proportionnel, il est possible d’utiliser des méthodes graphiques (en examinant la distribution des résidus de Schoenfeld) ou d’utiliser la procédure de test des résidus offerte par la commande cox.zph.

m3 <- coxph(st ~ sex + age, data=lung)
cox.zph(m3)

##               rho chisq     p
## sexFemale  0.1236 2.452 0.117
## age       -0.0275 0.129 0.719
## GLOBAL         NA 2.651 0.266

Exercices d’application

Régression logistique

Les données clslowbwt.dta (format Stata) contiennent des données sur des poids de naissance comparables à celles analysées plus haut. Cette fois, les observations ne sont toutes indépendantes car il y a plusieurs enfants par mère. On ne dispose pas non plus des variables telles que ht ou ui.

Décrire les données, en particulier le poids des bébés (bwt et low) en fonction des variables d’intérêt (lwt, race, smoke).
Estimer les paramètres d’un modèle prédisant low en fonction de ces 3 variables.
Comparer les résultats (en particulier les erreurs standard) avec ceux obtenus en utilisant un estimateur robuste de la variance (utiliser le package sandwich et la fonction coeftest ou la fonction robcov du package rms). La variable de cluster sera évidemment l’identifiant des mères.

Régression logistique pénalisée

Les données leukemia.rda contiennent le statut du patient (première colonne) et un ensemble de mesures normalisées d’expression génique (colonne 2 à 3052). L’objectif est de comparer une approche de sélection de variable par tests univariés et par modèle de régerssion pénalisée.

Réaliser un test de Student pour chaque variable X* et collecter le degré de signification du test. Si le résultat du test est stocké dans une variable, il est possible d’accéder à celui-ci en utilisant l’objet $p.value. Afficher la distribution de l’ensemble des p-valeurs sous forme d’histogramme. Combien y’a t-il de tests significatifs avec et sans correction pour les tests multiples ?
À l’aide du package glmnet, réaliser une régression logistique avec une pénalisation de type lasso (alpha=1). Comparer les variables sélectionnées avec celles retenues en (1).

Données de survie (1)

Les données compliance2.dta (format Stata) contiennent les données d’une étude sur la mortalité de sujets ayant accepté ou non de participer à un programme de dépistage de l’anévrisme de l’aorte abdominale. Au total il y a 555 participants. On souhaite comparer la survie des participants entre les deux groupes (particp). Les informations temporelles consistent en une date d’entrée (randate) et une date de dernier point (enddate) et il faudra donc définir un intervalle (time=, time2=) à partir de ces deux variables.

Construire la courbe de Kaplan Meier pour les deux groupes de cet échantillon et calculer la valeur du 10e percentile de survie avec son intervalle de confiance à 95 % pour les deux groupes de sujets.
Estimer les paramètres d’un modèle de Cox en (a) ajustant ou (b) stratifiant sur l’âge codé en 5 classes (ranagr).

Données de survie (2)

Dans un essai contre placebo sur la cirrhose biliaire, la D-penicillamine (DPCA) a été introduite dans le bras actif sur une cohorte de 312 patients. Au total, 154 patients ont été randomisés dans le bras actif (variable traitement, rx, 1=Placebo, 2=DPCA). Un ensemble de données telles que l’âge, des données biologiques et signes cliniques variés incluant le niveau de bilirubine sérique (bilirub) sont disponibles dans le fichier pbc.rda. 1 Le status du patient est enregistré dans la variable status (0=vivant, 1=décédé) et la durée de suivi (years) représente le temps écoulé en années depuis la date de diagnostic.

Les instructions suivantes permettent de charger les données et d’encoder les variables catégorielles :

load("data/pbc.rda")
pbc$rx <- factor(pbc$rx, levels = 1:2, labels = c("Placebo", "DPCA"))
pbc$sex <- factor(pbc$sex, levels = 0:1, labels = c("M","F"))

Voici un tableau descriptif des variables d’intérêt :

summary(status ~ rx + sex + cut2(age, g = 4), data = pbc, fun = table)

## status     N= 312 
## 
## +----------------+-----------+---+---+---+
## |                |           |N  |0  |1  |
## +----------------+-----------+---+---+---+
## |rx              |Placebo    |158| 93| 65|
## |                |DPCA       |154| 94| 60|
## +----------------+-----------+---+---+---+
## |sex             |M          | 36| 14| 22|
## |                |F          |276|173|103|
## +----------------+-----------+---+---+---+
## |cut2(age, g = 4)|[26.3,42.3)| 78| 59| 19|
## |                |[42.3,49.8)| 78| 52| 26|
## |                |[49.8,56.8)| 78| 41| 37|
## |                |[56.8,78.4]| 78| 35| 43|
## +----------------+-----------+---+---+---+
## |Overall         |           |312|187|125|
## +----------------+-----------+---+---+---+

Afficher la distribution des durées de suivi des 312 patients, en faisant apparaître distinctement les individus décédés. Calculer le temps médian (en années) de suivi pour chacun des deux groupes de traitement. Combien y’a t-il d’événements positifs au-delà de 10.5 années et quel est le sexe de ces patients ?
Calculer la médiane de survie et son intervalle de confiance à 95 % pour chaque groupe de sujets et afficher les courbes de survie correspondantes.
Effectuer un test du log-rank en considérant comme prédicteur le facteur rx. Comparer avec les résultats obtenus à partir d’un modèle de Cox.

Méthodes biostatistiques

Modèles de régression pour données non gaussiennes

Régression logistique

Chargement des données

Description des variables

Modélisation

Cas des données groupées

Régression de Poisson

Chargement des données

Description des variables

Modélisation

Modèle de Cox

Chargement des données

Description des variables

Modélisation

Exercices d’application

Régression logistique

Régression logistique pénalisée

Données de survie (1)

Données de survie (2)