Reputation: 59
I'm currently working on survival functions and especially on forest plots. I would like to plot the different levels of confidence intervals on the same line, with different colors or shades.
To do so, I use the forest function from the forestploter package. I use this package because it is greatly adaptable. However, I'm not able to have these different CIs. For the moment, I was able to get the different CIs but not on the same line.
Do anyone of you know how I could fix it ?
Here is an example of what I've been doing with the lung data from the survival package.
library(survival)
library(forestploter)
# Fit a Cox proportional hazards regression model
reg <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno, data = lung)
# Extract variable names, coefficients, and standard errors
variables <- names(coef(reg))
beta <- coef(reg)
se <- sqrt(diag(reg$var))
# Calculate hazard ratios (HR)
HR <- exp(beta)
# Calculate p-values
p <- 1 - pchisq((beta / se)^2, 1)
# Calculate confidence intervals (CI)
CI_95 <- exp(confint(reg))
CI_90 <- exp(confint(reg, level = 0.90))
# Format CI string for HR table
CI_hr_95 <- sprintf("%.2f (%.2f to %.2f)", HR, CI_95[, 1], CI_95[, 2])
# Create essai data.table with results
essai <- data.table(
variables, # Variable names
beta, # Coefficients
p = round(p, 3), # Rounded p-values
HR, # Hazard ratios
CI_95_lower = CI_95[, 1], # Lower CI for 95%
CI_95_upper = CI_95[, 2], # Upper CI for 95%
CI_90_lower = CI_90[, 1], # Lower CI for 90%
CI_90_upper = CI_90[, 2], # Upper CI for 90%
CI_hr_95 # Formatted HR CI string
)
essai$` `<- paste(rep(" ", 50), collapse = " ")
# Define forest plot theme
tm <- forest_theme(
base_size = 12, # Font size
base_family = "", # Font family (optional)
ci_pch = 10, # Point character for CI lines
ci_col = "black", # Color for CI lines
ci_alpha = 1, # Transparency of CI lines
ci_fill = NULL, # Fill color for CI regions (optional)
ci_lty = 1, # Line type for CI lines
ci_lwd = 2, # Line width for CI lines
ci_Theight = NULL, # Height for CI tick marks (optional)
legend_name = "CI threshold", # Legend label for CI lines
legend_value = c("90%", "95%") # Legend values for CI lines
)
# Prepare data for forest plot
forest_data <- essai[, c(1, 9, 10, 3)] # Select relevant columns
# Define forest plot elements
forest <- forest(
forest_data,
est = list(essai$HR, essai$HR), # Estimate values (HR for both lines)
lower = list(essai$CI_90_lower, essai$CI_95_lower), # Lower CI for each line
upper = list(essai$CI_90_upper, essai$CI_95_upper), # Upper CI for each line
sizes = 0.5, # Point size
ci_column = 3, # Column with CI information
ref_line = 1, # Reference line at HR=1
ticks_at = c(0, 1, 5.3), # X-axis ticks
xlim = c(0, 5.3), # X-axis limits
nudge_y = 0.3, # Vertical nudge for labels
# Even with nudge_y = 0, it is not on the same line
theme = tm, # Forest plot theme
arrow_lab = c("Longer WT", "Shorter WT") # Labels for HR lines (optional)
)
# Create the forest plot
plot(forest)
PS : For who might know this handbook, I would like to recreate the Figure 19.10 (p 468) from Regression Modeling Strategies: With Applications to Linear Models, Logistic Regression, and Survival Analysis (2015)
PS 2 : With my real data, using the forest function is great, because it makes easier to deal with label and to plot HR and CIs for interaction variables.
Thanks.
Upvotes: 0
Views: 91