Last chance! 50% off unlimited learning
Sale ends in
This function is a more specialized version of the termplot() function. It
creates a plot with the spline against hazard ratio. The plot can additianally have
indicator of variable density and have multiple lines.
plotHR(
models,
term = 1,
se = TRUE,
cntrst = ifelse(inherits(models, "rms") || inherits(models[[1]], "rms"), TRUE, FALSE),
polygon_ci = TRUE,
rug = "density",
xlab = "",
ylab = "Hazard Ratio",
main = NULL,
xlim = NULL,
ylim = NULL,
col.term = "#08519C",
col.se = "#DEEBF7",
col.dens = grey(0.9),
lwd.term = 3,
lty.term = 1,
lwd.se = lwd.term,
lty.se = lty.term,
x.ticks = NULL,
y.ticks = NULL,
ylog = TRUE,
cex = 1,
y_axis_side = 2,
plot.bty = "n",
axes = TRUE,
alpha = 0.05,
...
)# S3 method for plotHR
print(x, ...)
# S3 method for plotHR
plot(x, y, ...)
The function does not return anything
A single model or a list() with several models
The term of interest. Can be either the name or the number of the covariate in the model.
Boolean if you want the confidence intervals or not
By contrasting values you can have the median as a reference point making it easier to compare hazard ratios.
If you want a polygon as indicator for your confidence interval. This can also be in the form of a vector if you have several models. Sometimes you only want one model to have a polygon and the rest to be dotted lines. This gives the reader an indication of which model is important.
The rug is the density of the population along the spline variable. Often this is displayed as a jitter with bars that are thicker & more common when there are more observations in that area or a smooth density plot that looks like a mountain. Use "density" for the mountain view and "ticks" for the jitter format.
The label of the x-axis
The label of the y-axis
The main title of the plot
A vector with 2 elements containing the upper & the lower bound of the x-axis
A vector with 2 elements containing the upper & the lower bound of the y-axis
The color of the estimate line. If multiple lines you can have different colors by giving a vector.
The color of the confidence interval. If multiple lines you can have different colors by giving a vector.
The color of the density plot. Ignored if you're using jitter
The width of the estimated line. If you have more than one model then provide the function with a vector if you want to have different lines for different width for each model.
The typeof the estimated line, see lty. If you have more than one model then provide the function with a vector if you want to have different line types for for each model.
The line width of your confidence interval. This is ignored if you're using polygons for all the confidence intervals.
The line type of your confidence interval. This is ignored if you're using polygons for all the confidence intervals.
The ticks for the x-axis if you desire other than the default.
The ticks for the y-axis if you desire other than the default.
Show a logarithmic y-axis. Not having a logarithmic axis might seem easier to understand but it's actually not really a good idea. The distance between HR 0.5 and 2.0 should be the same. This will only show on a logarithmic scale and therefore it is strongly recommended to use the logarithmic scale.
Increase if you want larger font size in the graph.
The side that the y axis is to be plotted, see axis() for details
Type of box that you want. See the bty description in graphical parameters (par). If bty is one of "o" (the default), "l", "7", "c", "u", or "]" the resulting box resembles the corresponding upper case letter. A value of "n" suppresses the box.
A boolean that is used to identify if axes are to be plotted
The alpha level for the confidence intervals
Any additional values that are to be sent to the plot() function
Sent the `plotHR` object to plot
Ignored in plot
The function allows for plotting multiple splines in one graph. Sometimes you might want to show more than one spline for the same variable. This allows you to create that comparison.
Examples of a situation where I've used multiple splines in one plot is when I want to look at a variables behavior in different time periods. This is another way of looking at the proportional hazards assumption. The Schoenfeld residuals can be a little tricky to look at when you have the splines.
Another example of when I've used this is when I've wanted to plot adjusted and unadjusted splines. This can very nicely demonstrate which of the variable span is mostly confounded. For instance - younger persons may exhibit a higher risk for a procedure but when you put in your covariates you find that the increased hazard changes back to the basic
Reinhard Seifert, Max Gordon
org_par <- par(xaxs = "i", ask = TRUE)
library(survival)
library(rms)
library(dplyr)
library(Gmisc)
# Get data for example
n <- 1000
set.seed(731)
ds <- tibble(age = round(50 + 12 * rnorm(n), 1),
smoking = sample(c("Yes", "No"), n, rep = TRUE, prob = c(.2, .75)),
sex = sample(c("Male", "Female"), n, rep = TRUE, prob = c(.6, .4))) |>
# Build outcome
mutate(h = .02 * exp(.02 * (age - 50) + .1 *
((age - 50) / 10)^3 + .8 *
(sex == "Female") + 2 *
(smoking == "Yes")),
cens = 15 * runif(n),
dt = -log(runif(n)) / h,
e = if_else(dt <= cens, 1, 0),
dt = pmin(dt, cens),
# Add missing data to smoking
smoking = case_when(runif(n) < 0.05 ~ NA_character_,
TRUE ~ smoking)) |>
set_column_labels(age = "Age",
dt = "Follow-up time") |>
set_column_units(dt = "Year")
library(splines)
fit.coxph <- coxph(Surv(dt, e) ~ bs(age, 3) + sex + smoking, data = ds)
plotHR(fit.coxph, term = "age", plot.bty = "o", xlim = c(30, 70), xlab = "Age")
dd <- datadist(ds)
options(datadist = "dd")
fit.cph <- cph(Surv(dt, e) ~ rcs(age, 4) + sex + smoking, data = ds, x = TRUE, y = TRUE)
plotHR(fit.cph,
term = 1,
plot.bty = "L",
xlim = c(30, 70),
ylim = 2^c(-3, 3),
xlab = "Age"
)
plotHR(fit.cph,
term = "age",
plot.bty = "l",
xlim = c(30, 70),
ylog = FALSE,
rug = "ticks",
xlab = "Age"
)
unadjusted_fit <- cph(Surv(dt, e) ~ rcs(age, 4), data = ds, x = TRUE, y = TRUE)
plotHR(list(fit.cph, unadjusted_fit),
term = "age",
xlab = "Age",
polygon_ci = c(TRUE, FALSE),
col.term = c("#08519C", "#77777799"),
col.se = c("#DEEBF7BB", grey(0.6)),
lty.term = c(1, 2),
plot.bty = "l", xlim = c(30, 70)
)
par(org_par)
Run the code above in your browser using DataLab