# Can sim_slopes handle polynomial models?

I am trying to compare slopes of a simple model with a polynomial effect in R:

> fit = lm(y ~ poly(x, 2)*w, data = df)


here x is a continuous variable, and w a categorical effect with 2 levels

I tried to do so using sim_slopes from the interactions package, but I get a subscript out of bound error :

> sim_slopes(fit, pred=x, modx=w, johnson_neyman = FALSE)

Error in summat[if (make.names(pred) != pred) bt(pred) else pred, ] : subscript out of bounds


fitting a model where x is not polynomial works

I am assuming this error is because poly(x, 2) generates two variables (poly(x,2)1 and poly(x,2)2, the linear and orthogonal squared term), but in the original data there is only one x variable? I looked into the sim_slopes documentation and there is a mention about polynomial variables under data:

data: Optional, default is NULL. You may provide the data used to fit the model. This can be a better way to get mean values for centering and can be crucial for models with variable transformations in the formula (e.g., log(x)) or polynomial terms (e.g., poly(x, 2)). You will see a warning if the function detects problems that would likely be solved by providing the data with this argument and the function will attempt to retrieve the original data from the global environment.

However it is unclear how to provide the data. I tried to add the argument data=poly(x, 2) and different variations of this but it still gives me the same out of bounds error message

Anyone has a clue on how to use sim_slope with a polynomial variable?

here is an example to test it out

x <- runif(50, 0, 20)
y<- 0.5*x^2 - x  -1 + rnorm(50, 0, .1)
w <- sample( c("A", "B"), 50, replace=TRUE)
df<- data.frame(x, y, w)

fit=lm(y~poly(x,2)*w, data=df)
sim_slopes(fit, pred=x, modx=w, johnson_neyman = FALSE, data=df)

• Did you supply the dataset that contains the variables you are using to fit the model? If you are not using a dataset (your variables are just in the environment), make a dataset out of them and use it to fit the model and then supply that dataset to the function.
– Noah
Sep 6, 2023 at 23:35
• I did try using the argument data = <name of data frame> but it still doesn't work Sep 7, 2023 at 1:05
• supplying the data x or poly(x,2) does also not work (using my example in the post) Sep 7, 2023 at 1:13

This looks like a bug in the package, and you should bring it up with its author. For now, you can avoid poly(), which is likely the problematic function, by explicitly providing the variable and its square:

fit <- lm(y ~ (x + I(x^2)) * w, data = df)
sim_slopes(fit, pred = 'x', modx = 'w', johnson_neyman = FALSE, data = df)

# SIMPLE SLOPES ANALYSIS
#
# Slope of x when w = B:
#
#   Est.   S.E.    t val.      p
# ------- ------ --------- ------
#   10.69   0.00   2607.11   0.00
#
# Slope of x when w = A:
#
#   Est.   S.E.    t val.      p
# ------- ------ --------- ------
#   10.70   0.00   2793.12   0.00



This produces sensible results for me and they align with those from marginaleffects:

library("marginaleffects")
slopes(fit, variables = "x", newdata = datagrid(w = unique, x = mean))

# Term Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 % w    x
#    x     10.7    0.00383 2794   <0.001 Inf  10.7   10.7 A 11.7
#    x     10.7    0.00410 2607   <0.001 Inf  10.7   10.7 B 11.7


The marginaleffects results are consistent when using poly() instead:

fit <- lm(y ~ poly(x, 2) * w, data = df)
slopes(fit, variables = "x", newdata = datagrid(w = unique, x = mean))

# Term Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 % w    x
#    x     10.7    0.00383 2794   <0.001 Inf  10.7   10.7 A 11.7
#    x     10.7    0.00410 2607   <0.001 Inf  10.7   10.7 B 11.7

• Thanks that helps! I'll try both methods on my data! Sep 8, 2023 at 0:33