In raredd/forest: Forest plots

library('forest')
library('knitr')

knitr::opts_chunk$set(
  fig.align = 'center', fig.path = 'inst/etc/', fig.width = 8
)

forest plots

this is an experimental package and not to be trusted

---

Installation

# install.packages('devtools')
devtools::install_github('raredd/forest')

Basic usage and supported models/objects

x <- c(
  '[survival](https://cran.r-project.org/web/packages/survival)' = 'coxph',
  '[coxphf](https://cran.r-project.org/web/packages/coxphf)' = 'coxphf',
  '[cmprsk](https://cran.r-project.org/web/packages/cmprsk)' = 'crr',
  '[cmprsk2](https://github.com/raredd/cmprsk2)' = 'crr2',
  stats = 'formula',
  stats = 'glm',
  '[logistf](https://cran.r-project.org/web/packages/logistf)' = 'logistf'
)

knitr::kable(
  data.frame(
    object = sprintf('<code>%s</code>', x),
    package = sprintf('**%s**', names(x))
  )
)

library('survival')
forest(coxph(Surv(time, status) ~ I(age / 10) + factor(sex) + factor(ph.ecog), lung))
box('outer')

More options

lung2 <- within(lung, {
  sex <- factor(sex, 1:2, c('Male', 'Female'))
  ph.ecog <- factor(ph.ecog)
})

lbl <- c('10-year increase', levels(lung2$sex), levels(lung2$ph.ecog))
forest(
  coxph(Surv(time, status) ~ I(age / 10) + sex + ph.ecog, lung2),
  header = c('Age at enrollment', 'Sex', 'ECOG PS'),
  plotArgs = list(
    cex = 2, layout = 'unified', show_conf = TRUE, xlim = c(0, 10),
    labels = lbl, reset_par = FALSE, names = c('', 'N (%)', 'HR (95% CI)', 'p-value')
  )
)
title(xlab = 'Hazard ratio (95% CI)')
box('outer')

Multiple models

models <- list(
  'Model 1' = coxph(Surv(time, status) ~ age + sex + ph.ecog, lung2),
  'Model 2' = coxph(Surv(time, status) ~ age + sex, lung2),
  'Model 3' = coxph(Surv(time, status) ~ age, lung2)
)

forest2(
  models, cex = 2, col.sig = c('grey70', 'green4'),
  panel_size = c(2, 3, 2), xlim = c(0, 10)
)
box('outer')

Panel options

set.seed(1)
x <- forest(
  coxph(Surv(time, status) ~ I(age / 10) + sex + ph.ecog, lung2),
  plot = FALSE
)
y <- replicate(7, rnorm(20), simplify = FALSE)

plot(
  x, cex = 2, panel_size = c(4, 3, 2), xlim = c(-5, 5),
  left_panel = list(
    'Mean y' = sapply(y, function(yy) sprintf('%.3f', mean(yy))),
    'Median y' = sapply(y, function(yy) sprintf('%.3f', median(yy)))
  ),
  center_panel = {
    panel_box(y)
    axis(1, pos = 0.5)
  }
)
box('outer')

library('rawr') ## for tplot
hr_ci <- x$cleanfp_list$numeric[1:3]
plot(
  x, xlim = c(0, 10), reset_par = FALSE,
  center_panel = {
    panel_tplot(
      rev(asplit(hr_ci, 1)), type = 'd', cex = 3,
      pch = c(16, 1, 1), group.pch = FALSE,
      col = c(1, 2, 2), group.col = FALSE
    )
    axis(1, pos = 0.5)
  }
)
legend(
  6, 8.5, legend = c('Point estimate', '95% CI'),
  col = 1:2, pch = c(16, 1), bty = 'n', xpd = NA
)
box('outer')

Coerce raw data to a forest plot

set.seed(1)
x <- c(NA, 1, 2, NA, 1, 2, 3)
dat <- data.frame(
  N = x * 10,
  x = x,
  lower = x - 1,
  upper = x + 1,
  p.value = replace(runif(7), is.na(x), NA)
)
dat

x <- as.forest(
  x = dat$x, lower = dat$lower, upper = dat$upper, p.value = dat$p.value,
  labels = ifelse(is.na(dat$x), 'header', paste0('   ', x)),
  N = dat$N, P = dat$N / 100
)
x

plot(x, show_conf = TRUE)
box('outer')

Session info

within.list(sessionInfo(), loadedOnly <- NULL)

raredd/forest documentation built on Feb. 19, 2024, 9:22 p.m.