predict.coxph()
calcule le rapport de risque par rapport à la moyenne de l'échantillon pour toutes les variables prédictives . Les facteurs sont convertis en prédicteurs factices comme d'habitude dont la moyenne peut être calculée. Rappelons que le modèle Cox PH est un modèle linéaire pour le log-hasard ln h ( t ) :plnh(t)
lnh(t)=lnh0(t)+β1X1+⋯+βpXp=lnh0(t)+Xβ
Où est le danger de base non spécifié. De manière équivalente, le danger h ( t ) est modélisé comme h ( t ) = h 0 ( t ) ⋅ e β 1 X 1 + ⋯ + β p X p = h 0 ( t ) ⋅ e X β . Le rapport de risque entre deux personnes i et i ' avec des valeurs prédictivesh0(t)h(t)h(t)=h0(t)⋅eβ1X1+⋯+βpXp=h0(t)⋅eXβii′XiXi′t
hi(t)hi′(t)=h0(t)⋅eXiβh0(t)⋅eXi′β=eXiβeXi′β
ii′b1,…,bpβ1,…,βpeXibeXi′b
À titre d'exemple dans R, j'utilise les données de l'annexe de John Fox sur le modèle Cox-PH qui fournit un très joli texte d'introduction. Tout d'abord, nous récupérons les données et construisons un modèle Cox-PH simple pour le délai d'arrestation des prisonniers libérés ( fin
: facteur - a reçu une aide financière avec codage fictif "no"
-> 0, "yes"
-> 1 age
,: âge au moment de la libération, prio
: nombre de condamnations antérieures):
> URL <- "http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt"
> Rossi <- read.table(URL, header=TRUE) # our data
> Rossi[1:3, c("week", "arrest", "fin", "age", "prio")] # looks like this
week arrest fin age prio
1 20 1 no 27 3
2 17 1 no 18 8
3 25 1 no 19 13
> library(survival) # for coxph()
> fitCPH <- coxph(Surv(week, arrest) ~ fin + age + prio, data=Rossi) # Cox-PH model
> (coefCPH <- coef(fitCPH)) # estimated coefficients
finyes age prio
-0.34695446 -0.06710533 0.09689320
Maintenant, nous connectons les moyennes de l'échantillon pour nos prédicteurs dans le eX b formule:
meanFin <- mean(as.numeric(Rossi$fin) - 1) # average of financial aid dummy
meanAge <- mean(Rossi$age) # average age
meanPrio <- mean(Rossi$prio) # average number of prior convictions
rMean <- exp(coefCPH["finyes"]*meanFin # e^Xb
+ coefCPH["age"] *meanAge
+ coefCPH["prio"] *meanPrio)
Maintenant, nous branchons les valeurs des prédicteurs des 4 premières personnes dans le eX b formule.
r1234 <- exp(coefCPH["finyes"]*(as.numeric(Rossi[1:4, "fin"])-1)
+ coefCPH["age"] *Rossi[1:4, "age"]
+ coefCPH["prio"] *Rossi[1:4, "prio"])
Maintenant, calculez le risque relatif pour les 4 premières personnes par rapport à la moyenne de l'échantillon et comparez-le à la sortie de predict.coxph()
.
> r1234 / rMean
[1] 1.0139038 3.0108488 4.5703176 0.7722002
> relRisk <- predict(fitCPH, Rossi, type="risk") # relative risk
> relRisk[1:4]
1 2 3 4
1.0139038 3.0108488 4.5703176 0.7722002
Si vous avez un modèle stratifié, la comparaison en predict.coxph()
est par rapport aux moyennes des strates, cela peut être contrôlé via l' reference
option expliquée dans la page d'aide.
meanFin <- mean(as.numeric(Rossi$fin) - 1)
n'a pas beaucoup de sens, carfin
est catégorique. N'en avez-vous pas besoinmodeFin <- get_Mode(Rossi$fin)
dans ce cas?fin
est binaire, donc la représentation numérique du facteur n'a que les valeurs 1 et 2. La soustraction de 1 nous donne la variable codée factice avec les valeurs 0 et 1 qui apparaît également dans la matrice de conception. Notez que cela ne fonctionnera pas pour les facteurs avec plus de 2 niveaux. Il est certainement discutable que la moyenne des variables fictives ait un sens, mais c'est ce quipredict.coxph()
fait.