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 package
ale 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ż ?
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.
plot
icontrast
zamiastplot.Predict
icontrast.rms
. Użyłbymby
lublength
wewnątrzseq
zamiasttimes
i dałbymcontrast
dwie listy, żebyście dokładnie określili, co jest kontrastowane. Możesz także użyć cieniowaniaxYplot
dla pasm pewności.Odpowiedzi:
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.
źródło