Por que os valores de p geralmente são mais altos no modelo de risco proporcional de Cox do que na regressão logística?

Eu tenho aprendido sobre o modelo de risco proporcional de Cox. Eu tenho muita experiência montagem modelos de regressão logística e modelos para construir intuição Eu estive comparando ajuste usando coxpha partir do R "sobrevivência", com modelos de regressão logística ajuste usando glmcom family="binomial".

Se eu executar o 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"))

Recebo valores de p para idade de 0,0419 e 0,0254, respectivamente. Da mesma forma, se eu usar o sexo como preditor, com ou sem idade.

Acho isso intrigante porque acho que levar em consideração a quantidade de tempo decorrido ao ajustar o modelo daria mais poder estatístico do que apenas tratar a morte como um resultado binário, enquanto os valores de p pareceriam consistentes com um com menos poder estatístico. O que está acontecendo aqui?

O modelo de regressão logística assume que a resposta é um teste de Bernoulli (ou mais geralmente um binômio, mas por simplicidade, manteremos 0-1). Um modelo de sobrevivência assume que a resposta é tipicamente um momento para o evento (novamente, há generalizações disso que ignoraremos). Outra maneira de dizer isso é que as unidades estão passando por uma série de valores até que um evento ocorra. Não é que uma moeda seja lançada discretamente em cada ponto. (Isso poderia acontecer, é claro, mas você precisaria de um modelo para medidas repetidas - talvez um GLMM.)

Seu modelo de regressão logística considera cada morte como um lançamento de moeda que ocorreu nessa idade e surgiu como coroa. Da mesma forma, considera cada dado censurado como um único lançamento de moeda que ocorreu na idade especificada e surgiu. O problema aqui é que isso é inconsistente com o que os dados realmente são.

Aqui estão alguns gráficos dos dados e a saída dos modelos. (Observe que levanto as previsões do modelo de regressão logística para a previsão de vida, para que a linha corresponda ao gráfico de densidade condicional.)

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

Pode ser útil considerar uma situação em que os dados foram apropriados para uma análise de sobrevivência ou uma regressão logística. Imagine um estudo para determinar a probabilidade de um paciente ser readmitido no hospital dentro de 30 dias após a alta, sob um novo protocolo ou padrão de atendimento. No entanto, todos os pacientes são acompanhados até a readmissão e não há censura (isso não é muito realista), portanto, o tempo exato de readmissão pode ser analisado com a análise de sobrevida (a saber, um modelo de riscos proporcionais de Cox aqui). Para simular essa situação, usarei distribuições exponenciais com taxas 0,5 e 1 e usarei o valor 1 como ponto de corte para representar 30 dias:

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

Nesse caso, vemos que o valor de p do modelo de regressão logística ( 0.163) foi superior ao valor de p de uma análise de sobrevivência ( 0.005). Para explorar ainda mais essa idéia, podemos estender a simulação para estimar o poder de uma análise de regressão logística versus uma análise de sobrevivência e a probabilidade de que o valor p do modelo Cox seja menor que o valor p da regressão logística . Também usarei 1.4 como limite, para não prejudicar a regressão logística usando um ponto de corte abaixo do ideal:

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

Portanto, o poder da regressão logística é menor (cerca de 75%) do que a análise de sobrevivência (cerca de 93%) e 90% dos valores de p da análise de sobrevivência foram menores que os valores de p correspondentes da regressão logística. Levando em consideração os tempos de latência, em vez de apenas um valor menor ou maior que algum limite, gera mais poder estatístico como você havia intuído.