Почему значения р часто выше в модели пропорционального риска Кокса, чем в логистической регрессии?

Я узнал о модели пропорционального риска Кокса. У меня большой опыт подбора моделей логистической регрессии, и поэтому для построения интуиции я сравнивал модели, подходящие для использования coxphиз «выживания» R, с моделями логистической регрессии, подходящими для использования glmс family="binomial".

Если я запускаю код:

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

Я получаю р-значения для возраста 0,0419 и 0,0254 соответственно. Точно так же, если я использую секс в качестве предиктора, с возрастом или без.

Я нахожу это озадачивающим, потому что я думаю, что принятие во внимание количества времени, потраченного на подбор модели, даст больше статистической силы, чем просто рассмотрение смерти как бинарного результата, в то время как значения p будут казаться соответствующими тому, который имеет меньшую статистическую мощность. Что здесь происходит?

Модель логистической регрессии предполагает, что ответом является испытание Бернулли (или, в более общем случае, биномиальное, но для простоты мы оставим его 0-1). Модель выживания предполагает, что ответом обычно является время на событие (опять же, есть некоторые обобщения, которые мы пропустим). Еще один способ показать, что единицы проходят через ряд значений, пока не произойдет событие. Дело не в том, что монета фактически дискретно подбрасывается в каждой точке. ( Конечно, это может произойти, но тогда вам понадобится модель для повторных измерений - возможно, GLMM.)

Ваша модель логистической регрессии воспринимает каждую смерть как бросок монеты, который произошел в этом возрасте и пришел к хвосту. Кроме того, он рассматривает каждую цензурированную датум как бросок монеты, который произошел в указанном возрасте и выпал на голову. Проблема здесь в том, что это не соответствует тому, что на самом деле представляют собой данные.

Вот несколько графиков данных и вывод моделей. (Обратите внимание, что я переворачиваю прогнозы из модели логистической регрессии в предсказание того, чтобы быть живым, чтобы линия соответствовала графику условной плотности.)

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

Может быть полезно рассмотреть ситуацию, в которой данные были подходящими для анализа выживания или логистической регрессии. Представьте себе исследование, чтобы определить вероятность повторного поступления пациента в больницу в течение 30 дней после выписки в соответствии с новым протоколом или стандартом медицинской помощи. Тем не менее, все пациенты следят за реадмиссией, и цензура отсутствует (это не очень реалистично), поэтому точное время до реадмиссии можно проанализировать с помощью анализа выживаемости (а именно, модели пропорциональных рисков Кокса здесь). Чтобы смоделировать эту ситуацию, я буду использовать экспоненциальные распределения со ставками .5 и 1 и использовать значение 1 в качестве предельного значения для представления 30 дней:

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

В этом случае мы видим, что значение p из модели логистической регрессии ( 0.163) было выше, чем значение p из анализа выживания ( 0.005). Для дальнейшего изучения этой идеи мы можем расширить моделирование, чтобы оценить силу анализа логистической регрессии в сравнении с анализом выживаемости, и вероятность того, что значение p в модели Кокса будет ниже, чем значение p в логистической регрессии , Я также буду использовать 1.4 в качестве порога, чтобы не ставить в невыгодное положение логистическую регрессию, используя субоптимальное ограничение:

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

Таким образом, мощность логистической регрессии является ниже (около 75%) , чем анализ выживаемости (около 93%), и 90% р-значений из анализа выживаемости были ниже , чем соответствующие р-значения из логистической регрессии. Принимая во внимание время запаздывания, вместо того, чтобы быть меньше или превышать некоторый порог, вы получаете большую статистическую мощность, как вы интуитивно поняли.