¿Por qué son los p-valores, a menudo, superior en un riesgo proporcional de Cox, modelo que en la regresión logística?

Question

He estado aprendiendo acerca de los riesgos proporcionales de Cox modelo. Tengo un montón de experiencia en el ajuste de modelos de regresión logística, y así construir la intuición he estado comparando los modelos de ajuste utilizando coxph de la R de "supervivencia" con modelos de regresión logística ajuste utilizando glm con family="binomial".

Si ejecuto el código:

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

Puedo obtener los valores de p para la edad de 0.0419 y 0.0254 respectivamente. Del mismo modo, si puedo usar el sexo como un predictor, con o sin edad.

Esto me parece desconcertante, porque yo creo que el tomar la cantidad de tiempo transcurrido en cuenta a la hora de ajustar el modelo daría más poder estadístico que sólo el tratamiento de la muerte como un resultado binario, mientras que los valores de p parece coherente con uno que tiene menos potencia estadística. ¿Qué está pasando aquí?

Answer 1

El modelo de regresión logística se supone que la respuesta es un ensayo de Bernoulli (o más generalmente de un binomio, pero por simplicidad, vamos a mantenerlo 0-1). Un modelo de supervivencia supone que la respuesta es por lo general una vez al evento (de nuevo, hay generalizaciones de este que vamos a saltar). Otra manera de explicarlo es que las unidades que están pasando a través de una serie de valores hasta que se produce un evento. No es que una moneda es en realidad discretamente volcó en cada punto. (Que podría suceder, por supuesto, pero entonces se necesitaría un modelo de medidas repetidas—tal vez un GLMM.)

Su modelo de regresión logística que lleva a cada uno a la muerte como un tirón de la moneda ocurrida en el que la edad y llegó hasta la cola. Asimismo, considera que cada dato censurado como un solo tirón de la moneda que se produjo en la edad especificada y llegó hasta los jefes. El problema aquí es que es incoherente con lo que los datos que realmente son.

Aquí están algunas parcelas de los datos, y la salida de los modelos. (Nota que le doy la vuelta a las predicciones del modelo de regresión logística para la predicción de estar vivo, de modo que la línea que coincide con el condicional de la densidad de la trama.)

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

---

No estoy seguro de lo que todavía podría necesitar, pero puede ser útil considerar una situación en la que los datos fueron adecuados para el análisis de la supervivencia o de una regresión logística. Imagina un estudio para determinar la probabilidad de que un paciente va a ser readmitido en el hospital dentro de los 30 días de alta en un nuevo protocolo o estándar de cuidado. Sin embargo, todos los pacientes son seguidos a la readmisión, y no hay censura (esto no es muy realista), por lo que la hora exacta a la readmisión podría ser analizados con el análisis de supervivencia (viz., de riesgos proporcionales de Cox modelo aquí). Para simular esta situación, voy a usar las distribuciones exponenciales con tasas de .5 y 1, y utilizar el valor 1 como un corte para representar a los 30 días:

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

En este caso, vemos que el p-valor en el modelo de regresión logística (0.163) fue mayor que el valor de p a partir de un análisis de supervivencia (0.005). Para explorar esta idea, podemos extender la simulación para estimar el poder de un análisis de regresión logística frente a un análisis de supervivencia, y la probabilidad de que el p-valor en el modelo de Cox será menor que el valor de p a partir de la regresión logística. Yo también voy a utilizar 1.4 como el umbral, de manera que no me desventaja de la regresión logística mediante el uso de un óptimo de corte:

xs = seq(.1,5,.1)
xs[which.min(pexp(xs,.5)-pexp(xs,1))]  # 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

De modo que la potencia de la regresión logística es menor (alrededor de 75%) que para el análisis de supervivencia (93%), y el 90% de los p-valores de los análisis de supervivencia fueron inferiores a los correspondientes p-valores de la regresión logística. Tomando los tiempos de retardo en cuenta, en lugar de simplemente menor o mayor que un cierto umbral no ceder más poder estadístico como había intuido.