1
$\begingroup$

I am currently working with a MLR model comprising 1 numeric/continuous predictor variable (x1), several nominal categorical variables (x2 ... xi), and an interaction term between the continuous variable and one of the categorical variables (x1*x2). I wish to plot the relationship between y and x1 on a 2d scatterplot (ideally in ggplot2) with a line of best fit and confidence intervals. I believe that this is possible in theory, so long as the categorical variables are fixed at a pre-selected/reference level for plotting purposes. My understanding is that the model terms associated with the different levels of the categorical variables will just shift the intercept up/down the y-axis, and not fundamentally change the nature of the relationship between y and x1 (i.e., the slope). However, I have not been able to work out how to do generate such a plot in practice.

A reprex of a similar (toy) model is provided below:

# Generate data.frame
df<-
  data.frame(
    "y"=c(32,27,29,41,26,23,35,36,35,32,29,30,40,27,38,21,31,26,26,34,41,29,26,24),
    "x1"=c(28,32,36,40,44,48,52,56,60,64,68,72,72,68,64,60,56,52,48,44,40,36,32,28),
    "x2"=c("M","F"),
    "x3"=c("A","B","C"),
    "x4"=c("I","II","III","IV")
  )

df$x2<-
  as.factor(df$x2)

df$x3<-
  as.factor(df$x3)

df$x4<-
  as.factor(df$x4)

# Generate MLR model
lm<-
  lm(
    y~
      x1*x2+
      x3+
      x4,
    data=df
  )

summary(lm)

summary(lm) reads as follows: print-out of summary(lm)

Imagine that I wish to plot the relationship between y and x1 on a 2d scatterplot, and use x2=M, x3=A, and x4=III as the reference levels at which I wish to fix these covariates. How would I do so? I have tried manually calculating the predicted values for each of the data points if they were associated with these reference levels, and plotting them all, like so:

# Manually generate predictions
df$fixed<-
  32.424825+ # intercept
  df$x1*-0.003497+ # term for x1
  1*-4.525641+ # term for x2=M
  1*0+ # term for x3=A
  1*-3.666667+ # term for x4=III
  1*(df$x1*0.153846) # x1*x2 interaction term when x2=M

# Plot df$fixed vs df$x1
library(ggplot2)

p<-
  ggplot(
    data=df,
    mapping=aes(
      x=x1,
      y=fixed
      )
    )+
  geom_point()+
  geom_smooth(
    method="lm",
    se=T
  )

p

This approach has not worked for me, specifically as i) I wish to show the confidence intervals around the line of best fit, and ii) I have a very large dataset (~500k observations). In essence, I think that what I am doing here is passing ~500k fitted values to ggplot2, and then asking it to plot the line of best fit and associated confidence intervals. Unsurprisingly, there is essentially no uncertainty--given that the values are fitted, and that there are so many of them--so the confidence intervals are arbitrarily small. This would thus not appear to be the correct approach.

Is anyone aware of a method/package/function where I can plot y~x1 (ideally in ggplot2 graphics) with x2 .... xi held constant, in a way that still shows the uncertainty in the data (i.e., with confidence intervals)?

Thank you very much.

$\endgroup$
4
  • 1
    $\begingroup$ You might be able to do this with the ggeffects package. See its vignettes. $\endgroup$
    – dipetkov
    Commented Apr 7, 2023 at 16:12
  • 1
    $\begingroup$ @dipetkov thank you very much. I was able to use ggeffects to make these plots. Specifically, for anyone else who stumbles across this in future, I used ggpredict() and was able to plot using the ggpredict object directly through ggplot(), which proved much more powerful/flexible than the ggeffects::plot() function (which is nevertheless still helpful). What I actually wanted to plot was Marginal Effects--specifically, Marginal Effects at Representative Values. I did not know the terminology, hence why I struggled to find the correct tool for the job. $\endgroup$
    – PhelsumaFL
    Commented Apr 14, 2023 at 8:19
  • 1
    $\begingroup$ Thanks for the update. Aside: It's okay to answer one's own question and accept the answer; see here. Your question seemed to emphasize the plotting aspects; it wasn't clear that you have questions about the theory too. On that topic, ggeffects is a visualization interface to other packages that do the hard lifting of computing marginal effects such as emmeans. Admittedly, emmeans has a bit of a learning curve but its documentation and vignettes are very helpful. See here. $\endgroup$
    – dipetkov
    Commented Apr 14, 2023 at 8:31
  • $\begingroup$ @dipetkov thank you again. I think that I understand the theory behind the plots that I am generating, but wasn't clear on the terminology. Other people might be in a similar situation in future if I was. Hopefully this post will help them. I have answered the question myself as suggested. $\endgroup$
    – PhelsumaFL
    Commented Apr 14, 2023 at 8:43

1 Answer 1

1
$\begingroup$

I was able to generate these plots using ggeffects(), as suggested by dipetkov in the comments. An example is below. It is possible to alter the size of the confidence intervals with ci.lvl within the ggpredict() function, so I was able to create multiple CI bands simply by generating two objects of class("ggeffects")--one with CI set to .95 and one to .999--and plotting both with geom_ribbon(), using a lower alpha for plotting of the wider CI band.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.