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

9

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ż ?mi0

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:

cph () termplot2

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:

coxph () termplot2

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łę:

Działka kratowa

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łę:

Kontrastowa fabuła

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.

Max Gordon
źródło
3
Kilka komentarzy. Najpierw przeczytaj biostat.mc.vanderbilt.edu/Rrms, aby zobaczyć różnice między pakietami rms i Design. Po drugie, użyj plot () zamiast plot.Predict, aby zaoszczędzić pracę. Po trzecie, możesz łatwo wygenerować wykresy dla obu płci, np. Używając Predict (dopasowanie, wiek, płeć, zabawa = exp) # exp = anti-log; następnie wykreśl (wynik) lub wykreśl (wynik, ~ wiek | płeć). Nie używasz „x = NA” w Prognozie. rms używa grafiki kratowej, więc zwykłe parametry grafiki par i mfrow nie mają zastosowania. Zobacz przykłady w mojej ulotce z kursu rms na stronie biostat.mc . vanderbilt.edu/rms . Dla kontrastu.rms zapoznaj się z dokumentacją.
Frank Harrell,
1
Dziękuję bardzo za Twój wkład. Zaktualizowałem kod o lepsze przykłady i dodałem prof. Odpowiedź Thernau. PS Bardzo się cieszę, że planujesz nową wersję książki, poszerzając sekcję stronniczości punktu odcięcia, będzie bardzo przydatna jako odniesienie
Max Gordon
1
Możesz użyć ploti contrastzamiast plot.Predicti contrast.rms. Użyłbym bylub lengthwewnątrz seqzamiast timesi dałbym contrastdwie listy, żebyście dokładnie określili, co jest kontrastowane. Możesz także użyć cieniowania xYplotdla pasm pewności.
Frank Harrell,
1
Dzięki. Lubię używać fabuły. Przewiduję, bo wtedy dostaję odpowiednią pomoc w RStudio - coś, co w moim przypadku jest znacznie ważniejsze niż czas potrzebny na napisanie pełnej nazwy funkcji (używając autouzupełniania (tab)) tak naprawdę nie tracić tyle czasu).
Max Gordon,

Odpowiedzi:

5

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.

Frank Harrell
źródło
1
Dziękuję za Twoją odpowiedź. Czy uważasz, że powinienem wspomnieć o tym problemie prof. Terry Therneau, jeśli ma to być traktowane jako błąd / błędna interpretacja? Zajrzałem także do rozwiązań graficznych w pakiecie rms, nie do końca rozumiem użycie pliku kontrastu.rms w przypadku wykresów. Wydaje się, że Predict wykonuje termplot podobny wynik, ale nie mogę tego zrobić dokładnie tak, jak chcę ... zobacz moją aktualizację pytania.
Max Gordon,
2
Dobrze byłoby napisać do niego z zapytaniem i podziękować za przejazd na lotnisko, które dał mi kilka minut temu. Skomentuję powyżej pozostałe pytania.
Frank Harrell,