Pourquoi les valeurs de p sont-elles souvent plus élevées dans un modèle de risque proportionnel de Cox que dans la régression logistique?

J'ai appris le modèle de risque proportionnel de Cox. J'ai beaucoup d'expérience dans l'ajustement de modèles de régression logistique, et donc pour construire l'intuition, j'ai comparé des modèles ajustés à l'aide coxphde la «survie» R avec des modèles de régression logistique ajustés à l'aide glmde family="binomial".

Si je lance le code:

library(survival)
s = Surv(time=lung$time, event=lung$status - 1)
summary(coxph(s ~ age, data=lung))
summary(glm(status-1 ~ age, data=lung, family="binomial"))

J'obtiens des valeurs p pour l'âge de 0,0419 et 0,0254 respectivement. De même si j'utilise le sexe comme prédicteur, avec ou sans âge.

Je trouve cela déroutant car je pense que la prise en compte du temps écoulé lors de l'ajustement du modèle donnerait plus de puissance statistique que de simplement traiter la mort comme un résultat binaire, tandis que les valeurs de p sembleraient cohérentes avec une puissance statistique inférieure. Qu'est-ce qui se passe ici?

Le modèle de régression logistique suppose que la réponse est un essai de Bernoulli (ou plus généralement un binôme, mais pour simplifier, nous le garderons 0-1). Un modèle de survie suppose que la réponse est généralement un temps avant l'événement (encore une fois, il y a des généralisations de cela que nous allons ignorer). Une autre façon de dire que les unités passent par une série de valeurs jusqu'à ce qu'un événement se produise. Ce n'est pas qu'une pièce soit en fait discrètement retournée à chaque point. (Cela pourrait arriver, bien sûr, mais vous auriez alors besoin d'un modèle de mesures répétées, peut-être un GLMM.)

Votre modèle de régression logistique prend chaque décès comme un retournement de pièce qui s'est produit à cet âge et qui s'est produit. De même, il considère chaque donnée censurée comme un seul coup de pièce de monnaie qui s'est produit à l'âge spécifié et qui a surgi. Le problème ici est que cela ne correspond pas à ce que sont réellement les données.

Voici quelques tracés des données et la sortie des modèles. (Notez que je bascule les prédictions du modèle de régression logistique pour prédire être vivant afin que la ligne corresponde au tracé de densité conditionnelle.)

library(survival)
data(lung)
s = with(lung, Surv(time=time, event=status-1))
summary(sm <- coxph(s~age, data=lung))
# Call:
# coxph(formula = s ~ age, data = lung)
# 
#   n= 228, number of events= 165 
# 
#         coef exp(coef) se(coef)     z Pr(>|z|)  
# age 0.018720  1.018897 0.009199 2.035   0.0419 *
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
#     exp(coef) exp(-coef) lower .95 upper .95
# age     1.019     0.9815     1.001     1.037
# 
# Concordance= 0.55  (se = 0.026 )
# Rsquare= 0.018   (max possible= 0.999 )
# Likelihood ratio test= 4.24  on 1 df,   p=0.03946
# Wald test            = 4.14  on 1 df,   p=0.04185
# Score (logrank) test = 4.15  on 1 df,   p=0.04154
lung$died = factor(ifelse(lung$status==2, "died", "alive"), levels=c("died","alive"))
summary(lrm <- glm(status-1~age, data=lung, family="binomial"))
# Call:
# glm(formula = status - 1 ~ age, family = "binomial", data = lung)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -1.8543  -1.3109   0.7169   0.8272   1.1097  
# 
# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)  
# (Intercept) -1.30949    1.01743  -1.287   0.1981  
# age          0.03677    0.01645   2.235   0.0254 *
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
#     Null deviance: 268.78  on 227  degrees of freedom
# Residual deviance: 263.71  on 226  degrees of freedom
# AIC: 267.71
# 
# Number of Fisher Scoring iterations: 4
windows()
  plot(survfit(s~1))
windows()
  par(mfrow=c(2,1))
  with(lung, spineplot(age, as.factor(status)))
  with(lung, cdplot(age, as.factor(status)))
  lines(40:80, 1-predict(lrm, newdata=data.frame(age=40:80), type="response"),
        col="red")

Il peut être utile de considérer une situation dans laquelle les données étaient appropriées pour une analyse de survie ou une régression logistique. Imaginez une étude pour déterminer la probabilité qu'un patient soit réadmis à l'hôpital dans les 30 jours suivant son congé en vertu d'un nouveau protocole ou d'une nouvelle norme de soins. Cependant, tous les patients sont suivis jusqu'à la réadmission, et il n'y a pas de censure (ce n'est pas terriblement réaliste), donc le temps exact de réadmission pourrait être analysé avec une analyse de survie (à savoir, un modèle de risques proportionnels de Cox ici). Pour simuler cette situation, je vais utiliser des distributions exponentielles avec des taux 0,5 et 1, et utiliser la valeur 1 comme seuil pour représenter 30 jours:

set.seed(0775)  # this makes the example exactly reproducible
t1 = rexp(50, rate=.5)
t2 = rexp(50, rate=1)
d  = data.frame(time=c(t1,t2), 
                group=rep(c("g1","g2"), each=50), 
                event=ifelse(c(t1,t2)<1, "yes", "no"))
windows()
  plot(with(d, survfit(Surv(time)~group)), col=1:2, mark.time=TRUE)
  legend("topright", legend=c("Group 1", "Group 2"), lty=1, col=1:2)
  abline(v=1, col="gray")

with(d, table(event, group))
#      group
# event g1 g2
#   no  29 22
#   yes 21 28
summary(glm(event~group, d, family=binomial))$coefficients
#               Estimate Std. Error   z value  Pr(>|z|)
# (Intercept) -0.3227734  0.2865341 -1.126475 0.2599647
# groupg2      0.5639354  0.4040676  1.395646 0.1628210
summary(coxph(Surv(time)~group, d))$coefficients
#              coef exp(coef)  se(coef)        z    Pr(>|z|)
# groupg2 0.5841386  1.793445 0.2093571 2.790154 0.005268299

Dans ce cas, nous voyons que la valeur de p du modèle de régression logistique ( 0.163) était supérieure à la valeur de p d'une analyse de survie ( 0.005). Pour approfondir cette idée, nous pouvons étendre la simulation pour estimer la puissance d'une analyse de régression logistique par rapport à une analyse de survie, et la probabilité que la valeur p du modèle de Cox soit inférieure à la valeur p de la régression logistique . Je vais également utiliser 1.4 comme seuil, afin de ne pas désavantager la régression logistique en utilisant un seuil sous-optimal:

xs = seq(.1,5,.1)
xs[which.max(pexp(xs,1)-pexp(xs,.5))]  # 1.4

set.seed(7458)
plr = vector(length=10000)
psv = vector(length=10000)
for(i in 1:10000){
  t1 = rexp(50, rate=.5)
  t2 = rexp(50, rate=1)
  d  = data.frame(time=c(t1,t2), group=rep(c("g1", "g2"), each=50), 
                  event=ifelse(c(t1,t2)<1.4, "yes", "no"))
  plr[i] = summary(glm(event~group, d, family=binomial))$coefficients[2,4]
  psv[i] = summary(coxph(Surv(time)~group, d))$coefficients[1,5]
}
## estimated power:
mean(plr<.05)  # [1] 0.753
mean(psv<.05)  # [1] 0.9253
## probability that p-value from survival analysis < logistic regression:
mean(psv<plr)  # [1] 0.8977

Ainsi, la puissance de la régression logistique est inférieure (environ 75%) à l'analyse de survie (environ 93%), et 90% des valeurs p de l'analyse de survie étaient inférieures aux valeurs p correspondantes de la régression logistique. La prise en compte des temps de latence, au lieu d'être juste inférieure ou supérieure à un certain seuil, donne plus de puissance statistique comme vous l'aviez prévu.