Różny wykres prognostyczny od przeżycia coxph i rms cph

Stworzyłem własną, nieco ulepszoną wersję termplotu, której używam w tym przykładzie. Znajdziesz ją tutaj . Wcześniej pisałem na SO, ale im więcej o tym myślę, uważam, że to prawdopodobnie bardziej dotyczy interpretacji modelu proporcjonalnych zagrożeń Coxa niż faktycznego kodowania.

Problem

Kiedy patrzę na wykres współczynnika ryzyka, spodziewam się, że będę mieć punkt odniesienia, w którym przedział ufności wynosi naturalnie 0, i tak jest w przypadku, gdy używam cph () z, rms packageale nie kiedy używam coxph () z survival package. Czy poprawne zachowanie jest przez coxph (), a jeśli tak, to jaki jest punkt odniesienia? Również zmienna fikcyjna w coxph () ma interwał, a wartość jest inna niż ?$e^0$

Przykład

Oto mój kod testowy:

# Load libs
library(survival)
library(rms)

# Regular survival
survobj <- with(lung, Surv(time,status))

# Prepare the variables
lung$sex <- factor(lung$sex, levels=1:2, labels=c("Male", "Female"))
labels(lung$sex) <- "Sex"
labels(lung$age) <- "Age"

# The rms survival
ddist <- datadist(lung)
options(datadist="ddist")
rms_surv_fit <- cph(survobj~rcs(age, 4)+sex, data=lung, x=T, y=T)

Wykresy cph

Ten kod:

termplot2(rms_surv_fit, se=T, rug.type="density", rug=T, density.proportion=.05,
          se.type="polygon", yscale="exponential", log="y", 
          xlab=c("Age", "Sex"), 
          ylab=rep("Hazard Ratio", times=2),
          main=rep("cph() plot", times=2),
          col.se=rgb(.2,.2,1,.4), col.term="black")

daje ten wątek:

Wykresy Coxpha

Ten kod:

termplot2(surv_fit, se=T, rug.type="density", rug=T, density.proportion=.05, 
          se.type="polygon", yscale="exponential", log="y", 
          xlab=c("Age", "Sex"), 
          ylab=rep("Hazard Ratio", times=2),
          main=rep("coxph() plot", times=2),
          col.se=rgb(.2,.2,1,.4), col.term="black")

daje ten wątek:

Aktualizacja

Jak sugerował @Frank Harrell i po dostosowaniu się do sugestii w swoim ostatnim komentarzu otrzymałem:

p <- Predict(rms_surv_fit, age=seq(50, 70, times=20), 
             sex=c("Male", "Female"), fun=exp)
plot.Predict(p, ~ age | sex,
             col="black",
             col.fill=gray(seq(.8, .75, length=5)))

To dało bardzo ładną fabułę:

Ponownie spojrzałem na kontrast. Rms po komentarzu i wypróbowałem ten kod, który dał fabułę ... chociaż prawdopodobnie jest o wiele więcej do zrobienia :-)

w <- contrast.rms(rms_surv_fit, 
                  list(sex=c("Male", "Female"), 
                       age=seq(50, 70, times=20)))

xYplot(Cbind(Contrast, Lower, Upper) ~ age | sex, 
       data=w, method="bands")

Podaj tę fabułę:

AKTUALIZACJA 2

Prof. Thernau był na tyle uprzejmy, aby skomentować wątki o braku pewności siebie w talii:

Splajny wygładzające w Coxphie, podobnie jak w gamie, są znormalizowane, tak że suma (prognoza) = 0. Więc nie mam stałego pojedynczego punktu, dla którego wariancja jest bardzo mała.

Chociaż nie znam jeszcze GAM, wydaje się to odpowiadać na moje pytanie: wydaje się, że jest to kwestia interpretacji.

Myślę, że zdecydowanie powinien istnieć punkt, w którym przedział ufności ma zerową szerokość. Możesz także wypróbować trzeci sposób, polegający na użyciu wyłącznie funkcji rms. Pod plikiem pomocy dla pliku kontrast.rms znajduje się przykład uzyskania wykresu współczynnika ryzyka. Zaczyna się od komentarza # pokaż osobne szacunki według leczenia i płci. Musisz się zalogować, aby uzyskać współczynnik.