Tutorial on Forest Plot

Let’s say I have my data ready and I want to fit a multivariable logistic regression model and make a forest plot of the resulting odds ratios. I’m going to walk through my process step by step.

Step 1: Fit the model

Using my data called wbns3 and my outcome variable FOOD_INSEC, I will fit a multivariable logistic regression model with the following predictors: annual household income (PPINCIMP), internet access (XACSNET), family size (FAMSIZE), education level (PPEDUCAT), and age (scaled by 10 years, PPAGE_10).

model = wbns3 %>% 
  glm(
    formula = FOOD_INSEC ~ PPINCIMP + XACSNET + FAMSIZE + PPEDUCAT + PPAGE_10, 
    family = binomial
    )

Step 2: Tidy the model output

My next step is to tidy the model output and prepare it for plotting. I will use the broom.helpers package to tidy the model output.

Important

Hover over the numbers on the left side of the code chunk to see the explanations for each step of the code.

library(broom.helpers)
MLR_tidy0 = 
  tidy_and_attach(
    model,
    conf.int=T,
    exponentiate = T
    ) %>%
  tidy_remove_intercept() %>%
  tidy_add_reference_rows() %>%
  tidy_add_estimate_to_reference_rows() %>%
  tidy_add_term_labels() %>%
  mutate(
    var_label = case_match(
      var_label,
      "PPINCIMP" ~ "Annual\nhousehold\nincome",
      "XACSNET"  ~ "Internet\naccess",               
      "FAMSIZE"  ~ "Family\nsize",                   
      "PPEDUCAT" ~ "Education\nlevel",                
      "PPAGE_10" ~ ""
    ),
    label = case_match(
      label, 
      "PPAGE_10" ~ "Age (by 10 years)", 
      .default = label
    )
  ) %>%
  mutate(
    label = fct_rev(
      fct_inorder(label)
      )
    )

1

Take model to tidy the model output, get confidence intervals, and exponentiate the estimates to get odds ratios.

2

Remove the intercept from the tidy output.

3

Add rows for the reference categories of the categorical predictors. For example, for internet access, the reference category is “No” and it will show on our forest plot. You can chose to remove this, but it often helpful in multi-level categorical variables to show the reference category on the plot.

4

Add the odds ratio of 1 to the reference rows so that they can be plotted on the forest plot.

5

Add variable labels so we can show the variable names for categorical variables.

6

We need to use mutate to change a few of the variable names and labels.

7

Update the variable labels to be more descriptive. For categorical variables, we want to name the variable. For continuous variables, we want this to be blank, so we include a ” “.

8

Add line breaks using \n between words. This will help the variable labels fit better on the plot. You may need to try out different placements of the line breaks to get the best fit for your plot.

9

Update labels to be more presentable and descriptive. This update should only be for continuous variables.

12

Reorder the levels of label so they show in the order they appear in the data frame. This is important for the forest plot to show the variables in the order we want.

Step 3: Make a forest plot

Using the tidied data on the model output, I can now make a forest plot of the odds ratios and confidence intervals.

library(ggplot2)

ggplot(
  data = MLR_tidy0,
  aes(
    y=label,
    x=estimate,
    xmin=conf.low,
    xmax=conf.high
    )
  ) +
  
  geom_point(
    size = 2
    ) +
  
  geom_errorbarh(
    width = 0.2         
    ) +
  
  geom_vline(
    xintercept=1,
    color='#C2352F',
    linetype='dashed',
    ) +
  
  facet_grid(
    rows = vars(var_label),
    scales = "free",
    space='free_y',
    switch = "y"
    ) +
  
  labs(
    x = "OR (95% CI)",
    title = "Odds ratios of food insecurity",
    y = element_blank()
    ) +
  
  theme_classic() +
  
  theme(
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    title = element_text(size = 12),
    strip.placement = "outside",
    strip.text.y.left = element_text(size = 12, angle = 0),
    strip.background = element_blank()
    )

1

Start a ggplot to specify the data and aesthetics.

2

Assign the data frame MLR_tidy0 to the plot.

3

Set up the aesthetics for the plot.

4

Set the y-axis to be the label variable, which contains the labels for each predictor. This will allow us to show the predictor names on the y-axis of the plot.

5

Set the x-axis to be the estimate variable, which contains the odds ratios.

6

Set the xmin aesthetic to be the conf.low variable, which contains the lower bound of the confidence intervals for the odds ratios.

7

Set the xmax aesthetic to be the conf.high variable, which contains the upper bound of the confidence intervals for the odds ratios.

8

Add points for the odds ratios. The size argument controls the size of the points.

9

Add horizontal error bars for the confidence intervals. The height argument controls the width of the error bars.

10

Add a vertical dashed line at x=1 to indicate the null value for odds ratios. This helps to visually assess which predictors are statistically significant (i.e., those whose confidence intervals do not cross 1).

11

Set the x-intercept of the vertical line to 1, which is the null value for odds ratios.

12

Set the color of the vertical line to a red color (you can choose any color you like).

13

Set the line type of the vertical line to dashed. You can choose other line types such as “dotted” or “solid” if you prefer.

14

Facet the plot by the variable labels so that each predictor is shown in its own row. This is helpful for categorical variables that need to be gathered by the variable.

15

Specify that the rows of the facets should be based on the var_label variable, which contains the descriptive labels for each predictor. This allows us to group the categorical levels of each predictor together in the plot. For example, all levels of annual household income will be grouped together under the label “Annual household income”.

16

The scales = "free" argument allows each facet to have its own x-axis scale, which can be helpful if the odds ratios vary widely across predictors.

17

The space='free_y' argument allows the height of each facet to adjust based on the number of levels in the predictor.

18

The switch = "y" argument moves the facet labels to the left side of the plot for better readability.

19

Add labels for the axes and title of the plot.

20

Set the x-axis label to “OR (95% CI)” to indicate that the x-axis shows odds ratios and their 95% confidence intervals.

21

Set the title of the plot to “Odds ratios of food insecurity” (you can choose any title you like).

22

Set the y-axis label to be blank since we are showing the predictor names on the y-axis and don’t need an additional label.

23

Use a classic theme for the plot to improve its appearance. I didn’t realize this before, but if you put it after the other theme(), it will override the other theme settings. So make sure to put theme_classic() before the other theme() to keep your custom theme settings.

24

Customize the theme of the plot to improve its appearance.

25

Set the size of the axis titles to 14. For example, the x-axis title “OR (95% CI)” will be size 14.

26

Set the size of the axis text to 14. For example, the x-axis tick labels (the numbers on the x-axis) will be size 14.

27

Set the size of the plot title to 14.

28

Move the facet strip labels to the outside of the plot for better readability. This puts the labels like “Annual household income” on the left side of the income labels.

29

Set the size of the facet strip text to 14 and set the angle to 0 so that the text is horizontal. This makes the facet labels easier to read.

30

Remove the background of the facet strips for a cleaner look. Try commenting this out! It might be nice to have blocks so it’s easier to see which levels belong to which predictor, especially if you have a lot of predictors with many levels.