---
title: "Customize Forest Plots and Tables"
output:
  rmarkdown::html_vignette:
    highlight: pygments
vignette: >
  %\VignetteIndexEntry{Customize Forest Plots and Tables}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8.5,
  fig.height = 5.5
)
```

```{r setup}
library(ggforestplotR)
library(ggplot2)
```

This article focuses on utility and enhanced customization of forest plots and accompanied tables.

## Group rows and control strip placement

`facet` creates section panels, and `facet_strip_position` controls which
side gets the strip labels.

```{r facet-right}
coefs <- data.frame(
  term = c("Age", "BMI", "Smoking", "Stage II", "Stage III"),
  estimate = c(0.12, -0.10, 0.18, 0.30, 0.46),
  conf.low = c(0.03, -0.18, 0.04, 0.10, 0.18),
  conf.high = c(0.21, 0.02, 0.32, 0.50, 0.74),
  sample_size = c(120, 115, 98, 87, 83),
  p_value = c(0.04, 0.15, 0.29, 0.001, 0.075),
  section = c("Clinical", "Clinical", "Clinical", "Tumor", "Tumor")
)


ggforestplot(
  coefs,
  facet = "section",
  facet_strip_position = "right",
  striped_rows = TRUE
)
```

## Subgroup analyses

### From a tibble

There are two ways to plot subgroup analyses using `ggforestplot()`. From a tibble, you can simply set `subgroup = "subgroup_variable_name"`. This functionality still works when some variables have subgroups and others don't, demonstrated here.

```{r mixed-subgroups}
mixed_coefs <- tibble::tribble(
  ~term,    ~subgroup, ~estimate, ~conf.low, ~conf.high,
  "Age",    NA_character_, 1.03,      1.01,       1.05,
  "White",  "Race",        1.01,      0.95,       1.07,
  "Black",  "Race",        0.89,      0.80,       0.99,
  "BMI",    NA_character_, 0.97,      0.94,       1.00,
  "Female", "Sex",         0.96,      0.89,       1.04,
  "Male",   "Sex",         0.98,      0.92,       1.06
)

ggforestplot(
  mixed_coefs,
  term = "term",
  subgroup = "subgroup",
  estimate = "estimate",
  conf.low = "conf.low",
  conf.high = "conf.high",
  exponentiate = TRUE,
  striped_rows = TRUE
) +
  add_forest_table()
```

### Continuous-by-categorical

The second way is to plot subgroup analyses from a fitted model. From your fitted model, call `tidy_forest_model()` and specify your `subgroup` and optionally your `focal` variable. `subgroup = "auto"` chooses the subgroup computations based on the model, otherwise you can specify the subgroup. `p_method` allows you to append overall p-values to the parent-row, or you can choose `p_method = "level"` if you want subgroup-specific p-values. Internally, the function calls on `marginaleffects` to compute interaction estimates.

```{r fitted-subgroups}
fit <- lm(wt ~ mpg*as.factor(cyl) + hp, data = mtcars)

fit |>
  tidy_forest_model(subgroup = "auto", focal = "mpg", p_method = "overall") |>
    ggforestplot(striped_rows = T) +
    add_forest_table(columns = c("term", "estimate", "p.value"))
```

### Categorical-by-categorical
```{r subgroup-categorical}
fit2 <- lm(wt ~ as.factor(gear)*as.factor(cyl) + hp + am, data = mtcars)

fit2 |>
  tidy_forest_model(subgroup = "gear", focal = "cyl", p_method = "level") |>
  ggforestplot(striped_rows = T) +
  theme(legend.position = "top") +
  add_forest_table(columns = c("term", "estimate", "p.value"))
```

## Distinct variable separation

Use `separate_groups` and `separate_lines` when you want a more distinct visual separation between variables. This is especially useful for categorical variables with many levels. `separate_groups` automatically appends the variable name to the level.

```{r separators}
block_coefs <- data.frame(
  term = c("race_black", "race_white", "race_other", "age", "bmi"),
  label = c("Black", "White", "Other", "Age", "BMI"),
  estimate = c(0.24, 0.08, -0.04, 0.12, -0.09),
  conf.low = c(0.10, -0.04, -0.18, 0.03, -0.17),
  conf.high = c(0.38, 0.20, 0.10, 0.21, -0.01),
  variable_block = c("Race", "Race", "Race", "Age", "BMI")
)

ggforestplot(
  block_coefs,
  label = "label",
  separate_groups = "variable_block",
  separate_lines = TRUE,
  striped_rows = TRUE
) +
  scale_y_discrete(limits = rev(c("BMI", "Age", "Race: White", 
                                  "Race: Black", "Race: Other")))
```

## Add a side table

`add_forest_table()` allows you to attach model information to the coefficient plot. The table can be added to either the left or right side and allows for some customization. You should **always** add the table **LAST**, after styling your plot because the function calls on `patchwork` internally. `patchwork` requires specific syntax to customize plots and is generally more difficult to get working correctly.

You can choose which columns from your dataframe to include in the table using the `columns` argument, and can change the labels using `column_labels`. If some of the term labels need to be changed, use `term_labels` to assign them new values. Some of the column labels are automatically assigned if no value is provided.

Notice how we are explicitly naming *p.value* column? This is necessary in most cases because aliases are not yet incorporated (but they will be...I promise I'm getting to it).

```{r left-side-table}
ggforestplot(
  coefs,
  facet = "section",
  facet_strip_position = "right",
  p.value = "p_value",
  striped_rows = TRUE,
  term_labels = c("Smoking" = "Smoking status")
) +
  add_forest_table(
    columns = c("term", "sample_size", "estimate", "p_value"),
    column_labels = c("term" = "Variable", "sample_size" = "N",
                      "estimate" = "Beta (95% CI)", "p_value" = "P-value")
  )
```

## Customize the table

`add_forest_table` also lets you change some minor styling elements of the forest table.

```{r}
ggforestplot(
  coefs,
  n = "sample_size",
  p.value = "p_value",
  striped_rows = TRUE
) +
  add_forest_table(
    position = "left",
    grid_lines = T,
    grid_line_linetype = 2,
    grid_line_colour = "red"
  )
```


## Split tables

`add_split_table()` can be used to create more traditional looking forest plots. You can choose which summary information goes to which side. Like `add_forest_table()`, it should be added after any plot-level styling. 

Use the `estimate_fmt` argument to change how your estimates are displayed. You can also control digits via `estimate_digits`, `interval_digits`, and `p_digits`.

```{r split-table}
ggforestplot(
  coefs,
  n = "sample_size",
  p.value = "p_value",
  striped_rows = TRUE
) +
  scale_x_continuous(limits = c(-.8,.8)) +
  add_split_table(
    left_columns = c("term","n"),
    right_columns = c("estimate","p"),
    column_labels = c("estimate" = "Beta [95% CI]"),
    estimate_fmt = "{estimate} [{conf.low}, {conf.high}]",
    estimate_digits = 2,
    interval_digits = 3,
    p_digits = 2
  ) 
```

## Plotting other types of model coefficients

You can use `exponentiate = TRUE` for models on the log-odds scale (or similar).

```{r logistic-regression-data}
data(CO2)

l1 <- glm(Treatment ~ conc + uptake + Type, family = binomial(link = "logit"), 
    data = CO2)
```

```{r logistic-regression, warning=FALSE}

ggforestplot(l1, exponentiate = TRUE, striped_rows = T, term_labels = c("TypeMississippi" = "Mississippi")) +
  add_forest_table(position = "left", 
                   columns = c("term", "estimate"))
```

We can do this for survival models as well.

```{r survival-analysis-data}
lung <- survival::lung

lung <- lung |>  
  dplyr::mutate(
    status = dplyr::recode(status, `1` = 0, `2` = 1)
  )

s1 <- survival::coxph(Surv(time, status) ~ sex + age + ph.karno + pat.karno, data = lung)
```

```{r survival-analysis-plot}
ggforestplot(s1, exponentiate = T, striped_rows = T) +
  add_forest_table()
```

## Compare multiple estimates

The `group` argument is handy when comparing estimates from several models.

```{r comparison}
comparison_coefs <- data.frame(
  term = rep(c("Age", "BMI", "Smoking", "Stage II", "Stage III"), 2),
  estimate = c(0.12, -0.10, 0.18, 0.30, 0.46, 0.08, -0.05, 0.24, 0.40, 0.58),
  conf.low = c(0.03, -0.18, 0.04, 0.10, 0.18, 0.00, -0.13, 0.10, 0.20, 0.30),
  conf.high = c(0.21, -0.02, 0.32, 0.50, 0.74, 0.16, 0.03, 0.38, 0.60, 0.86),
  model = rep(c("A", "B"), each = 5)
)

ggforestplot(
  comparison_coefs,
  group = "model",
  striped_rows = TRUE,
  dodge_width = 0.5
) +
  theme(legend.position = "top") +
  scale_color_manual(values = c("#1F968BFF", "#453781FF")) +
  labs(color = "Model") +
  add_forest_table(
    column_labels = c("term" = "Term", 
                      "model" = "Model", 
                      "estimate"  = "Estimate (95% CI)")
    )
```

You can also use `bind_forest_models()` to plot estimates from several fit models at once.

```{r bind-models, fig.width=8, fig.height=4.5, fig.dpi=300}
fit1 <- lm(mpg ~ cyl, data = mtcars)
fit2 <- lm(mpg ~ cyl + disp, data = mtcars)
fit3 <- lm(mpg ~ cyl + disp + wt, data = mtcars)

bound_models <- bind_forest_models(list(fit1,fit2,fit3), 
                                   model_labels = c("Unadjusted", 
                                                    "Adjusted", 
                                                    "Fully Adjusted"))

ggforestplot(bound_models, striped_rows = T, p.value = "p.value") +
  scale_x_continuous(limits = c(-6,1)) +
  theme(legend.position = "top") +
  scale_color_manual(values = c("#1F968BFF", "#453781FF", "#FDE725FF")) +
  add_forest_table(columns = c("term", "model","estimate", "p.value"),
                   p_digits = 4,
                   )
```