How to use plot_model to represent the coefficient of an interaction between a quadratic variable and a dummy in R?

Question

I am trying to plot an interaction between a squared variable and a dummy from a multi-level logistic regression in R. Whilst I have no problem plotting these interactions with a non-quadratic variable, I have no idea about how to do it with interactions.

I can plot perfectly these interactions:

plot_model(model, type="eff", 
            pred.type = "re", 
            terms = c("first_variable[all]","dummy_variable"),  show.data = F,ci.lvl=.95)

But whenever I try to plot this

plot_model(model, type="eff", 
            pred.type = "re", 
            terms = c("I(first_variable^2)[all]","dummy_variable"),  show.data = F,ci.lvl=.95)

I always get this response: Error: Some of the specified terms were not found in the model. Maybe misspelled?

I've tried spelling first_variable^2, as.is(first_variable^2), but nothing seems to work. And it is a variable contained in my regression: (glmer(y~first_variable+second_variable+I(first_variable^2)+I(first_variable^2):second_variable+(1|level2)+(1|level3), data = df3, family = binomial, nAGQ=1))

Any idea on how can I get around it?

Answer 1

I recently made a comment about my self-made function for plotting non-linear effects, but I just learned that sjPlot can do it just fine now (it couldn't at the time I made it).

Your code is not working when you use terms = c("I(first_variable^2)[all]")

Just use terms = c("first_variable[all]") and you should get the output you need. Unless you're trying to isolate the interaction with the quadratic term without incorporating the linear component?

In the example below, I get a plot of the interaction with a quadratic term:

hdp <- read.csv("https://stats.idre.ucla.edu/stat/data/hdp.csv")
m <- glmer(remission ~ Sex*RBC + Sex*I(RBC^2) +
             (1 | DID), data = hdp, family = binomial(link = "logit"))

sjPlot::plot_model(m, type="eff", pred.type = "re", terms = c("RBC[all]", "Sex"), show.data = F, ci.lvl = .95)

I also get the same output when I plot:

sjPlot::plot_model(m, type="pred", pred.type = "fe", terms = c("RBC[all]", "Sex"), show.data = F, ci.lvl = .95)