2015-02-07 16 views
5

Mam następujący model:Jak wykreślić model hazardu Coxa z wypustami

coxph(Surv(fulength, mortality == 1) ~ pspline(predictor)) 

gdzie fulength jest czas trwania obserwacji (w tym śmiertelność), jest czynnikiem predykcyjnym śmiertelności.

Wyjście powyższego polecenia jest następująca:

      coef se(coef) se2 Chisq DF p  
pspline(predictor), line 0.174 0.0563 0.0562 9.52 1.00 0.002 
pspline(predictor), nonl      4.74 3.09 0.200 

Jak mogę wykreślić ten model tak, aby uzyskać piękny krzywego linii z 95% zespołów ufności i współczynnik ryzyka na osi y? Co mam nastawione na coś podobnego do tego:

enter image description here

+0

http : //stat.ethz.ch/R-manual/R-devel/library/graphics/html/curve.html http://www.r-bloggers.com/plotting-95-confidence-bands-in-r- 2/ – efrem

+0

Używam pakietów rms/Hmisc Franka Harrell'a, które prawdopodobnie są w stanie dostarczyć coś bardzo podobnego do tego, chociaż nie wiem o działce po prawej stronie. Jestem statystycznie obrażony, wykreślając rozkład normalny pod wynikami mężczyzn i kobiet. Nie wiem, czy RMS obsługuje psplines, ponieważ Frank woli ograniczone sześcienne splajny, ale jeśli zamieścisz jakieś dane, z chęcią spróbuję. –

+0

Dzięki BondedDust. jak zainstalowałeś te pakiety? Podczas próby instalacji otrzymuję komunikaty o błędach, takie jak ta instalacja pakietu "TH.data" miała niezerowy kod wyjścia – Oposum

Odpowiedz

4

To jest, gdy pojawi się po uruchomieniu pierwszego przykładu w CPH z RMS pakietu:

n <- 1000 
set.seed(731) 
age <- 50 + 12*rnorm(n) 
label(age) <- "Age" 
sex <- factor(sample(c('Male','Female'), n, 
       rep=TRUE, prob=c(.6, .4))) 
cens <- 15*runif(n) 
h <- .02*exp(.04*(age-50)+.8*(sex=='Female')) 
dt <- -log(runif(n))/h 
label(dt) <- 'Follow-up Time' 
e <- ifelse(dt <= cens,1,0) 
dt <- pmin(dt, cens) 
units(dt) <- "Year" 
dd <- datadist(age, sex) 
options(datadist='dd') 
S <- Surv(dt,e) 

f <- cph(S ~ rcs(age,4) + sex, x=TRUE, y=TRUE) 
cox.zph(f, "rank")    # tests of PH 
anova(f) 
plot(Predict(f, age, sex)) # plot age effect, 2 curves for 2 sexes 

enter image description here

Ponieważ kombinacja pakietów rms/Hmisc korzysta z wykresów sieciowych, adnotacja z marginalną funkcją gęstości wiekowej musiała zostać wykonana przy pomocy funkcji sieci. Z drugiej strony, jeśli chcesz zmienić wartość reakcji na względną zagrożenia można po prostu dodać argumentu dla „fun = Exp” do Przewidzieć zadzwonić i relable wykres, aby uzyskać:

png(); plot(Predict(f, age, sex, fun=exp), ylab="Relative Hazard");dev.off() 

enter image description here

+0

Udało mi się zainstalować 'Hmisc', ale nie' rms'. pojawia się błąd: 'BŁĄD: instalowanie obiektów Rd nie powiodło się dla pakietu 'quantreg' BŁĄD: zależność 'TH.data' nie jest dostępna dla pakietu 'multcomp' BŁĄD: zależności 'quantreg', 'multcomp' nie są dostępne dla pakietu 'rms'' – Oposum

+0

Wygląda na to, że używane repozytorium ma uszkodzoną lub brakującą wersję multcomp i quantreg. Po prostu próbowałem zainstalować binarną wersję Mac z multcomp z repozytorium Berkeley i nie miałem błędów. –

+0

Wykryliśmy alternatywny sposób instalacji brakujących pakietów. Przechodzę przez wszystkie kroki, ale kiedy próbuję wykreślić 'plot (Predict (f, predictor, fun = exp), ylab =" Relative Hazard ")' Otrzymuję błąd 'Error in value.chk (at, which (name == n), NA, np, lim): zmienny predykator nie ma limitów zdefiniowanych przez datadysta' – Oposum