Trying to create a forest-plot in R using ggplot2 but running into some errors.

Data:

structure(list(explanatory_variables = structure(1:15, levels = c("Age", 
"Disease duration", "Smoking", "BMI", "HLA-B27", "Uveitis", "TNFi start year", 
"csDMARD", "CRP", "Psoriasis", "Arthritis", "NSAID", "Enthesitis", 
"TNFi type", "IBD"), class = "factor"), unadj_coef = c(0.93, 
0.92, 1.13, 1, 1.02, 1.32, 0.93, 0.83, 0.76, 1.42, 1.77, 0.04, 
1.26, 0.93, 1.3), adj_coef = c(0.91, 0.91, 1.13, 1, 1.02, 1.33, 
0.92, 0.83, 0.76, 1.42, 1.77, 0.04, 1.26, 0.93, 1.3), pct_change = c(-1.8, 
-0.9, 0.6, -0.5, 0.4, 0.3, -0.2, -0.2, 0.1, 0.1, 0.1, -0.1, -0.1, 
0, 0)), row.names = c(NA, -15L), class = "data.frame")

Code I have tried so far

# Sort data.frame by pct change
basdai_data <- basdai_data[order(abs(basdai_data$pct_change), decreasing = TRUE),]
basdai_data

# Reorder levels of explanatory variable
basdai_data$explanatory_variables  <- factor(basdai_data$explanatory_variables , levels = as.character(basdai_data$explanatory_variables ))
breaks <- as.character(levels(basdai_data$explanatory_variables))
breaks
basdai_data$explanatory_variables

# Create forest plot
ggplot(basdai_data, aes(x = pct_change, y = explanatory_variables))  
  geom_vline(xintercept = 0, linetype = "dashed")  
  scale_y_discrete(breaks = breaks, labels = breaks)  
  scale_x_continuous(breaks = seq(-3, 3, 1))  
  labs(x = "Percentage change", y = "Explanatory variable")  
  theme_pubr()  
  theme(text = element_text(size = 15, family = "Calibri"),
        axis.title = element_text(size = 20))

This is my error:

"Error in ans[ypos] <- rep(yes, length.out = len)[ypos] : replacement has length zero
In addition: Warning message:
In rep(yes, length.out = len) : 'x' is NULL so the result will be NULL"

My expectation was a forest-plot with on the left side of the y-axis the names of the explanatory variables, on the x-axis the percentage change, at 0 a dashed vertical line. However, instead I received the error above.

My questions are actually two-folded: how do I fix the error and as a bonus, would it be possible to add 4 columns of data ("explanatory_variable", "unadjusted coef.", "adjusted coef.", "pct. change" before the forest-plot? That would be perfect. The data is also present in the data-set.

CodePudding user response：

Your percentage change is not in the order you specified because you are sorting by the absolute value of percentage change. You can simply do:

basdai_data <- basdai_data[order(basdai_data$pct_change, decreasing = TRUE),]

You are getting the error simply because you have not included any data layers. If you add a geom_point, the result is:

ggplot(basdai_data, aes(x = pct_change, y = explanatory_variables))  
  geom_vline(xintercept = 0, linetype = "dashed")  
  geom_point()  
  scale_y_discrete(breaks = breaks, labels = breaks)  
  scale_x_continuous(breaks = seq(-3, 3, 1))  
  labs(x = "Percentage change", y = "Explanatory variable")  
  theme_pubr()  
  theme(text = element_text(size = 15, family = "Calibri"),
        axis.title = element_text(size = 20))

To get a table on the left side of the plot, you could do something like:

p1 <- ggplot(basdai_data, aes(x = pct_change, y = explanatory_variables))  
  geom_vline(xintercept = 0, linetype = "dashed")  
  geom_point()  
  scale_y_discrete(breaks = breaks, labels = breaks)  
  scale_x_continuous(breaks = seq(-3, 3, 1))  
  labs(x = "Percentage change", y = "Explanatory variable")  
  theme_pubr()  
  theme(text = element_text(size = 15, family = "Calibri"),
        axis.title = element_text(size = 20),
        axis.line.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

p2 <- ggplot(basdai_data, aes(x = 1, y = explanatory_variables))  
  geom_text(aes(label = explanatory_variables), hjust = 0)  
  geom_text(aes(x = 3, label = unadj_coef))  
  geom_text(aes(x = 4, label = adj_coef))  
  geom_text(aes(x = 5, label = pct_change))  
  scale_x_continuous(NULL, breaks = c(1, 3, 4, 5), limits = c(1, 5.5),
                     labels = c('Variable', 'Unadjusted coef', 'Adj. coef',
                                '% change'), position = 'top')  
  theme_void()  
  theme(axis.text.x.top = element_text(hjust = c(0, 0.5, 0.5, 0.5),
                                       face = 2))

library(patchwork)

p2   p1