Dlaczego wartości p są często wyższe w modelu proporcjonalnego hazardu Coxa niż w regresji logistycznej?

Dowiedziałem się o modelu proporcjonalnego hazardu Coxa. Mam dużo doświadczenia okucia modele regresji logistycznej, a więc budować modele intuicji Byłem porównujące dopasowuje się stosując coxphod R „przetrwania” ze modele regresji logistycznej dopasowuje się stosując glmprzy family="binomial".

Jeśli uruchomię kod:

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

Dostaję wartości p dla wieku odpowiednio 0,0419 i 0,0254. Podobnie jeśli używam seksu jako predyktora, z wiekiem lub bez.

Uważam to za zagadkowe, ponieważ uważam, że uwzględnienie czasu, który upłynął przy dopasowaniu modelu, dałoby więcej mocy statystycznej niż tylko traktowanie śmierci jako wyniku binarnego, podczas gdy wartości p wydawałyby się spójne z tym, że ma on mniejszą moc statystyczną. Co tu się dzieje?

Model regresji logistycznej zakłada, że ​​odpowiedzią jest próba Bernoulliego (lub bardziej ogólnie dwumianowa, ale dla uproszczenia utrzymamy ją 0-1). Model przetrwania zakłada, że ​​reakcja jest zazwyczaj czasem na zdarzenie (ponownie, istnieją uogólnienia tego, które pomijamy). Innym sposobem jest to, że jednostki przechodzą przez szereg wartości, aż do wystąpienia zdarzenia. Nie jest tak, że moneta jest dyskretnie rzucana w każdym punkcie. ( Może się to oczywiście zdarzyć, ale wtedy potrzebujesz modelu do powtarzanych pomiarów - być może GLMM).

Twój model regresji logistycznej traktuje każdą śmierć jako rzut monetą, który miał miejsce w tym wieku i pojawił się ogon. Podobnie, uważa każdy cenzurowany układ odniesienia za pojedynczą monetę, która wystąpiła w określonym wieku i pojawiła się w głowach. Problem polega na tym, że jest to niezgodne z rzeczywistymi danymi.

Oto niektóre wykresy danych i dane wyjściowe modeli. (Zauważ, że przerzucam przewidywania z modelu regresji logistycznej na przewidywanie, że żyje, aby linia pasowała do wykresu gęstości warunkowej.)

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

Pomocne może być rozważenie sytuacji, w której dane były odpowiednie do analizy przeżycia lub regresji logistycznej. Wyobraź sobie badanie mające na celu określenie prawdopodobieństwa, że ​​pacjent zostanie ponownie przyjęty do szpitala w ciągu 30 dni od wypisu na podstawie nowego protokołu lub standardu opieki. Jednak wszyscy pacjenci są monitorowani o readmisję i nie ma cenzury (to nie jest strasznie realistyczne), więc dokładny czas na readmisję można przeanalizować za pomocą analizy przeżycia (tj. Tutaj model ryzyka proporcjonalnego Coxa). Aby zasymulować tę sytuację, użyję rozkładów wykładniczych ze współczynnikami .5 i 1 oraz użyję wartości 1 jako wartości granicznej reprezentującej 30 dni:

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

W tym przypadku widzimy, że wartość p z modelu regresji logistycznej ( 0.163) była wyższa niż wartość p z analizy przeżycia ( 0.005). Aby dalej badać ten pomysł, możemy rozszerzyć symulację, aby oszacować moc analizy regresji logistycznej w porównaniu z analizą przeżycia, oraz prawdopodobieństwo, że wartość p z modelu Coxa będzie niższa niż wartość p z regresji logistycznej . Użyję również wartości progowej 1,4, aby nie zaszkodzić regresji logistycznej poprzez użycie suboptymalnego odcięcia:

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

Zatem moc regresji logistycznej jest niższa (około 75%) niż analiza przeżycia (około 93%), a 90% wartości p z analizy przeżycia było niższe niż odpowiadające im wartości p z regresji logistycznej. Biorąc pod uwagę czasy opóźnienia, zamiast tylko mniejszego lub większego niż pewien próg, daje więcej mocy statystycznej, jak sobie wyobrażałeś.