# How to interpret interactions involving quadratic terms in a GLM?

I'm trying to interpret the output of a GLM from a tutorial that models species abundances as a function of environmental data and species traits.

The function traitglm comes from the mvabund package and in the author's words:

"It takes three arguments: L; R; and Q, and uses these to construct a regression model for abundance (L) as a function of additive terms for each variable in R and Q, quadratic terms for any continuous covariates (such as Weber’s length), and interaction terms between environmental and trait variables."

Here's the code for the model:

library(mvabund)
ft=traitglm(antTraits$$abund,antTraits$$env,antTraits$$traits,method="manyglm") ft$$fourth # print fourth corner terms


And here's the author's note on how to interpret the coefficients:

"Because all predictors are standardised, you can interpret the size of coefficients as a measure of importance. As interaction terms, the fourth coefficients each have an interpretation as the amount by which a unit (1 sd) change in the trait variable changes the slope of the relationship between abundance and a given environmental variable. (But not that there are quadratic terms in the model too, which makes it a little messier to interpret, a plot would be a good idea.) So for example, one of the strongest interactions seems to be between Pilosity and Canopy cover - when there is more canopy cover, ants tend to have lower pilosity scores (less hairy)."

I'm confused by the last sentence, because the coefficient is based on a quadratic function of canopy cover which is a continuous variable.

Here's the full set of results in case you can't run the model:

results <- structure(c(-0.00534511447773003, -0.0447629598020786, 0.106170620522886,
0.156827334590098, 0.114907563126611, 0.104301091980924, 0.0495441802223576,
0.0883899778038507, -0.0337263993908883, -0.0561717751546072,
0.278956272132365, 0.216643904693244, -0.061217876467139, 0.0548056188010987,
0.0357985850266346, 0.0865376639056529, -0.084668676986989, -0.0889573039665969,
0.210179645850993, 0.243034306663744, 0.143412679809558, 0.0343629954232624,
0.0328636340019757, -0.0484200427818454, 0.166704472957665, -0.119947322629352,
0.321139481278101, 0.273490092239953, 0.198330953622823, -0.0192768490596816,
-0.00916876109019902, -0.0148181631062852, -0.0259056584320059,
0.0396076363521851, 0.0390951350868973, 0.0436559367092447, 0.0674481033620884,
0.0790378396927535, 0.0290644391611984, -0.0548418395761887), dim = c(8L,
5L), dimnames = list(c("Femur.length", "No.spines", "Pilosity1",
"Pilosity2", "Pilosity3", "Polymorphism1", "Polymorphism2", "Webers.length"
), c("Bare.ground", "Canopy.cover", "Shrub.cover", "Volume.lying.CWD",
"Feral.mammal.dung")))

• The ft\$fourth here differs from that in the tutorial, which reports values for the reference levels ("0") of Pilosity and Polymorphism that don't show up in your results. (Maybe intercept was omitted.) You can specify "the amount by which a unit (1 sd) change in the trait variable [pilosity] changes the slope of the relationship between abundance and a given environmental variable [Canopy.cover]," but with a quadratic you need to know the baseline around which the 1 sd change was evaluated, and that matrix entry would change with that baseline choice.
– EdM
Commented Jul 7, 2023 at 18:18