2015-07-21 83 views
0

1)如何在以下示例中将y轴更改为“比值比”,“死亡概率”和“死亡率”?如何绘制带有危险无线电,死亡概率或y轴死亡率的受限立方样条?

2)在以下示例中,如何将y轴更改为“fit2”的“风险比”?

library(Hmisc) 
library(survival) 
library(rms) 

data(pbc) 
d <- pbc 
rm(pbc) 
d$died <- ifelse(d$status == 2, 1, 0) 
d$status <- ifelse(d$status != 0, 1, 0) 

ddist <- datadist(d) 
options(datadist='ddist') 

fit <- lrm(status ~ rcs(age, 4), data=d) 
(an <- anova(fit)) 
plot(Predict(fit), anova=an, pval=TRUE) 

fit2 <- cph(Surv(time, status) ~ rcs(age, 4), data=d) 
(an2 <- anova(fit2)) 
plot(Predict(fit2), anova=an, pval=TRUE) 

我期待着您的帮助!

更新1 继BondedDust的答案,我又写道如下:

# probability 
getProbability <- function(x) { 
    exp(x)/(1+exp(x))*100 
} 

fit <- lrm(status ~ rcs(age, 4), data=d) 
(an <- anova(fit)) 
plot(Predict(fit, fun=getProbability), anova=an, pval=TRUE, ylab="Probability of death [%]") 

# overall probability to die 
table(d$status) 
round(table(d$status)[[2]]/sum(table(d$status))*100, digits=1) # = 44.5% 

由于总体概率死了,是44.5%的预测概率的计算所得的情节似乎是正确的,以我为非统计学家,不是吗?

回答

2

如果你想在赔率比率,那么你需要添加一个fun= -argument转换到的几率比规模:

plot(Predict(fit,fun=exp), anova=an, pval=TRUE, ylab="Odds ratio") 

我不知道我知道你changing to the "probability of mortality", and "mortality rate" for "fit"的意思。反逻辑函数为exp(x)/(1+exp(x)),但为了从系数中构建事件或速率的估计值,您需要合并截距项。也许你可以从拟合对象的拟合值组件中提取一些有用的东西来满足你的作业问题的要求。

> names(fit) 
[1] "freq"    "sumwty"   "stats"    "fail"    
[5] "coefficients"  "var"    "u"     "deviance"   
[9] "est"    "non.slopes"  "linear.predictors" "penalty.matrix" 
[13] "info.matrix"  "weights"   "call"    "Design"   
[17] "scale.pred"  "terms"    "assign"   "na.action"   
[21] "fail"    "interceptRef"  "nstrata"   
> str(fit$linear.predictors) 
Named num [1:418] -0.187 -0.235 0.41 -0.24 -0.538 ... 
- attr(*, "names")= chr [1:418] "1" "2" "3" "4" ... 

作为用于转化为对数比值的系数,以比值比同样的方法适用于转化数危害风险比:

plot(Predict(fit2, fun=exp), anova=an, pval=TRUE, ylab="Hazard ratio")