Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Interaction term interpretation and prediction #147

Closed
solenemarion opened this issue Oct 24, 2022 · 7 comments
Closed

Interaction term interpretation and prediction #147

solenemarion opened this issue Oct 24, 2022 · 7 comments

Comments

@solenemarion
Copy link

Hi there,
I am running an HMSC model with an interaction term between multiple covariates and I wonder if there is an easy way to interpret HMSC output for interaction? For example, I was not able to use the prediction functions for my interaction covariates.
Any suggestion?
Thanks

@gtikhonov
Copy link
Member

You mean the prediction functions do not work at all with interaction covariates, or you are struggling to make them working in a manner that you want them to? Could you be a bit more specific?
An even better way would be to provide a minimal working example demonstrating your problem. Not necessarily with you data - quite the contrary actually, since simulated data that sufficiently represents your modelling task/questions would be the best.

@solenemarion
Copy link
Author

solenemarion commented Oct 25, 2022

thanks for your quick answer, what I mean is that I am not sure how to interpret my beta parameter estimates and the beta plot when using an interaction term between 2 covariates. So I thought I will use the prediction function in the same way as in normal regression modelling but I tried to use the "constructGradient" function with the focal variable being the interaction covariate but it doesn't seem to work. I am not sure if the construct gradient is the best way to interpret interaction terms but I am just not sure about what other way will be the best... sorry if it is a very stupid question! thanks

For example:

set.seed(1)
n = 50
x = rnorm(n)
x2 = rnorm(n)

beta1 = 0
beta2 = 1
L = beta1 + beta2*x

y1 = L + rnorm(n, sd = 1)

Y = as.matrix(y1)
XData = data.frame(x = x, x2=x2)

m = Hmsc(Y = Y, XData = XData, XFormula = ~ x*x2)

nChains = 2
thin = 5
samples = 1000
transient = 500*thin
verbose = 500*thin

m = sampleMcmc(m, thin = thin, samples = samples,
                      transient = transient, nChains = nChains,
                      verbose = verbose)

Gradient = constructGradient(m, focalVariable = "x:x2") # where here I want to visualise the prediction of the interaction term.

@jarioksa
Copy link
Collaborator

jarioksa commented Oct 25, 2022

constructGradient really cannot do this. It produces a data frame which has variables needed to build the model matrix. That model frame contains variables that the formula needs to build the model matrix which is used internally. For instance, for factor variables it will contain the original factors, but not their contrasts (which are used internally in model matrix). For terms like ~ poly(x, 2) it will contain only x and the polynomial is built internally. Same for your interaction term: your formula ~x*x2 contains only two variables x and x2 and these are internally expanded to terms of model matrix ((Intercept), x, x2, x:x2).

There is a clear reason for this. When you set variables to arbitrary values in model frame, the internally calculated variables should reflect these arbitrary values. That is, term x:x2 should be the product of your fixed arbitrary values x and x2. There is no way of doing this other way round: there is an infinite number of pairs of x and x2 that give you the same product x:x2. You don't have variable x:x2, but only variables x and x2, and x:x2 is calculated from these variables internally.

You cannot inspect the interaction alone, and therefore your approach was doomed originally (even if some magic would make it work in constructGradient). Interaction is an ... uhh ... interaction and it only makes sense when inspected together with main effects. My feeling is that constructGradient is too simple a tool for effective inspection of interactions, but you should either compare the results against similar model without interaction, or craft yourself a more appropriate design, such as a grid of x and x2 and plot the results as a surface or perspective plot. Even that can be difficult to understand and interpret.

There is an option of setting non.focalVariables to a constant values (non.focalVariables = list(variable = list(type = 3, value))) in which case the plot against x is linearly related to the plot against x:x2, but you should do this to several values of x2 how the interaction influences the results.

Numerical interactions terms can be difficult to interpret, but they may be the best choice in this case.

@solenemarion
Copy link
Author

OK, thanks for this! yes, I thought construct gradient was not the way to go but thought I will ask in case I missed something.

I guess my question will be then how to interpret output from plot beta with my interaction term? I am guessing the plot beta here is irrelevant and the comparison method with a model without interaction term is maybe the way to go (as suggested above)?

But then I still have the same issue of interpretation, if the model with interaction is different/"better" than the other model, then what is the explanation of this? For example, if my variables are "Distance to feature" and "Elevation" if my model with an interaction term between the two covariates fits better than the model without an interaction term, then how do I interpret the results of this model: with more distance, higher elevation or lower elevation, we detect fewer/more species?
Thanks

@solenemarion
Copy link
Author

I think my issue is more about the interpretation of the beta results itself because the interaction makes it hard to interpret. For example, that is one plot I got with interaction terms between multiple covariates. However, I am not sure what is the best interpretation of this result. For example, for the interaction term Distance:Management (mgt), I get a lot of negative responses but I am not sure how to interpret that: is it when both Distance and management increase that you do get not detection of species?
Thanks for the help
Solene
image

@ovaskain
Copy link
Collaborator

ovaskain commented Oct 26, 2022 via email

@solenemarion
Copy link
Author

Great thanks! That helps a lot.
I will give it a go with some plots!
Cheers
Solene

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants