Warum sind die p-Werte in einem Cox-Proportional-Hazard-Modell häufig höher als in einer logistischen Regression?

Ich habe etwas über das Cox-Proportional-Hazard-Modell gelernt. Ich habe viel Erfahrung mit der Anpassung von logistischen Regressionsmodellen, und um die Intuition zu verbessern, habe ich Modelle, die mit coxphdem R "Survival" passen, glmmit logistischen Regressionsmodellen verglichen, die mit dem R "Survival" passen family="binomial".

Wenn ich den Code ausführe:

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"))

Ich erhalte p-Werte für das Alter von 0,0419 bzw. 0,0254. Ebenso, wenn ich Sex als Prädiktor benutze, mit oder ohne Alter.

Ich finde das rätselhaft, weil ich denke, dass die Berücksichtigung der verstrichenen Zeit beim Anpassen des Modells mehr statistische Aussagekraft liefert als nur die Behandlung des Todes als binäres Ergebnis, während die p-Werte mit einem Wert übereinstimmen, der weniger statistische Aussagekraft hat. Was geht hier vor sich?

Das logistische Regressionsmodell geht davon aus, dass es sich bei der Antwort um einen Bernoulli-Versuch handelt (oder allgemeiner um einen Binomialversuch, aber der Einfachheit halber behalten wir ihn bei 0-1). Ein Überlebensmodell geht davon aus, dass die Antwort in der Regel eine Zeit bis zum Ereignis ist (es gibt wiederum Verallgemeinerungen davon, die wir überspringen werden). Eine andere Möglichkeit ist, dass Einheiten eine Reihe von Werten durchlaufen , bis ein Ereignis eintritt. Es ist nicht so, dass eine Münze an jedem Punkt diskret geworfen wird. (Das könnte natürlich passieren, aber dann brauchen Sie ein Modell für wiederholte Maßnahmen - vielleicht ein GLMM.)

Ihr logistisches Regressionsmodell nimmt jeden Todesfall als Münzwurf auf, der in diesem Alter stattgefunden hat und einen Endpunkt erreicht hat. Ebenso betrachtet es jedes zensierte Datum als einen einzelnen Münzwurf, der im angegebenen Alter stattfand und auftauchte. Das Problem hierbei ist, dass dies nicht mit den tatsächlichen Daten übereinstimmt.

Hier sind einige Diagramme der Daten und die Ausgabe der Modelle. (Beachten Sie, dass ich die Vorhersagen aus dem logistischen Regressionsmodell in die Vorhersage der Lebendigkeit umlege, sodass die Linie mit dem Diagramm der bedingten Dichte übereinstimmt.)

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")

Es kann hilfreich sein, eine Situation zu betrachten, in der die Daten für eine Überlebensanalyse oder eine logistische Regression geeignet waren. Stellen Sie sich eine Studie vor, um die Wahrscheinlichkeit zu bestimmen, dass ein Patient innerhalb von 30 Tagen nach seiner Entlassung nach einem neuen Protokoll oder nach einem neuen Versorgungsstandard wieder in das Krankenhaus eingeliefert wird. Alle Patienten werden jedoch bis zur Rückübernahme nachverfolgt, und es gibt keine Zensur (dies ist nicht besonders realistisch), sodass der genaue Zeitpunkt der Rückübernahme mithilfe einer Überlebensanalyse (hier ein Cox-proportionales Gefährdungsmodell) analysiert werden könnte. Um diese Situation zu simulieren, verwende ich Exponentialverteilungen mit Raten von 0,5 und 1 und verwende den Wert 1 als Grenzwert für 30 Tage:

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

In diesem Fall sehen wir , dass der p-Wert aus dem logistischen Regressionsmodell ( 0.163) war höher als der p-Wert aus einer Überlebensanalyse ( 0.005). Um diese Idee weiter zu untersuchen, können wir die Simulation erweitern, um die Leistung einer logistischen Regressionsanalyse im Vergleich zu einer Überlebensanalyse und die Wahrscheinlichkeit, dass der p-Wert aus dem Cox-Modell niedriger als der p-Wert aus der logistischen Regression ist, abzuschätzen . Ich werde auch 1,4 als Schwellenwert verwenden, damit ich die logistische Regression nicht durch einen suboptimalen Cutoff benachteilige:

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

Die Potenz der logistischen Regression ist also niedriger (ca. 75%) als die Überlebensanalyse (ca. 93%), und 90% der p-Werte aus der Überlebensanalyse waren niedriger als die entsprechenden p-Werte aus der logistischen Regression. Wenn Sie die Verzögerungszeiten berücksichtigen, erhalten Sie, anstatt nur einen Schwellenwert zu unterschreiten oder zu überschreiten, mehr statistische Leistung, als Sie gedacht hatten.