2015-02-07 86 views
5

一個Cox風險模型,我有以下模型:如何繪製與花鍵

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

其中fulength是隨訪(包括死亡)的持續時間,預測是死亡率的預測。

上面的命令的輸出是這樣的:

      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 

如何我可以繪製這個模型,使我得到95%的置信帶和危險比在y軸上的漂亮曲線玲瓏線?我所瞄準的是類似這樣的東西:

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 2015-02-07 18:24:31

+0

我使用Frank Harrell的rms/Hmisc軟件包,雖然我不知道右邊的情節,但可能能夠提供非常類似於輸出的內容。根據男性和女性的結果繪製正態分佈,我在統計上被冒犯了。我不知道rms是否支持psplines,因爲Frank喜歡三次樣條的限制,但是如果你發佈了一些數據,我很樂意嘗試一下。 – 2015-02-07 18:30:54

+0

感謝BondedDust。你是如何安裝這些軟件包的? 當我嘗試安裝時,出現如下錯誤消息:'TH.data'軟件包的安裝具有非零退出狀態 – Oposum 2015-02-07 18:35:54

回答

4

這是當你當你在RMS封裝的貼裝運行的第一個例子:

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

由於rms/Hmisc軟件包組合使用格點圖,具有邊際年齡密度特徵的註釋需要使用格點函數完成。在另一方面,如果你想改變的響應值相對危險,你可以再補充一個「有趣= EXP」的說法在預測打電話relable圖獲得:

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

enter image description here

+0

我能夠安裝'Hmisc'而不是'rms'。我得到的錯誤是這樣的: '錯誤:安裝Rd對象失敗的包'quantreg' 錯誤:依賴'TH.data'不可用於包'multcomp' 錯誤:依賴關係'quantreg','multcomp'不是可用於打包'rms'' – Oposum 2015-02-07 19:58:21

+0

聽起來就像您使用的存儲庫有multcomp和quantreg的損壞或丟失版本。我只是嘗試從Berkeley回購安裝二進制Mac版本的multcomp,並沒有錯誤。 – 2015-02-07 20:02:00

+0

找出安裝缺失軟件包的替代方法。我可以通過所有的步驟,但是當我嘗試繪製plot(Predict(f,predictor,fun = exp),ylab =「Relative Hazard」)時,我會得到錯誤值error.chk(at,which == n),NA,np,lim):變量預測器沒有由datadist定義的限制 – Oposum 2015-02-07 20:27:47