title: How can I use lm() with poly() and ggeffect() to plot model effects when variable in quadratic term is not being plotted
body: <p>I am trying to plot model effects using the ggeffects package.  Usually this is no problem. However, I have a model in which there is an interaction between a quadratic term and a first order term (Ex: var_a^2:var_b).  In my actual model there are two quadratic terms and one is involved in an interaction with a third term.  The example below shows the issue with two variables to keep the example as simple as possible.</p>
  <p>If I use the poly() function as follows everything works fine; the model works fine and I can plot model effects. However, I would like to be able to remove he var_a^2:var_b interaction from the model.</p>
  <pre><code># Create a data frame with explanatory variables
  dat &lt;- expand.grid(var_a = c(-1, 0, 1),
                     var_b = c(-1, 0, 1))
  
  # Replicate so there plenty of data points for modelling
  dat &lt;- bind_rows(dat, dat)
  
  # Create response value
  dat$resp &lt;- dat$var_a + 4 * dat$var_a^2 + dat$var_a * dat$var_b + rnorm(nrow(dat), 0, 0.2)
  
  
  mod1 &lt;- lm(resp ~ poly(var_a, 2) * var_b, data = dat)
  summary(mod1)
  
  plot(ggeffect(mod1, c('var_a', 'var_b')))
  plot(ggeffect(mod1, c('var_b')))
  </code></pre>
  <p>Here's the plot of the var_a:var_b interaction
  <a href="https://i.stack.imgur.com/7xCcJ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/7xCcJ.png" alt="a:b interaction plot" /></a></p>
  <p>Here's a plot of effect var_b alone.  I wouldn't normally plot var_b without var_a since they are involved in an interaction, but this is just to show my issue.  It works here, but not below.
  <a href="https://i.stack.imgur.com/opKPc.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/opKPc.png" alt="Effect b plot" /></a></p>
  <p>As stated above, I would like to remove the var_a^2:var_b interaction.  As far as I know, to do that I need to separate the first and second order poly() terms as done below.</p>
  <pre><code># This is the same as mod 1 without the var_a^2:var_b interaction
  # I've split up the poly() terms so I can remove any unwanted terms
  mod2 &lt;- lm(resp ~ poly(var_a, 1) * var_b + poly(var_a, 2)[, 2], data = dat)
  summary(mod2)
  </code></pre>
  <p>I can plot the var_a:var_b effect and get the essentially the same result as above.</p>
  <pre><code>plot(ggeffect(mod2, c('var_a', 'var_b')))
  </code></pre>
  <p>But I can no longer plot the var_b effect alone.</p>
  <pre><code>plot(ggeffect(mod2, 'var_b'))
  </code></pre>
  <p><a href="https://i.stack.imgur.com/PiVPj.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/PiVPj.png" alt="output from mod2 b" /></a></p>
  <p>For what it's worth, I believe the ggeffect() function is not working because the effects package on which it relies sets all variables other than the one being plotted (var_b) to a constant value which does not allow the quadratic effect of var_a to be calculated.</p>
  <p>Is there a way to use the poly() function in a way that allows me to remove unwanted interactions from the model while still being able to use ggeffect() to plot model effects?</p>
  <p>I tried to add ggeffects and effects tags to this question but don't have enough points to create a new tag. If anyone would like to do that, please do.</p>