multiforest

This package is a modification (with a very limited scope) of the very nice package forestmodel that is available here: Github

multiforest will take a dataset and run uni- and multivariate poisson/binomial regressions on a binary outcome variable. These plots can take time to construct manually and one nice thing of using this package is that it makes it easy to do compare across sensitivity analyses (like how does results change when this group is excluded, or definition of outcome variable is changed)

Installation

devtools::install_github("LarsHernandez/multiforest")

Data

library(forestmodel)
library(multiforest)

library(tidyverse)
library(rlang)
library(labelled)
library(survival)
library(pec)
library(knitr)
library(kableExtra)

First we need some data, i take some random clinical trial data from the package pec and define 5 year survival. It’s important only to leave the data that will be contained in the models.

data(Pbc3, package = "pec")

pb <- Pbc3 %>%
  transmute(
    dead5yr = if_else(days < (365 * 5) & status == 2, T, F),
    tment   = if_else(tment == 0, "Placebo", "CyA"),    
    sex     = if_else(sex == 1, "Male", "Female"),
    age = case_when(
      between(age, 0, 50) ~ "18 - 50",
      between(age, 50, 60) ~ "50 - 60",
      between(age, 60, 140) ~ "60 - 75"),
    weight = case_when(
      between(weight, 0, 55) ~ "< 55kg",
      between(weight, 55, 65) ~ "55 - 65kg",
      between(weight, 65, 100) ~ "65 - 100kg"),
    stage = case_when(
      stage < 3 ~ "I-II",
      stage == 3 ~ "III",
      stage == 4 ~ "IV",
      is.na(stage) ~ "IV"),
    unit = case_when(
      unit == 1 ~ "Hvidovre",
      unit == 2 ~ "London",
      unit == 3 ~ "Copenhagen",
      unit == 4 ~ "Barcelona",
      unit == 5 ~ "Munich",
      unit == 6 ~ "Lyon"),
    gibleed = if_else(gibleed == 1, "Yes", "No")) %>%
  mutate(unit  = fct_relevel(unit, "London"),
         tment = fct_relevel(tment, "Placebo"))

kable(head(pb))

dead5yr	tment	sex	age	weight	stage	unit	gibleed
TRUE	CyA	Male	60 - 75	65 - 100kg	IV	London	No
FALSE	CyA	Female	18 - 50	55 - 65kg	III	London	No
FALSE	CyA	Female	18 - 50	65 - 100kg	III	London	No
TRUE	CyA	Female	60 - 75	\< 55kg	IV	London	Yes
TRUE	CyA	Female	60 - 75	\< 55kg	IV	London	No
TRUE	CyA	Male	60 - 75	55 - 65kg	IV	London	No

Before plotting we can make headers a bit nicer by applying a label from the package labelled to plot we use the function mforestmodel where we have to specify limits and the dependent variable.

var_label(pb) <- list(tment = "Treatment",
                      sex = "Sex",
                      age = "Age group",
                      weight = "Weight group",
                      stage = "Cancer stage",
                      unit = "Hospital",
                      gibleed = "Gastrointestinal bleeding")

mforestmodel(pb, dependent="dead5yr", lim=c(-2.4,2.4))

There are also a few options included for modifying the plot. For more options just use the raw function code from gíthub. The final plot is a ggplot object, and therefore normal options of ggplot can be supplied, like the labs argument.

mforestmodel(pb, dependent="dead5yr", lim=c(-2.4,2.4),
             pala="#2171b5", palb="#9ecae1", 
             #legend_position = c(0.55,0.93), 
             spaces = c(0.015,0.22,0.2,0.005,0.2,0.02),
             header = c("Group","Patients","RR of death (95% CI)","Univariate","P.val","Multivariate","P.val"),
             p_ci = 1, p_eps=0.01) + 
  labs(title="Relative risk of death within 5 years of diagnosis",
       subtitle="PBC3 - randomized clinical trial between 1 Jan. 1983 and 1 Jan. 1987")

Here is another example with some other data, and we can switch from poisson (RR) to a binomial/logistic model (OR)

data(cgd, package = "survival")

cg <- cgd %>% 
  select(treat, sex, hos.cat, status, weight, inherit, height) %>%
  mutate(
    hos.cat = as.character(hos.cat),
    bmi = weight / (height/100)^2,
    weight = case_when(
      between(weight, 0, 55) ~ "< 55kg",
      between(weight, 55, 65) ~ "55 - 65kg",
      between(weight, 65, 100) ~ "65 - 100kg"),
    bmi = case_when(
      between(bmi, 0, 17.4) ~ "< 17.5 kg/m\U00B2",
      between(bmi, 17.5, 23.9) ~ "17.5 - 24.9 kg/m\U00B2",
      between(bmi, 24, 220) ~ "> 25 kg/m\U00B2")) %>% 
  select(-weight, -height) %>% 
  mutate(bmi = fct_relevel(bmi, "< 17.5 kg/m\U00B2", "17.5 - 24.9 kg/m\U00B2"))
  
var_label(cg) <- list(treat = "Treatment",
                      sex = "Sex",
                      hos.cat = "Hospital region",
                      inherit = "Inherited",
                      bmi = "Body Mass Index")

mforestmodel(cg, dependent="status", lim=c(-3,2), est_family = "binomial")