I'm using a Generalized Additive Model (GAM) to assess the effect of longitudinal changes in proportion of forests on abundance of skunks. My GAM has the following structure where X are the X coordinates of the centroids of trapping sites.

mod <- gam(nb_unique ~ s(x,prop_forest), offset=log_trap_eff, family=nb(theta=NULL, link="log"), data=succ_capt_skunk, method = "REML", select = TRUE) 

Family: Negative Binomial(13.446) 
Link function: log 

nb_unique ~ s(x, prop_forest)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.02095    0.03896  -51.87   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                   edf Ref.df Chi.sq  p-value    
s(x,prop_forest) 3.182     29  17.76 0.000102 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.37   Deviance explained =   49%
-REML = 268.61  Scale est. = 1         n = 58

Should I test for the simple effects of covariates X and prop_forest in the GAM when the interaction is significant ?


There are a couple of issues here:

  1. The s() smoother is a 2d spline that assume the same degree of smoothness in both directions. For 2-d smooths of variables measured on different scales, as you have here, you should use a tensor product smooth: te(x, prop_forest)

  2. The output, I would argue, doesn't say that the interaction is significant. It only says that the effect captured/modelled by the 2-d spline is inconsistent with the null function. If you want to do an ANOVA-like decomposition into main effects and interaction you need the special ti() smoothers.

    For the current model then you would fit:

    nb_unique ~ s(x) + s(prop_forest) + ti(x, prop_forest)

    and use the info from summary() to check if the interaction is needed (the row for the ti(x, prop_forest) part). Or you can compute the simpler model

    nb_unique ~ s(x) + s(prop_forest)

    and compare the models with a generalised likelihood ratio test (via the anova() method (anova(mod1, mod2)), although this is less reliable than the method used to compute p-values in the summary() output IIRC, or using AIC via the AIC().

    Assuming that the ti() interaction term is significant, then you should refit that model replacing the three terms (two s() main effects and one ti() interaction) with a single te(x, prop_forest) term. The reason for this is that even though the models are equivalent, the way they are implemented implies more penalty matrices, and hence more smoothness parameters to be estimated, for the ANOVA-decomposition model formulation (the one with the ti() term).

So really your question is a little premature given the second point; what you show doesn't indicate a significant interaction per se.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.