Analyse de survie
Lazare Sân-ib DA
18/12/2024
Nous allons utiliser deux packages R :
• survival pour le calcul d’analyses de survie
• survminer pour résumer et visualiser les résultats de l’analyse de survie
library("survival")
library("survminer")## Loading required package: ggplot2## Loading required package: ggpubr##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myelomaNous utiliserons les données sur le cancer du poumon disponibles dans le package survival.
Préparation des données
data(cancer, package="survival")
head(lung)## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 0QUELQUES DETAILS SUR LES VARIABLES DE LA BASE
• inst : code de l’établissement
• time : Temps de survie en jours
• status : statut de censure 1=censuré, 2=mort
• age : âge en années
• sex : Masculin=1 Féminin=2
• ph.ecog : score de performance ECOG (0=bon 5=mort)
• ph.karno : score de performance de Karnofsky (mauvais=0-bon=100) noté par un médecin
• pat.karno : score de performance de Karnofsky évalué par le patient
• repas.cal : Calories consommées aux repas
• wt.loss : Perte de poids au cours des six derniers mois
Caalcul le temps de survie
Dans le cas où nous avons des données avec de dates de debut et de fin ,la prémière etape consiste à s’assurer qu’elles sont du type date
library("tibble")
date_ex <-
tibble(
sx_date = c("2007-06-22", "2004-02-13", "2010-10-27"),
last_fup_date = c("2017-04-15", "2018-07-04", "2016-10-31")
)
date_ex## # A tibble: 3 × 2
## sx_date last_fup_date
## <chr> <chr>
## 1 2007-06-22 2017-04-15
## 2 2004-02-13 2018-07-04
## 3 2010-10-27 2016-10-31Création d’un objet de survie
surv_objet<-Surv(lung$time, lung$status)[1:10]
surv_objet## [1] 306 455 1010+ 210 883 1022+ 310 361 218 166ESTIMATION DE KAPLAN-MEIER
La fonction survfit() dans le package de survival peut être utilisée pour calculer l’estimation de survie de Kaplan-Meier. Ses principaux arguments incluent :
• un objet de survie créé à l’aide de la fonction Surv()
• et l’ensemble de données contenant les variables.
Fonction de survie
Son premier argument est formula. Le côté gauche de cette formule spécifie les informations sur les temps de survie à l’aide de la fonction Surv(), et le côté droit est utilisé pour spécifier les variables de regroupement. Argument data spécifie le bloc de données qui contient les variables d’intérêt(dans notre cas lung).
fit <- survfit(Surv(time, status) ~ sex, data = lung)
print(fit)## Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
##
## n events median 0.95LCL 0.95UCL
## sex=1 138 112 270 212 310
## sex=2 90 53 426 348 550MEDIAN ET MOYENNE ESTIMES
print(survfit(Surv(time, status) ~ sex, data = lung),print.rmean = TRUE)## Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
##
## n events rmean* se(rmean) median 0.95LCL 0.95UCL
## sex=1 138 112 326 22.9 270 212 310
## sex=2 90 53 461 34.7 426 348 550
## * restricted mean with upper limit = 1022Les temps de survie médians pour chaque groupe représentent le moment auquel la probabilité de survie, S(t), est de 0,5
QUANTILE
On peut utiliser la méthode quantile() pour calculer les temps de suivi correspondants auxquels la probabilité de survie prend une valeur spécifique. Par exemple, nous voulons trouver à combien de jours la probabilité de survie est égale à 0,7 et à combien de jours elle est égale à 0,6 ;le code est :
quantile(fit, probs = 1 - c(0.7, 0.6))## $quantile
## 30 40
## sex=1 166 202
## sex=2 285 348
##
## $lower
## 30 40
## sex=1 135 176
## sex=2 201 285
##
## $upper
## 30 40
## sex=1 197 269
## sex=2 351 444NB Dans l’argument probs de quantile()nous devons spécifier un moins nos probabilités de survie cibles ; en effet, la fonction selon la convention de la fonction de distribution cumulative (CDF) et la CDF est égale à un moins la probabilité de survie. De plus, notez également que nous obtenons les bornes inférieure et supérieure des intervalles de confiance à 95% pour les quantiles que nous avons demandés.
differentes composantes accessibles
d <- data.frame(time = fit$time,
n.risk = fit$n.risk,
n.event = fit$n.event,
n.censor = fit$n.censor,
surv = fit$surv,
upper = fit$upper,
lower = fit$lower
)
head(d)## time n.risk n.event n.censor surv upper lower
## 1 11 138 3 0 0.9782609 1.0000000 0.9542301
## 2 12 135 1 0 0.9710145 0.9994124 0.9434235
## 3 13 134 2 0 0.9565217 0.9911586 0.9230952
## 4 15 132 1 0 0.9492754 0.9866017 0.9133612
## 5 26 131 1 0 0.9420290 0.9818365 0.9038355
## 6 30 130 1 0 0.9347826 0.9768989 0.8944820En utilisant la méthode summary() et son argument times, nous pouvons obtenir les probabilités de survie à des moments de suivi spécifiques.
• n : nombre total de sujets dans chaque courbe.
• time : les points temporels sur la courbe.
• n.risk : le nombre de sujets à risque au temps t
• n.event : le nombre d’événements qui se sont produits à l’instant t.
• n.censor : le nombre de sujets censurés, qui sortent de l’ensemble de risques, sans événement, à l’instant t.
summary(fit)$time## [1] 11 12 13 15 26 30 31 53 54 59 60 65 71 81 88 92 93 95
## [19] 105 107 110 116 118 131 132 135 142 144 147 156 163 166 170 175 176 177
## [37] 179 180 181 183 189 197 202 207 210 212 218 222 223 229 230 239 246 267
## [55] 269 270 283 284 285 286 288 291 301 303 306 310 320 329 337 353 363 364
## [73] 371 387 390 394 428 429 442 455 457 460 477 519 524 533 558 567 574 583
## [91] 613 624 643 655 689 707 791 814 883 5 60 61 62 79 81 95 107 122
## [109] 145 153 166 167 182 186 194 199 201 208 226 239 245 268 285 293 305 310
## [127] 340 345 348 350 351 361 363 371 426 433 444 450 473 520 524 550 641 654
## [145] 687 705 728 731 735 765summary(fit)$surv #estimation de kaplan## [1] 0.97826087 0.97101449 0.95652174 0.94927536 0.94202899 0.93478261
## [7] 0.92753623 0.91304348 0.90579710 0.89855072 0.89130435 0.87681159
## [13] 0.86956522 0.86231884 0.84782609 0.84057971 0.83333333 0.82608696
## [19] 0.81884058 0.81159420 0.80434783 0.79710145 0.78985507 0.78260870
## [25] 0.76811594 0.76086957 0.75362319 0.74637681 0.73913043 0.72463768
## [31] 0.70289855 0.69565217 0.68840580 0.68108233 0.67375887 0.66643540
## [37] 0.65178847 0.64446500 0.62981807 0.62249460 0.61499467 0.60730724
## [43] 0.59952125 0.59173526 0.58394926 0.57616327 0.56837728 0.56048316
## [49] 0.55247625 0.54423034 0.53598442 0.52760967 0.51923491 0.51072286
## [55] 0.50221082 0.49369877 0.48503739 0.47637600 0.46755423 0.45873245
## [61] 0.44991067 0.44108889 0.43189954 0.42251042 0.41290791 0.40330540
## [67] 0.39370289 0.38410038 0.37449787 0.35529285 0.34569034 0.33608783
## [73] 0.32648533 0.31688282 0.30728031 0.29767780 0.28741304 0.27714829
## [79] 0.26688354 0.25620820 0.24553286 0.23437227 0.22321169 0.21205110
## [85] 0.20089052 0.18972994 0.17856935 0.16740877 0.15624818 0.14508760
## [91] 0.13392701 0.12276643 0.11160584 0.10044526 0.08928468 0.07812409
## [97] 0.06696351 0.05357081 0.03571387 0.98888889 0.97777778 0.96666667
## [103] 0.95555556 0.94444444 0.93333333 0.92208835 0.91070455 0.89932074
## [109] 0.87655313 0.86516932 0.85378551 0.84240171 0.83053689 0.81867208
## [115] 0.80663278 0.79459349 0.77051490 0.75808724 0.74523830 0.73216395
## [121] 0.71860536 0.70451505 0.68952537 0.67420259 0.65852346 0.64284433
## [127] 0.62636114 0.60987795 0.59339476 0.57691157 0.56042839 0.54344571
## [133] 0.52646303 0.50891426 0.48934064 0.46976701 0.45019339 0.43061976
## [139] 0.41104614 0.38941213 0.36777812 0.34325958 0.31205416 0.28084875
## [145] 0.24964333 0.21843791 0.18723250 0.15602708 0.12482167 0.08321444str(summary(fit)) #plus de détails sur le contenu de la liste## List of 20
## $ n : int [1:2] 138 90
## $ time : num [1:150] 11 12 13 15 26 30 31 53 54 59 ...
## $ n.risk : num [1:150] 138 135 134 132 131 130 129 128 126 125 ...
## $ n.event : num [1:150] 3 1 2 1 1 1 1 2 1 1 ...
## $ n.censor : num [1:150] 0 0 0 0 0 0 0 0 0 0 ...
## $ surv : num [1:150] 0.978 0.971 0.957 0.949 0.942 ...
## $ std.err : num [1:150] 0.0124 0.0143 0.0174 0.0187 0.0199 ...
## $ cumhaz : num [1:150] 0.0217 0.0291 0.0441 0.0516 0.0593 ...
## $ std.chaz : num [1:150] 0.0126 0.0146 0.018 0.0195 0.021 ...
## $ strata : Factor w/ 2 levels "sex=1","sex=2": 1 1 1 1 1 1 1 1 1 1 ...
## $ type : chr "right"
## $ logse : logi TRUE
## $ conf.int : num 0.95
## $ conf.type : chr "log"
## $ lower : num [1:150] 0.954 0.943 0.923 0.913 0.904 ...
## $ upper : num [1:150] 1 0.999 0.991 0.987 0.982 ...
## $ t0 : num 0
## $ call : language survfit(formula = Surv(time, status) ~ sex, data = lung)
## $ table : num [1:2, 1:9] 138 90 138 90 138 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "sex=1" "sex=2"
## .. ..$ : chr [1:9] "records" "n.max" "n.start" "events" ...
## $ rmean.endtime: num [1:2] 1022 1022
## - attr(*, "class")= chr "summary.survfit"Courbe de survie
Nous utiliserons la fonction ‘ggsurvplot()‘dans le package Survminer R pour produire les courbes de survie des deux groupes de sujets
ggsurvplot(fit,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"))## Warning in geom_segment(aes(x = 0, y = max(y2), xend = max(x1), yend = max(y2)), : All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
## All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
## All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
## All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
INTERPRETATION DU GRAPHIQUE L’axe horizontal (axe des x) représente le temps en jours et l’axe vertical (axe des y) montre la probabilité de survie ou la proportion de personnes survivantes. Les lignes représentent les courbes de survie des deux groupes. Une chute verticale dans les courbes indique un événement. La coche verticale sur les courbes signifie qu’un patient a été censuré à ce moment.
• Au temps zéro, la probabilité de survie est de 1,0 (ou 100 % des participants sont en vie).
• Au temps 250, la probabilité de survie est d’environ 0,55 (ou 55%) pour le sexe=1 et de 0,75 (ou 75%) pour le sexe=2.
• La survie médiane est d’environ 270 jours pour le sexe=1 et de 426 jours pour le sexe=2, suggérant une bonne survie pour le sexe=2 par rapport au sexe=1.Les temps de survie médians pour chaque groupe représentent le moment auquel la probabilité de survie, S(t), est de 0,5. Le temps de survie médian pour le sexe = 1 (groupe masculin) est de 270 jours, contre 426 jours pour le sexe = 2 (féminin). Il semble y avoir un avantage de survie pour les femmes atteintes d’un cancer du poumon par rapport aux hommes. Cependant, pour évaluer si cette différence est statistiquement significative, il faut un test statistique formel, un sujet qui est abordé dans les sections suivantes
Fonctions de risques cumulatif
La fonction de risque cumulatif et la fonction de survie sont liées par la relation suivante : \(𝑆(𝑡) = 𝑒𝑥𝑝(−𝐻(𝑡))\). Il correspond au nombre d’événements qui seraient attendus pour chaque individu au temps t si l’événement était un processus répétable .Lorsque les deux courbes sont proportionnelles l’une à l’autre (c’est-à-dire qu’elles s’éloignent régulièrement l’une de l’autre) on dit qu’il y a une différence de survie significative entre les deux groupes.
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
fun = "cumhaz")
COMPARAISON DES COURBES DE SURVIES
Test de ‘log-rank‘
Le test du log-rank est la méthode la plus largement utilisée pour comparer deux ou plusieurs courbes de survie (c’est-à-dire pour tester l’hypothèse si les fonctions de survie de différents groupes de sujets diffèrent de manière statistiquement significative). L’hypothèse nulle est qu’il n’y a pas de différence de survie entre les deux groupes.
Le test du log-rankz est un test non paramétrique, qui ne fait aucune hypothèse sur les distributions de survie.La fonction 𝑠𝑢𝑟𝑣𝑑𝑖𝑓𝑓() dans le package de survival peut être utilisée pour calculer le test log-rank comparant deux ou plusieurs courbes de survie. Pour tester avec le test du log-rank s’il existe des différences dans les taux de survie dans l’ensemble de données Lung entre les hommes et les femmes, nous utilisons le code :
surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung)
surv_diff## Call:
## survdiff(formula = Surv(time, status) ~ sex, data = lung)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138 112 91.6 4.55 10.3
## sex=2 90 53 73.4 5.68 10.3
##
## Chisq= 10.3 on 1 degrees of freedom, p= 0.001Le test du log-rank pour la différence de survie donne une valeur p de p = 0,001, indiquant que les groupes de sexe diffèrent significativement en termes de survie. Pour evaluer la même hypothèse on peut utiliser le test de Peto & Peto Gehan-Wilcoxon, nous utilisons survdiff()à nouveau la fonction, mais maintenant nous définissons l’argument rho à 1
peto_peto <- survdiff(Surv(time, status == 2) ~ sex, data = lung, rho = 1)
peto_peto## Call:
## survdiff(formula = Surv(time, status == 2) ~ sex, data = lung,
## rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138 70.4 55.6 3.95 12.7
## sex=2 90 28.7 43.5 5.04 12.7
##
## Chisq= 12.7 on 1 degrees of freedom, p= 4e-04La conclusion reste la même avec une pvalue faible.
Modèle de cox
La fonction 𝑐𝑜𝑥𝑝ℎ()dans le package de survival peut être utilisée pour calculer le modèle de régression à risques proportionnels de Cox dans R. 𝑐𝑜𝑥𝑝ℎ(𝑓𝑜𝑟𝑚𝑢𝑙𝑎, 𝑑𝑎𝑡𝑎, 𝑚𝑒𝑡ℎ𝑜𝑑).Nous ajusterons la régression de Cox en utilisant les covariables suivantes : âge, sexe, ph.ecog et wt.loss.
Régression de Cox univariée
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
##
## coef exp(coef) se(coef) z p
## sex -0.5310 0.5880 0.1672 -3.176 0.00149
##
## Likelihood ratio test=10.63 on 1 df, p=0.001111
## n= 228, number of events= 165summary(res.cox)## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
##
## n= 228, number of events= 165
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex -0.5310 0.5880 0.1672 -3.176 0.00149 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 0.588 1.701 0.4237 0.816
##
## Concordance= 0.579 (se = 0.021 )
## Likelihood ratio test= 10.63 on 1 df, p=0.001
## Wald test = 10.09 on 1 df, p=0.001
## Score (logrank) test = 10.33 on 1 df, p=0.001INTERPRETATION DES RESULTATS
Signification statistique. La colonne marquée « z » donne la valeur statistique de Wald. Il correspond au rapport de chaque coefficient de régression à son erreur type (𝑧 = 𝑐𝑜𝑒𝑓/𝑠𝑒(𝑐𝑜𝑒𝑓)). La statistique de wald évalue si le bêta (𝛽) d’une variable donnée est statistiquement significativement différent de 0. D’après la sortie ci-dessus,nous pouvons conclure que les variables sexe ont des coefficients statistiquement très significatifs.
Coefficients de régression. La deuxième caractéristique à noter dans les résultats du modèle de Cox est le signe des coefficients de régression (𝑐𝑜𝑒𝑓). Un signe positif signifie que le danger (risque de décès) est plus élevé, et donc le pronostic pire, pour les sujets ayant des valeurs plus élevées de cette variable. La variable sexe est codée sous forme de vecteur numérique. 1: mâle, 2: femelle. Le résumé R du modèle de Cox donne le hazard ratio (HR) pour le deuxième groupe par rapport au premier groupe, c’est-à-dire les femmes par rapport aux hommes. Le coefficient bêta pour le sexe = -0,53 indique que les femmes ont un risque de décès plus faible (taux de survie plus faible) que les hommes, dans ces données.
Rapports de risque. Les coefficients exponentiels (𝑒𝑥𝑝(𝑐𝑜𝑒𝑓) = 𝑒𝑥𝑝(−0, 53) =0, 59), également appelés hazard ratios, donnent l’ampleur de l’effet des covariables. Par exemple, le fait d’être une femme (sexe =2) réduit le risque d’un facteur de 0,59, ou 41 %. Être une femme est associé à un bon pronostic.
Intervalles de confiance des hazard ratios. La sortie sommaire donne également des intervalles de confiance supérieurs et inférieurs à 95 % pour le hazard ratio (exp(coef)), la borne inférieure à 95 % = 0,4237, la borne supérieure à 95 % = 0,816.
Signification statistique globale du modèle. Enfin, la sortie donne des valeurs p pour trois tests alternatifs pour la signification globale du modèle: le test du rapport de vraisemblance, le test de Wald et les statistiques de logrank de score. Ces trois méthodes sont asymptotiquement équivalentes. Pour un N suffisamment grand, ils donneront des résultats similaires. Pour les petits N, ils peuvent différer quelque peu. Le test du rapport de vraisemblance a un meilleur comportement pour les échantillons de petite taille, il est donc généralement préféré.
ANALYSE DE REGRESSION DE COX MULTIVARIEE
res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
summary(res.cox)## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)
##
## n= 227, number of events= 164
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.011067 1.011128 0.009267 1.194 0.232416
## sex -0.552612 0.575445 0.167739 -3.294 0.000986 ***
## ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0111 0.9890 0.9929 1.0297
## sex 0.5754 1.7378 0.4142 0.7994
## ph.ecog 1.5900 0.6289 1.2727 1.9864
##
## Concordance= 0.637 (se = 0.025 )
## Likelihood ratio test= 30.5 on 3 df, p=1e-06
## Wald test = 29.93 on 3 df, p=1e-06
## Score (logrank) test = 30.5 on 3 df, p=1e-06EXPLICATION DU MODELE
La valeur de p pour les trois tests globaux (probabilité, Wald et score) est significative, ce qui indique que le modèle est significatif. Ces tests évaluent l’hypothèse nulle omnibus selon laquelle tous les bêtas (𝛽) sont nuls. Dans l’exemple ci-dessus, les statistiques de test sont en étroite concordance et l’hypothèse nulle omnibus est fermement rejetée.Dans l’analyse de Cox multivariée, les covariables sexe et ph.ecog restent significatives (p< 0,05).
Cependant, l’âge de la covariable n’est pas significatif (p=0,23, ce qui est supérieur à 0,05). La valeur p pour le sexe est de 0,000986, avec un hazard ratio HR = exp(coef) = 0,58, indiquant une forte relation entre le sexe des patients et une diminution du risque de décès. Les hazard ratios des covariables peuvent être interprétés comme des effets multiplicatifs sur le danger. Par exemple, le fait de maintenir les autres covariables constantes, étant une femme (sexe = 2), réduit le risque d’un facteur de 0,58, ou 42 %. Nous concluons qu’être une femme est associé à un bon pronostic. De même, la valeur p pour ph.ecog est de \(4,45e-05\),avec un hazard ratio HR = 1,59, indiquant une forte relation entre la valeur ph.ecog et un risque accru de décès. En maintenant les autres covariables constantes, une valeur plus élevée de ph.ecog est associée à une faible survie. En revanche, la valeur de p pour l’âge est maintenant p = 0,23. Le hazard ratio HR = exp(coef) = 1,01, avec un intervalle de confiance à 95 % de 0,99 à 1,03. Étant donné que l’intervalle de confiance pour la FC comprend 1, ces résultats indiquent que l’âge contribue moins à la différence de HR après ajustement pour les valeurs ph.ecog et le sexe du patient, et seulement la tendance vers la signification.
Par exemple, en maintenant les autres covariables constantes, un âge supplémentaire induit un risque quotidien de décès d’un facteur 𝑒𝑥𝑝(𝛽)= 1,01, ou 1%, ce qui n’est pas une contribution significative.
MODELES DE DEFAILLANCE ACCELERE
La fonction qui correspond aux modèles AFT(Accelerated Failure Times) du package de survie est survreg(). Son premier argument est une formule et a une syntaxe similaire à la fonction survfit(). L’ argument dist spécifie la distribution des temps de survie ( Remarque : l’ argument dist spécifie la distribution des temps de survie et non les temps de survie du journal). La distribution par défaut (c’est-à-dire si vous ne spécifiez pas l’ argument dist vous-même) est la distribution de Weibull. Comme pour les autres fonctions d’ajustement de modèle dans R, la summary() fonction renvoie une sortie détaillée du modèle ajusté. Voici le code pour la distribution de Weibull :
fit_weibull <- survreg(Surv(time, status) ~sex + age+ph.ecog, data = lung)
summary(fit_weibull)##
## Call:
## survreg(formula = Surv(time, status) ~ sex + age + ph.ecog, data = lung)
## Value Std. Error z p
## (Intercept) 6.27344 0.45358 13.83 < 2e-16
## sex 0.40109 0.12373 3.24 0.0012
## age -0.00748 0.00676 -1.11 0.2690
## ph.ecog -0.33964 0.08348 -4.07 4.7e-05
## Log(scale) -0.31319 0.06135 -5.11 3.3e-07
##
## Scale= 0.731
##
## Weibull distribution
## Loglik(model)= -1132.4 Loglik(intercept only)= -1147.4
## Chisq= 29.98 on 3 degrees of freedom, p= 1.4e-06
## Number of Newton-Raphson Iterations: 5
## n=227 (1 observation deleted due to missingness)Pour ajuster le même modèle mais avec la distribution exponentielle, le code est :
fit_exp <- survreg(Surv(time, status) ~sex + age+ph.ecog, data = lung,
dist = "exponential")
summary(fit_exp)##
## Call:
## survreg(formula = Surv(time, status) ~ sex + age + ph.ecog, data = lung,
## dist = "exponential")
## Value Std. Error z p
## (Intercept) 6.37342 0.62076 10.27 < 2e-16
## sex 0.50906 0.16716 3.05 0.00232
## age -0.01022 0.00918 -1.11 0.26555
## ph.ecog -0.40502 0.11270 -3.59 0.00033
##
## Scale fixed at 1
##
## Exponential distribution
## Loglik(model)= -1143.6 Loglik(intercept only)= -1156.1
## Chisq= 25.12 on 3 degrees of freedom, p= 1.5e-05
## Number of Newton-Raphson Iterations: 4
## n=227 (1 observation deleted due to missingness)On distingue également les distributions log-normale, log-logistic..
QUESTION FONDAMENTALE :
si on souhaite modéliser avec une structure plus complexe devrions nous choisir un modèle paramétrique ajusté avec la fonction survreg ou un modèle de cox ajusté avec coxph. Si nous souhaitons utiliser le modèle pour la prédiction nous devons utiliser la fonction survreg car coxph n’extrapole pas au-delà de la dernière observation. «De combien le risque de décès diminue-t-il si un nouveau traitement médical est administré à un patient ?»:utilisation de Coxph «Quelle proportion de patients mourront dans 2 ans d’après les résultats des données d’une expérience qui n’a duré que 4 mois» utilisation de Survreg