I managaged to align the y-axis title using the margin argument instead of hjust. hjust did not work as expected, due to using the argument angle = 0. This is because justification works relative to the text direction.
To add the p-value labels I then used geom_text.
Plot produced:
Code used to produce the plot:
###################################################################################
# CREATE PLOTS-------------------------------------------------------------------#
###################################################################################
# Packages
library(tidyverse) # for maipulating data frames
# Amended forest plot function
forestplot_trend <- function (df, name = name, estimate = estimate, se = se, pvalue = NULL,
ptrend = NULL, colour = NULL, shape = NULL, logodds = FALSE,
psignif = 0.05, facet_var = NULL, ci = 0.95, ...)
{
# Ensure input is a data frame and logodds is a logical value
stopifnot(is.data.frame(df))
stopifnot(is.logical(logodds))
# Convert input column names to quosures for tidy evaluation
name <- enquo(name)
estimate <- enquo(estimate)
se <- enquo(se)
pvalue <- enquo(pvalue)
ptrend <- enquo(ptrend)
facet_var <- enquo(facet_var)
colour <- enquo(colour)
shape <- enquo(shape)
# Capture additional arguments passed to the function
args <- list(...)
# Compute the confidence interval constant (Z-score for normal distribution)
const <- stats::qnorm(1 - (1 - ci)/2)
# Prepare data frame: factorize `name` column, compute confidence intervals, and format labels
df <- df %>%
dplyr::mutate(
`:=`(!!name, factor(!!name, levels = !!name %>% unique() %>% rev(), ordered = TRUE)), # Reverse factor levels
.xmin = !!estimate - const * !!se, # Lower CI bound
.xmax = !!estimate + const * !!se, # Upper CI bound
.filled = TRUE, # Default filled state
.label = sprintf("%.2f", !!estimate) # Format estimate values as strings
)
# If logodds is TRUE, exponentiate the estimates and confidence intervals
if (logodds) {
df <- df %>% mutate(
.xmin = exp(.data$.xmin),
.xmax = exp(.data$.xmax),
`:=`(!!estimate, exp(!!estimate))
)
}
# If p-value is provided, mark significant values based on threshold (psignif)
if (!rlang::quo_is_null(pvalue)) {
df <- df %>% dplyr::mutate(.filled = !!pvalue < !!psignif)
}
# Initialize ggplot with estimate on x-axis and name on y-axis and ptrend values
g <- ggplot2::ggplot(df, aes(x = !!estimate, y = !!name)) +
geom_text(
aes(x = Inf, y = !!name, label = sprintf("%.3f", !!ptrend)),
hjust = 1.3, vjust = 0.5, family = "Arial"
)
# Apply a forest plot theme and add background elements
g <- g + theme_forest() + scale_colour_ng_d() + scale_fill_ng_d() +
geom_stripes() +
geom_vline(xintercept = ifelse(logodds, 1, 0), linetype = "solid", size = 0.4, colour = "black") # Reference line
# Adjust the x-axis to give more space for the labels
g <- g + coord_cartesian(xlim = c(min(df[[as_label(estimate)]]) - 0.5, max(df[[as_label(estimate)]]) + 0.5)) # Increase xlim range
# Add effect estimates with confidence intervals
g <- g + geom_effect(
ggplot2::aes(
xmin = .data$.xmin, xmax = .data$.xmax,
colour = !!colour, shape = !!shape, filled = .data$.filled
),
position = ggstance::position_dodgev(height = 0.5)
) +
ggplot2::scale_shape_manual(values = c(21L, 22L, 23L, 24L, 25L)) +
guides(colour = guide_legend(reverse = TRUE), shape = guide_legend(reverse = TRUE))
# Add optional title, subtitle, caption, x-label, and y-label if provided
if ("title" %in% names(args)) {
g <- g + labs(title = args$title)
}
if ("subtitle" %in% names(args)) {
g <- g + labs(subtitle = args$subtitle)
}
if ("caption" %in% names(args)) {
g <- g + labs(caption = args$caption)
}
if ("xlab" %in% names(args)) {
g <- g + labs(x = args$xlab)
}
if (!"ylab" %in% names(args)) {
args$ylab <- ""
}
g <- g + labs(y = args$ylab)
# Set axis limits if provided
if ("xlim" %in% names(args)) {
g <- g + coord_cartesian(xlim = args$xlim)
}
if ("ylim" %in% names(args)) {
g <- g + ylim(args$ylim)
}
# Adjust x-axis tick marks if provided and logodds is FALSE
if ("xtickbreaks" %in% names(args) & !logodds) {
g <- g + scale_x_continuous(breaks = args$xtickbreaks)
}
return(g)
}
# Create the forest plot function with added labels
forest_plot_data <- function(proteins, data, main_comparison) {
# Define the desired order for Carrier_comparison based on main_comparison
comparison_levels <- if (main_comparison == "C") {
c("CvsB", "AvsB")
} else {
c("AvsB", "CvsB")
}
# Filter data for selected proteins
data <- data %>%
filter(Assay %in% proteins)
# Ensure the Carrier_comparison is a factor with the desired level order
data$Carrier_comparison <- factor(data$Carrier_comparison, levels = comparison_levels)
# Plot
plot <- forestplot_trend(
df = data,
name = Assay, # Use combined label for y-axis
estimate = BETA,
se = SE,
pvalue = P_uncorrected,
ptrend = QM_P_trend, # Add QM_P_trend values
psignif = 0.01, # fill estimate if significant at this level
ci = 0.95,
logodds = FALSE,
xlab = "Beta (95% confidence intervals)",
ylab = "Analyte",
colour = Tertile,
shape = Tertile,
position = position_dodge(width = 0.7) # Increase dodge width for better separation
) +
facet_grid(. ~ Carrier_comparison, scales = "free", space = "free") +
theme(
text = element_text(family = "Arial"), # Change "Arial" to your preferred font
strip.text = element_text(size = 16, color = "black", hjust = 0.5, vjust = 5.5, margin = margin(t = 15, b = 15)),
axis.text.x = element_text(size = 12, color = "black"),
axis.text.y = element_text(size = 12, color = "black"),
axis.title.x = element_text(size = 14, color = "black", vjust = 0.5),
axis.title.y = element_text(size = 14, color = "black", face = "bold", angle = 0, vjust = 1.00, margin = margin(r = -50)), # Rotated y-axis title
legend.title = element_text(size = 14, color = "black"),
legend.text = element_text(size = 12, color = "black"), # Customize legend labels
legend.key.size = unit(1.5, "lines")
)
# Add "p-value" label above p-values in each panel
plot <- plot +
geom_text(
data = data.frame(Carrier_comparison = comparison_levels,
label = "bolditalic('p') * bold('-value')"), # Correct expression for bold and italic
aes(x = Inf, y = Inf, label = label),
hjust = 1.1, vjust = 1.5, size = 4, family = "Arial", color = "black",
parse = TRUE # Use parse = TRUE to interpret the label as an expression
)
return(plot)
}
# Create list of proteins
proteins <- levels(data$Assay)
# Generate plots
CvsB_forest_plot <- forest_plot_data(proteins, data, main_comparison = "C")
# Display the plot
CvsB_forest_plot
