Analyse de survie
Lazare Sân-ib DA
11 juillet 2023
Nous allons utiliser le jeu de données fecondite fourni par l’extension questionr. Ce jeu de données comporte trois tableaux de données : menages , femmes et enfants .
Nous souhaitons étudier ici la survie des enfants entre la naissance et l’âge de 5 ans. Dans un premier temps, nous comparerons la survie des jeunes filles et des jeunes garçons. Dans un second temps, nous procéderons à une analyse multivariée en prenant en compte les variables suivantes : • sexe de l’enfant • milieu de résidence • niveau de vie du ménage • structure du ménage • niveau d’éducation de la mère • âge de la mère à la naissance de l’enfant • enfin, une variable un peu plus compliquée, à savoir si le rang de naissance de l’enfant (second, troisième, quatrième, etc.) est supérieur au nombre idéal d’enfants selon la mère.
Nous allons préparer les données selon deux approches : soit en utilisant l’extension data.table , soit en utilisant l’extension dplyr.
library(questionr, quietly = TRUE)
data(fecondite)
lookfor(menages)## pos variable label col_type values
## 1 id_menage Identifiant du ménage dbl
## 2 taille Taille du ménage (nombre ~ dbl
## 3 sexe_chef Sexe du chef de ménage dbl+lbl [1] homme
## [2] femme
## 4 structure Structure démographique d~ dbl+lbl [0] pas d'adulte
## [1] un adulte
## [2] deux adultes de sexe ~
## [3] deux adultes de même ~
## [4] trois adultes ou plus~
## [5] adultes sans lien de ~
## 5 richesse Niveau de vie (quintiles) dbl+lbl [1] très pauvre
## [2] pauvre
## [3] moyen
## [4] riche
## [5] très richelookfor(femmes)## pos variable label col_type values
## 1 id_femme Identifiant de l'enquêtée dbl
## 2 id_menage Identifiant du ménage dbl
## 3 poids Poids statistique dbl
## 4 date_entretien Date de passation du questio~ date
## 5 date_naissance Date de naissance date
## 6 age Âge révolu (en années) à la ~ dbl
## 7 milieu Milieu de résidence dbl+lbl [1] urbain
## [2] rural
## 8 region Région de résidence dbl+lbl [1] Nord
## [2] Est
## [3] Sud
## [4] Ouest
## 9 educ Niveau d'éducation dbl+lbl [0] aucun
## [1] primaire
## [2] secondaire
## [3] supérieur
## 10 travail A un emploi ? hvn_lbl_ [0] non
## [1] oui
## [9] manquant
## 11 matri Statut matrimonial dbl+lbl [0] célibataire
## [1] mariée
## [2] en concubinage
## [3] veuve
## [4] divorcée
## [5] séparée
## 12 religion Religion dbl+lbl [1] musulmane
## [2] chrétienne
## [3] protestante
## [4] sans religion
## [5] autre
## 13 journal Lit la presse ? dbl+lbl [0] non
## [1] oui
## 14 radio Ecoute la radio ? dbl+lbl [0] non
## [1] oui
## 15 tv Regarde la télévision ? dbl+lbl [0] non
## [1] oui
## 16 nb_enf_ideal Nombre idéal d'enfants hvn_lbl_ [96] Ne sait pas
## [99] manquant
## 17 test A déjà fait un test de dépis~ hvn_lbl_ [0] non
## [1] oui
## [9] manquantPréparation des données avec data.table
Tout d’abord, regardons sous quel format elles sont stockées.
class(menages)## [1] "tbl_df" "tbl" "data.frame"describe(menages)## [1814 obs. x 5 variables] tbl_df tbl data.frame
##
## $id_menage: Identifiant du ménage
## numeric: 1 2 3 4 5 6 7 8 9 10 ...
## min: 1 - max: 1814 - NAs: 0 (0%) - 1814 unique values
##
## $taille: Taille du ménage (nombre de membres)
## numeric: 7 3 6 5 7 6 15 6 5 19 ...
## min: 1 - max: 31 - NAs: 0 (0%) - 30 unique values
##
## $sexe_chef: Sexe du chef de ménage
## labelled double: 2 1 1 1 1 2 2 2 1 1 ...
## min: 1 - max: 2 - NAs: 0 (0%) - 2 unique values
## 2 value labels: [1] homme [2] femme
##
## $structure: Structure démographique du ménage
## labelled double: 4 2 5 4 4 4 5 2 5 5 ...
## min: 1 - max: 5 - NAs: 0 (0%) - 5 unique values
## 6 value labels: [0] pas d'adulte [1] un adulte [2] deux adultes de sexe opposé [3] deux adultes de même sexe [4] trois adultes ou plus avec lien de parenté [5] adultes sans lien de parenté
##
## $richesse: Niveau de vie (quintiles)
## labelled double: 1 2 2 1 1 3 2 5 4 3 ...
## min: 1 - max: 5 - NAs: 0 (0%) - 5 unique values
## 5 value labels: [1] très pauvre [2] pauvre [3] moyen [4] riche [5] très richeLes tableaux de données sont au format tibble (c’est-à-dire sont de la classe tbl_df ) et les variables catégorielles sont du type haven_labelled. Ceformat correspond au format de données si on les avait importées depuis SPSS avec l’extension haven.
En premier lieu, il nous faut convertir les tableaux de données au format data.table , ce qui peut se faire avec la fonction setDT . Par ailleurs, nous allons également charger en mémoire l’extensionlabelled pour la gestion des vecteurs labellisés.
library(labelled)
library(data.table)
setDT(menages)
setDT(femmes)
setDT(enfants)En premier lieu, il nous faut calculer la durée d’observation des enfants, à savoir le temps passé entre la date de naissance et la date de passation de l’entretien . Pour récupérer des variables du fichier femmes dans le fichier enfants , nous allons procéder à une fusion de table Pour le calcul de la durée d’observation, nous allons utiliser le package lubridate . Nous effectuerons l’analyse en mois (puisque l’âge au décès est connuen mois). Dès lors, la durée d’observation sera calculée en mois.
enfants <- merge(
enfants,
femmes[, .(id_femme, date_entretien)],
by = "id_femme",
all.x = TRUE
)
# duree observation en mois
library(lubridate, quietly = TRUE)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:data.table':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionenfants[, duree_observation := time_length(interval(date_naissance, date_entretien), unit = "months")]Regardons maintenant comment les âges au décès ont été collectés.
freq(enfants$age_deces)## n % val%
## 0 62 3.9 43.7
## 1 5 0.3 3.5
## 2 6 0.4 4.2
## 3 5 0.3 3.5
## 4 4 0.3 2.8
## 5 5 0.3 3.5
## 6 5 0.3 3.5
## 7 4 0.3 2.8
## 8 3 0.2 2.1
## 9 6 0.4 4.2
## 10 1 0.1 0.7
## 11 5 0.3 3.5
## 12 10 0.6 7.0
## 13 3 0.2 2.1
## 14 1 0.1 0.7
## 16 1 0.1 0.7
## 19 2 0.1 1.4
## 21 1 0.1 0.7
## 24 10 0.6 7.0
## 36 2 0.1 1.4
## 48 1 0.1 0.7
## NaN 1442 91.0 NALes âges au décès sont ici exprimés en mois révolus. Les décès à un mois révolu correspondent à desdécès entre 1 et 2 mois exacts. Nous allons ici adopter la troisième approche en considérant que les décès se répartissent de manière uniforme au sein d’un même mois. Nous aurons donc recours à la fonction runif qui permets de générer des valeurs aléatoires entre 0 et 1 selon une distribustion uniforme.
enfants[, age_deces_impute := age_deces + runif(.N)]Pour définir notre objet de survie, il nous faudra deux variables. Une première, temporelle, indiquant la durée à laquelle survient l’évènement étudié (ici le décès) pour ceux ayant vécu l’évènement et la durée d’observation pour ceux n’ayant pas vécu l’évènement (censure à droite). Par ailleurs, une seconde variable indiquant si les individus ont vécu l’évènement (0 pour non, 1 pour oui). Or, ici, la variable survie est codée 0 pour les décès et 1 pour ceux ayant survécu. Pour plus de détails, voir l’aide de la fonction Surv .
enfants[, deces := 0]
enfants[survie == 0, deces := 1]
var_label(enfants$deces) <- "Est décédé ?"
val_labels(enfants$deces) <- c(non = 0, oui = 1)
enfants[, time := duree_observation]
enfants[deces == 1, time := age_deces_impute]Occupons-nous maintenant des variables explicatives que nous allons inclure dans l’analyse. Tout d’abord, ajoutons à la table enfants les variables nécessaires des tables femmes et menages . Notons qu’il nous faudra importer id_menage de la table femmes pour pouvoir fusionner ensuite la table enfants avec la table menages . Par ailleurs, pour éviter une confusion sur la variable date_naissance, nous renommons à la volée cette variable de la table femmes en date_naissance_mere.
enfants <- merge(
enfants,
femmes[, .(
id_femme, id_menage, milieu, educ,
date_naissance_mere = date_naissance, nb_enf_ideal
)],
by = "id_femme",
all.x = TRUE
)
enfants <- merge(
enfants,
menages[, .(id_menage, structure, richesse)],
by = "id_menage",
all.x = TRUE
)Les variables catégorielles sont pour l’heure sous formes de vecteurs labellisés. Or, dans un modèle, il est impératif de les convertir en facteurs pour qu’elles soient bien traitées comme des variables catégorielles (autrement elles seraient traitées comme des variables continues). On aura donc recours à la fonction to_factor de l’extensionlabelled.
enfants[, sexe := to_factor(sexe)]
enfants[, richesse := to_factor(richesse)]Regardons plus attentivement, la variable structure.
freq(enfants$structure)## n % val%
## [0] pas d'adulte 0 0.0 0.0
## [1] un adulte 45 2.8 2.8
## [2] deux adultes de sexe opposé 505 31.9 31.9
## [3] deux adultes de même sexe 71 4.5 4.5
## [4] trois adultes ou plus avec lien de parenté 764 48.2 48.2
## [5] adultes sans lien de parenté 199 12.6 12.6Tout d’abord, la modalité «pas d’adulte» n’est pas représentée dans l’échantillon. On aura donc recours à l’argument drop_unused_labels pour ne pas conserver cette modalité. Par ailleurs, nous considérons que la situation familiale à partir de laquelle nous voudrons comparer les autres dans notre modèle, donc celle qui doit être considérée comme la modalité de référence, est celle du ménage nucléaire. Cette modalité («deux adultes de sexe opposé») n’étant pas la première, nous aurons recours à la fonction relevel .
enfants[, structure := to_factor(structure, drop_unused_labels = TRUE)]
enfants[, structure := relevel(structure, "deux adultes de sexe opposé")]Regardons la variable educ.
freq(enfants$educ)## n % val%
## [0] aucun 1084 68.4 68.4
## [1] primaire 368 23.2 23.2
## [2] secondaire 115 7.3 7.3
## [3] supérieur 17 1.1 1.1La modalité «supérieur» est peu représentée dans notre échantillon. Nous allons la fusionner avec la modalité «secondaire»
enfants[, educ2 := educ]
enfants[educ == 3, educ2 := 2]
val_label(enfants$educ2, 2) <- "secondaire ou plus"
val_label(enfants$educ2, 3) <- NULL
enfants[, educ2 := to_factor(educ2)]
freq(enfants$educ2)## n % val%
## aucun 1084 68.4 68.4
## primaire 368 23.2 23.2
## secondaire ou plus 132 8.3 8.3Calculons maintenant l’âge de la mère à la naissance de l’enfant et découpons le en groupes d’âges
enfants[, age_mere_naissance := time_length(
interval(date_naissance_mere, date_naissance),
unit = "years"
)]
enfants$gpage_mere_naissance <- cut(
enfants$age_mere_naissance,
include.lowest = TRUE, right = FALSE,
breaks=c(13, 20, 30, 50)
)
levels(enfants$gpage_mere_naissance) <- c(
"19 ou moins", "20-29", "30 et plus"
)
enfants$gpage_mere_naissance <- relevel(enfants$gpage_mere_naissance, "20-29")
freq(enfants$gpage_mere_naissance)## n % val%
## 20-29 888 56.1 56.1
## 19 ou moins 257 16.2 16.2
## 30 et plus 439 27.7 27.7Reste à calculer si le rang de naissance de l’enfant est supérieur au nombre idéal d’enfants tel que défini par la mère. On aura recours à la fonction rank appliquée par groupe (ici calculé séparément pour chaque mère). L’argument ties.method permet d’indiquer comment gérer les égalités (ici les naissances multiples, e.g. les jumeaux). Comme nous voulons comparer le rang de l’enfant au nombre idéal d’enfants,nous allons retenir la méthode “max” pour obtenir, dans le cas présent, le nombre total d’enfants déjà nés1 Avant de calculer un rang, il est impératif de trier préalablement le tableau
setorder(enfants, id_femme, date_naissance)
enfants[, rang := rank(date_naissance, ties.method = "max"), by = id_femme]
enfants[, rang_apres_ideal := "non"]
# note: unclass() requis en raison d'un bug non corrigé dans haven empéchant de comparer haven_labelled_spss et integer
enfants[rang > unclass(nb_enf_ideal), rang_apres_ideal := "oui"]
enfants[, rang_apres_ideal := factor(rang_apres_ideal)]
enfants[, rang_apres_ideal := relevel(rang_apres_ideal, "non")]Préparation des données avec dplyr
Tout d’abord, regardons sous quel format elles sont stockées
data(fecondite)
class(menages)## [1] "tbl_df" "tbl" "data.frame"describe(menages)## [1814 obs. x 5 variables] tbl_df tbl data.frame
##
## $id_menage: Identifiant du ménage
## numeric: 1 2 3 4 5 6 7 8 9 10 ...
## min: 1 - max: 1814 - NAs: 0 (0%) - 1814 unique values
##
## $taille: Taille du ménage (nombre de membres)
## numeric: 7 3 6 5 7 6 15 6 5 19 ...
## min: 1 - max: 31 - NAs: 0 (0%) - 30 unique values
##
## $sexe_chef: Sexe du chef de ménage
## labelled double: 2 1 1 1 1 2 2 2 1 1 ...
## min: 1 - max: 2 - NAs: 0 (0%) - 2 unique values
## 2 value labels: [1] homme [2] femme
##
## $structure: Structure démographique du ménage
## labelled double: 4 2 5 4 4 4 5 2 5 5 ...
## min: 1 - max: 5 - NAs: 0 (0%) - 5 unique values
## 6 value labels: [0] pas d'adulte [1] un adulte [2] deux adultes de sexe opposé [3] deux adultes de même sexe [4] trois adultes ou plus avec lien de parenté [5] adultes sans lien de parenté
##
## $richesse: Niveau de vie (quintiles)
## labelled double: 1 2 2 1 1 3 2 5 4 3 ...
## min: 1 - max: 5 - NAs: 0 (0%) - 5 unique values
## 5 value labels: [1] très pauvre [2] pauvre [3] moyen [4] riche [5] très richeNous allons charger en mémoire l’extension labelled pour la gestion des vecteurs labellisés en plus de dplyr.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(labelled)Nous allons calculée la duré des observations.
enfants_original <- enfants
library(lubridate)
enfants <- enfants %>%
left_join(
femmes %>% select(id_femme, date_entretien),
by = "id_femme"
) %>%
copy_labels_from(enfants_original) %>%
mutate(duree_observation = time_length(
interval(date_naissance, date_entretien),
unit = "months"
))enfants$date_entretien[enfants$duree_observation < 0] <-
enfants$date_entretien[enfants$duree_observation < 0] %m+% months(1)
enfants <- enfants %>%
mutate(duree_observation = time_length(
interval(date_naissance, date_entretien),
unit = "months"
))Les âges au décès sont ici exprimés en mois révolus. Les décès à un mois révolu correspondent à des décès entre 1 et 2 mois exacts.
Nous allons ici adopter la troisième approche en considérant que les décès se répartissent de manière uniforme au sein d’un même mois. Nous aurons donc recours à la fonction runif qui permets de générer des valeurs aléatoires entre 0 et 1 selon une distribustion uniforme.
enfants <- enfants %>%
dplyr::mutate(age_deces_impute = age_deces + runif(n()))Pour définir notre objet de survie, il nous faudra deux variables. Une première, temporelle, indiquant la durée à laquelle survient l’évènement étudié (ici le décès) pour ceux ayant vécu l’évènement et la durée d’observation pour ceux n’ayant pas vécu l’évènement (censure à droite). Par ailleurs, une seconde variable indiquant si les individus ont vécu l’évènement (0 pour non, 1 pour oui). Or, ici, la variable survie est codée 0 pour les décès et 1 pour ceux ayant survécu. Pour plus de détails, voir l’aide de la fonction Surv
enfants <- enfants %>%
mutate(deces = if_else(survie == 0, 1, 0)) %>%
set_variable_labels(deces = "Est décédé ?") %>%
set_value_labels(deces = c(non = 0, oui = 1)) %>%
mutate(time = if_else(deces == 1, age_deces_impute, duree_observation))Occupons-nous maintenant des variables explicatives que nous allons inclure dans l’analyse. Tout d’abord, ajoutons à la table enfants les variables nécessaires des tables femmes et menages . Notons qu’il nous faudra importer id_menage de la table femmes pour pouvoir fusionner ensuite la table enfants avec la table menages . Par ailleurs, pour éviter une confusion sur la variable date_naissance, nous renommons à la volée cette variable de la table femmes en date_naissance_mere.
enfants <- enfants %>%
left_join(
select(femmes,
id_femme, id_menage, milieu, educ,
date_naissance_mere = date_naissance, nb_enf_ideal
),
by = "id_femme"
) %>%
left_join(
select(menages, id_menage, structure, richesse),
by = "id_menage"
) %>%
copy_labels_from(enfants_original) %>%
copy_labels_from(femmes) %>%
copy_labels_from(menages)Les variables catégorielles sont pour l’heure sous formes de vecteurs labellisés. Or, dans un modèle, il est impératif de les convertir en facteurs pour qu’elles soient bien traitées comme des variables catégorielles (autrement elles seraient traitées comme des variables continues). On aura donc recours à la fonction to_factor de l’extensionlabelled.
enfants <- enfants %>%
mutate(sexe = to_factor(sexe), richesse = to_factor(richesse))Regardons plus attentivement, la variable structure
freq(enfants$structure)## n % val%
## [0] pas d'adulte 0 0.0 0.0
## [1] un adulte 45 2.8 2.8
## [2] deux adultes de sexe opposé 505 31.9 31.9
## [3] deux adultes de même sexe 71 4.5 4.5
## [4] trois adultes ou plus avec lien de parenté 764 48.2 48.2
## [5] adultes sans lien de parenté 199 12.6 12.6enfants <- enfants %>%
mutate(structure = relevel(
to_factor(structure, drop_unused_labels = TRUE),
"deux adultes de sexe opposé"
))
#Regardons la variable educ.
freq(enfants$educ)## n % val%
## [0] aucun 1084 68.4 68.4
## [1] primaire 368 23.2 23.2
## [2] secondaire 115 7.3 7.3
## [3] supérieur 17 1.1 1.1La modalité «supérieur» est peu représentée dans notre échantillon. Nous allons la fusionner avec la modalité «secondaire»
enfants <- enfants %>%
mutate(educ2 = ifelse(educ == 3, 2, educ)) %>%
set_value_labels(educ2 = c(
aucun = 0,
primaire = 1,
"secondaire ou plus" = 2
)) %>%
mutate(educ2 = to_factor(educ2))
freq(enfants$educ2)## n % val%
## aucun 1084 68.4 68.4
## primaire 368 23.2 23.2
## secondaire ou plus 132 8.3 8.3calcul de l’âge de la mere à la maissance et presentons sous forme d’interval d’âge
enfants <- enfants %>%
mutate(
age_mere_naissance = time_length(
interval(date_naissance_mere, date_naissance),
unit = "years"
),
gpage_mere_naissance = cut(
age_mere_naissance,
include.lowest = TRUE, right = FALSE,
breaks=c(13, 20, 30, 50)
)
)
levels(enfants$gpage_mere_naissance) <- c(
"19 ou moins", "20-29", "30 et plus"
)
enfants$gpage_mere_naissance <- relevel(enfants$gpage_mere_naissance, "20-29")
freq(enfants$gpage_mere_naissance)## n % val%
## 20-29 888 56.1 56.1
## 19 ou moins 257 16.2 16.2
## 30 et plus 439 27.7 27.7Reste à calculer si le rang de naissance de l’enfant est supérieur au nombre idéal d’enfants tel que défini par la mère
enfants <- enfants %>%
arrange(id_femme, date_naissance) %>%
group_by(id_femme) %>%
mutate(
rang = rank(date_naissance, ties.method = "max"),
rang_apres_ideal = ifelse(rang > as.integer(nb_enf_ideal), "oui", "non"),
rang_apres_ideal = factor(rang_apres_ideal, levels = c("non", "oui"))
)Kaplan-Meier
La courbe de survie de Kaplan-Meier s’obtient avec la fonction survfit de l’extensionsurvival.
library(survival)
km_global <- survfit(Surv(time, deces) ~ 1, data = enfants)
km_global## Call: survfit(formula = Surv(time, deces) ~ 1, data = enfants)
##
## n events median 0.95LCL 0.95UCL
## [1,] 1584 142 NA NA NAPour la représenter, on pourra avoir recours à la fonction ggsurvplot de l’extension survminer .
library(survminer, quietly = TRUE)## Warning: package 'survminer' was built under R version 4.2.3## Warning: package 'ggpubr' was built under R version 4.2.3##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myelomaggsurvplot(km_global)
On peut facilement représenter à la place la courbe cumulée des évènements (l’inverse de la courbe de survie) et la table des effectifs en fonction du temps.
ggsurvplot(km_global, fun = "event", risk.table = TRUE, surv.scale = "percent")
Pour comparer deux groupes (ici les filles et les garçons), il suffit d’indiquer la variable de comparaison à survfit .
km_sexe <- survfit(Surv(time, deces) ~ sexe, data = enfants)
km_sexe## Call: survfit(formula = Surv(time, deces) ~ sexe, data = enfants)
##
## n events median 0.95LCL 0.95UCL
## sexe=masculin 762 94 NA NA NA
## sexe=féminin 822 48 NA NA NALa fonction survdiff permets de calculer le test du logrank afin de comparer des courbes de survie. La mortalité infanto-juvénile diffère-t-elle significativement selon le sexe de l’enfant?
survdiff(Surv(time, deces) ~ sexe, data = enfants)## Call:
## survdiff(formula = Surv(time, deces) ~ sexe, data = enfants)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sexe=masculin 762 94 66.3 11.6 21.8
## sexe=féminin 822 48 75.7 10.2 21.8
##
## Chisq= 21.8 on 1 degrees of freedom, p= 3e-06ggsurvplot(km_sexe, conf.int = TRUE, risk.table = TRUE, pval = TRUE, data = enfants)
Modèle de Cox
Un modèle de Cox se calcule aisément avec coxph {survival}.
mod1 <- coxph(
Surv(time, deces) ~ sexe + milieu + richesse +
structure + educ2 + gpage_mere_naissance + rang_apres_ideal,
data = enfants
)
mod1## Call:
## coxph(formula = Surv(time, deces) ~ sexe + milieu + richesse +
## structure + educ2 + gpage_mere_naissance + rang_apres_ideal,
## data = enfants)
##
## coef exp(coef)
## sexeféminin -0.810893 0.444461
## milieu 0.652254 1.919862
## richessepauvre -0.075035 0.927711
## richessemoyen 0.325651 1.384931
## richesseriche 0.348630 1.417124
## richessetrès riche 0.469238 1.598775
## structureun adulte -0.145170 0.864875
## structuredeux adultes de même sexe 0.608983 1.838561
## structuretrois adultes ou plus avec lien de parenté 0.052616 1.054025
## structureadultes sans lien de parenté -0.128580 0.879343
## educ2primaire -0.032828 0.967705
## educ2secondaire ou plus -0.206516 0.813413
## gpage_mere_naissance19 ou moins -0.316501 0.728694
## gpage_mere_naissance30 et plus -0.005125 0.994888
## rang_apres_idealoui 1.393110 4.027357
## se(coef) z p
## sexeféminin 0.177832 -4.560 5.12e-06
## milieu 0.269618 2.419 0.0156
## richessepauvre 0.250611 -0.299 0.7646
## richessemoyen 0.248014 1.313 0.1892
## richesseriche 0.298173 1.169 0.2423
## richessetrès riche 0.428549 1.095 0.2735
## structureun adulte 0.600295 -0.242 0.8089
## structuredeux adultes de même sexe 0.376593 1.617 0.1059
## structuretrois adultes ou plus avec lien de parenté 0.196633 0.268 0.7890
## structureadultes sans lien de parenté 0.305361 -0.421 0.6737
## educ2primaire 0.205727 -0.160 0.8732
## educ2secondaire ou plus 0.367273 -0.562 0.5739
## gpage_mere_naissance19 ou moins 0.268056 -1.181 0.2377
## gpage_mere_naissance30 et plus 0.191673 -0.027 0.9787
## rang_apres_idealoui 0.603099 2.310 0.0209
##
## Likelihood ratio test=38.27 on 15 df, p=0.0008224
## n= 1584, number of events= 142De nombreuses variables ne sont pas significatives. Voyons si nous pouvons, avec la fonction step ,améliorer notre modèle par minimisation de l’AIC ou Akaike Information Criterion.
mod2 <- step(mod1)## Start: AIC=2026.67
## Surv(time, deces) ~ sexe + milieu + richesse + structure + educ2 +
## gpage_mere_naissance + rang_apres_ideal
##
## Df AIC
## - structure 4 2021.6
## - richesse 4 2022.3
## - educ2 2 2023.0
## - gpage_mere_naissance 2 2024.2
## <none> 2026.7
## - rang_apres_ideal 1 2028.4
## - milieu 1 2030.8
## - sexe 1 2046.8
##
## Step: AIC=2021.63
## Surv(time, deces) ~ sexe + milieu + richesse + educ2 + gpage_mere_naissance +
## rang_apres_ideal
##
## Df AIC
## - richesse 4 2016.6
## - educ2 2 2017.8
## - gpage_mere_naissance 2 2019.1
## <none> 2021.6
## - rang_apres_ideal 1 2023.2
## - milieu 1 2025.2
## - sexe 1 2041.9
##
## Step: AIC=2016.61
## Surv(time, deces) ~ sexe + milieu + educ2 + gpage_mere_naissance +
## rang_apres_ideal
##
## Df AIC
## - educ2 2 2012.9
## - gpage_mere_naissance 2 2014.5
## <none> 2016.6
## - rang_apres_ideal 1 2017.8
## - milieu 1 2018.3
## - sexe 1 2037.2
##
## Step: AIC=2012.89
## Surv(time, deces) ~ sexe + milieu + gpage_mere_naissance + rang_apres_ideal
##
## Df AIC
## - gpage_mere_naissance 2 2010.8
## <none> 2012.9
## - rang_apres_ideal 1 2014.0
## - milieu 1 2015.3
## - sexe 1 2033.4
##
## Step: AIC=2010.83
## Surv(time, deces) ~ sexe + milieu + rang_apres_ideal
##
## Df AIC
## <none> 2010.8
## - rang_apres_ideal 1 2012.1
## - milieu 1 2013.5
## - sexe 1 2031.2On peut obtenir facilement les coefficients du modèle avec l’excellente fonction tidy de l’extension broom. Ne pas oublier de préciser exponentiate = TRUE . En effet, dans le cas d’un modèle de Cox, l’exponentiel des coefficients corresponds auratio des risques instantannésouhazard ratio (HR) en anglais.
library(broom, quietly = TRUE)
tidy(mod2, exponentiate = TRUE)## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 sexeféminin 0.443 0.177 -4.59 0.00000445
## 2 milieu 1.52 0.199 2.10 0.0362
## 3 rang_apres_idealoui 3.57 0.584 2.18 0.0294#Pour représenter ces rapports de risque, on peut ici encore avoir recours à la fonction ggcoef_model de l’extension GGally.Vérification de la validité du modèle
Un modèle de Cox n’est valable que sous l’hypothèse de la proportionnalité des risques relatifs. Selon cette hypothèse les résidus de Schoenfeld ne dépendent pas du temps. Cette hypothèse peut être testée avec la fonction cox.zph .
test <- cox.zph(mod2)
test## chisq df p
## sexe 0.3872 1 0.53
## milieu 0.0617 1 0.80
## rang_apres_ideal 0.7628 1 0.38
## GLOBAL 1.2360 3 0.74Une valeur de p inférieure à 5 % indique que l’hypothèse n’est pas vérifiée. Il apparaît que p est supérieur à 5 % globalement et pour chaque variable prise individuellement. Notre modèle est donc valide.
Il est possible de représenter la distribution des résidus de Schoenfeld à l’aide de ggcoxzph de l’extensionsurvminer, afin de voir si leur répartition change au cours du temps.
ggcoxzph(test)## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values