Add horizontal line to an "effects plot" with different intercepts per panel

Question

I'm using the effects package to plot the effects from an ordinal logistic regression and I'd like to add a different horizontal line (represented the expected baseline probabilities if the predictor has no effect) to the multi-panel plot with a different baseline probability for each panel of the plot.

Here's an example dataset:

# Simulate some data
data = data.frame(predictor_x = rnorm(100, mean = 5, sd = 1),
                    response_y = as.factor(rbinom(100, size = 2, prob = 0.5)))
# Make a model
library(MASS)
model = polr(response_y ~ predictor_x, data = data, Hess= TRUE)
# Plot
library(effects)
p1 = plot(Effects(focal.predictor = "predictor_x", model))
# View
p1
# Try to add a red dashed line for expected probabilities
p1 + abline(h = c(0.25, 0.5, 0.25), col = "red", lty = 3) #failed

Given the plots, I'd like to add a horizontal line at y = .25, for the upper plot, y = .5, for the middle plot, and y = .25, for the lower plot.

Also, I know the + is usually ggplot format so I'm also not sure if that will work either or if you can/how to add more information (like a line) to a saved plot object?

Answer 1

The plot method for Effect objects uses the lattice plotting system, which is neither base R nor ggplot based. However, the ggeffect package can effectively produce the same plot using the ggplot framework, which is a bit easier to work with. In your case, we can do:

library(MASS)
library(effects)
library(ggeffects)

data <- data.frame(predictor_x = rnorm(100, mean = 5, sd = 1),
                  response_y = as.factor(rbinom(100, size = 2, prob = 0.5)))

model <- polr(response_y ~ predictor_x, data = data, Hess= TRUE)

p1 <- plot(ggeffect(model))$predictor_x

p1$facet$params$nrow <- 3

p1 + 
  geom_hline(data = data.frame(response.level = c("X0", "X1", "X2"),
                               y = c(0.25, 0.5, 0.25)),
             aes(yintercept = y), color = "red3", linetype = 2) +
  geom_rug(data = within(data, response.level <- paste0("X", response_y)),
                         aes(x = predictor_x, y = 0.5), sides = "b")