1) How can I change the y-axis to "odds ratio", "probability of mortality", and "mortality rate" for "fit" in the following example?
2) How can I change the y-axis to "hazard ratio" for "fit2" in the following example?
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)
I am looking forward to your help!
Update 1 Following BondedDust's answer I worte the following:
# 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%
Since the overall probability to die is 44.5% the calculation of the predicted probability and the resulting plot seems to be correct to me as non-statistician, isn't it?
If you want the Odds ratio then you need to add a fun=-argument to transform to the odds ratio scale:
plot(Predict(fit,fun=exp), anova=an, pval=TRUE, ylab="Odds ratio")
I'm not sure I know what you mean by changing to the "probability of mortality", and "mortality rate" for "fit". The inverse logit function is exp(x)/(1+exp(x)) but in order to construct an estimate for events or rates from the coefficients, you would need to incorporate the intercept term. Perhaps you can extract something useful from the fitted values component of the fit object to satisfy the requirements of your homework problem.
> 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" ...
The same method as was used for transforming the coefficients for log-odds to odds ratio works for transforming log-hazard to Hazard ratio:
plot(Predict(fit2, fun=exp), anova=an, pval=TRUE, ylab="Hazard ratio")