why the results from the joint_tests function (emmeans package) do not show one of the interactions of the model?

user_15108978 0
Viewed 283

I run a GLMM_adaptive model (I am doing a resource selection function) and I am using the joint_tests function (emmeans package) to compute joint tests of the terms in the model. The problem is that one of the interactions does not appear in the results.

The model is:

mod.hinc <- mixed_model(fixed = Used ~  scale(ndvi) * season * vegfactor + 
                      scale(ndvi^2) + scale(distance^2) + scale(distance) * season, 
                    random = ~ 1 | id, data = hin.c,
                    family = binomial(link="logit"))

After running the model I run the joint_tests function:

install.packages("emmeans")
library(emmeans)
joint_tests(mod.hinc)

And this is the result:

 joint_tests(mod.hinc)
 model term            df1 df2 F.ratio p.value
 ndvi                    1 Inf  36.465  <.0001
 season                  3 Inf  22.265  <.0001
 vegfactor               4 Inf   4.548  0.0011
 distance                1 Inf  33.939  <.0001
 ndvi:season             3 Inf  13.826  <.0001
 ndvi:vegfactor          4 Inf   8.500  <.0001
 season:vegfactor       12 Inf   6.544  <.0001
 ndvi:season:vegfactor  12 Inf   5.165  <.0001

I cannot find the reason why the interaction scale(distance)*season does not appear in the results.

Any help on that issue is welcome. I can provide more details about the model if is required.

Thank you very much in advance.

Juan

emmeans mixed-models
1 Answers

The short answer is that distance:season is not shown because it came up with zero d.f. for the associated interaction contrasts. You could verify this by running joint_tests(mod.hinc, show0df = TRUE).

Why it has 0 d.f. is less clear. However, that is not the only problem here. You have to be extremely careful with numeric predictors when using joint_tests(); it does not do a model ANOVA; instead, as documented, it constructs a reference grid from the fitted model and performs joint tests of interaction contrasts related to the predictors. With numeric predictors, the results depend on the reference grid used.

In this particular instance, the model includes quadratic effects of ndvi and distance; however, the default reference grid is constructed using the range of the covariates -- only two distinct values. Thus, we can pick up the effects of the overall linear trends, but not the curvature effects implied by the quadratic terms. That's why only 1 d.f. of those factors' main effects are tested. There are really 2 d.f. in the effects of ndvi and distance. In order to capture all of those effects, we need to have at least three distinct values of these covariates in the reference grid. One way (not the only way) to accomplish that is to reduce the covariates to their means, plus or minus 1 SD -- which can be accomplished via this code:

meanpm1sd <- function(x)
    c(mean(x) - sd(x), mean(x), mean(x) + sd(x))

joint_tests(mod.hinc, cov.reduce = meanpm1sd)

This will yield a different set of joint tests that likely will include 2-d.f. tests of ndvi and distance. But I don't know if you will still have some interactions missing due to zero-d.f. dimensionalities.

You can look directly at the estimates being tested in detail if you have any questions about what those effects are. For example, for season:distance,

### construct the needed reference grid once and for all
RG <- ref_grid(mod.h1nc, cov.reduce = meanpm1sd)   

EMM <- emmeans(RG, ~ season * distance)
CON <- contrast(EMM, interaction = "consec")

EMM   ### see estimates
CON   ### see interaction contrasts
test(CON, joint = TRUE)

I hope this helps shed some light on what is going on.

user_39615660
Related