Mengapa nilai-p sering lebih tinggi dalam model bahaya proporsional Cox daripada dalam regresi logistik?

Saya telah belajar tentang model bahaya proporsional Cox. Saya memiliki banyak pengalaman pas model regresi logistik, dan model sehingga untuk membangun intuisi saya telah membandingkan cocok menggunakan coxphdari R "survival" dengan model regresi logistik cocok menggunakan glmdengan family="binomial".

Jika saya menjalankan kode:

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

Saya mendapatkan nilai-p untuk usia masing-masing 0,0419 dan 0,0254. Demikian pula jika saya menggunakan seks sebagai alat prediksi, dengan atau tanpa usia.

Saya menemukan ini membingungkan karena saya akan berpikir bahwa mengambil jumlah waktu berlalu ketika pas model akan memberikan lebih banyak kekuatan statistik daripada hanya memperlakukan kematian sebagai hasil biner, sedangkan nilai p akan tampak konsisten dengan yang memiliki lebih sedikit kekuatan statistik. Apa yang terjadi disini?

Model regresi logistik mengasumsikan responnya adalah uji coba Bernoulli (atau lebih umum binomial, tetapi untuk kesederhanaan, kami akan menyimpannya 0-1). Sebuah model survival mengasumsikan bahwa responnya biasanya waktu untuk peristiwa (sekali lagi, ada generalisasi dari ini yang akan kita lewati). Cara lain untuk mengatakannya adalah bahwa unit melewati serangkaian nilai sampai suatu peristiwa terjadi. Bukannya sebuah koin sebenarnya dibalik secara terpisah pada setiap titik. (Itu bisa terjadi, tentu saja, tetapi kemudian Anda akan membutuhkan model untuk tindakan berulang-ulang — mungkin GLMM.)

Model regresi logistik Anda mengambil setiap kematian sebagai flip koin yang terjadi pada usia itu dan muncul ekor. Demikian juga, ia menganggap setiap datum yang disensor sebagai flip koin tunggal yang terjadi pada usia yang ditentukan dan muncul kepala. Masalahnya di sini adalah bahwa itu tidak konsisten dengan apa sebenarnya data tersebut.

Berikut adalah beberapa plot data, dan output dari model. (Perhatikan bahwa saya membalik prediksi dari model regresi logistik menjadi prediksi yang hidup sehingga garis tersebut cocok dengan plot kepadatan bersyarat.)

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

Mungkin bermanfaat untuk mempertimbangkan situasi di mana data tersebut sesuai untuk analisis survival atau regresi logistik. Bayangkan sebuah penelitian untuk menentukan probabilitas bahwa seorang pasien akan diterima kembali ke rumah sakit dalam waktu 30 hari setelah dikeluarkan berdasarkan protokol atau standar perawatan baru. Namun, semua pasien diikuti untuk masuk kembali, dan tidak ada sensor (ini tidak terlalu realistis), sehingga waktu yang tepat untuk pendaftaran kembali dapat dianalisis dengan analisis kelangsungan hidup (yaitu, model bahaya proporsional Cox di sini). Untuk mensimulasikan situasi ini, saya akan menggunakan distribusi eksponensial dengan harga 0,5 dan 1, dan menggunakan nilai 1 sebagai cutoff untuk mewakili 30 hari:

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

Dalam hal ini, kita melihat bahwa p-value dari model regresi logistik ( 0.163) adalah lebih tinggi dari p-value dari analisis survival ( 0.005). Untuk mengeksplorasi ide ini lebih lanjut, kita dapat memperluas simulasi untuk memperkirakan kekuatan analisis regresi logistik vs analisis survival, dan probabilitas bahwa nilai-p dari model Cox akan lebih rendah daripada nilai-p dari regresi logistik . Saya juga akan menggunakan 1.4 sebagai ambang, sehingga saya tidak merugikan regresi logistik dengan menggunakan cutoff suboptimal:

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

Jadi kekuatan regresi logistik adalah rendah (sekitar 75%) dari analisis survival (sekitar 93%), dan 90% dari p-nilai dari analisis survival lebih rendah dari p-nilai yang sesuai dari regresi logistik. Memperhitungkan waktu jeda, alih-alih hanya kurang dari atau lebih besar dari beberapa ambang batas, menghasilkan lebih banyak kekuatan statistik seperti yang telah Anda intuisi.