1
$\begingroup$

I have fitted a generalized additive model with an interaction term using the gam function from the mgcv package in R. I would like to perform the default Type II analysis used by the Anova() function from the car package to check which variable is significantly associated with the outcome. It is my understanding that car::Anova() is a useful function for any type of model where a single predictor is involved in multiple terms (e.g., non-linear terms or interactions). However, when I run car::Anova() on my interaction model, I receive an error message. I would like to know if car::Anova() is the right test to use for generalized additive models, and if not, what alternative tests I should consider? Many thanks

Here is my model

mod1 <-
  gam(
    disease_severity ~  te(min_rh,  daily_minimum_temperature, k = 4) +  te(max_ws, rain_per_rainy_day, k = 5),
    family = betar(),
    method = "REML",
    data = dat_season
  )

summary(mod1)

Here is the output:

Family: Beta regression(8.84) 
Link function: logit 

Formula:
disease_severity ~ te(min_rh, daily_minimum_temperature, k = 4) + 
    te(max_ws, rain_per_rainy_day, k = 5)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.03709    0.12465  -0.298    0.766

Approximate significance of smooth terms:
                                       edf Ref.df Chi.sq p-value    
te(min_rh,daily_minimum_temperature) 3.000   3.00  39.62  <2e-16 ***
te(max_ws,rain_per_rainy_day)        6.608   6.93  78.30  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.867   Deviance explained = 92.9%
-REML = -45.944  Scale est. = 1         n = 41

I'd like to know which predictors is significantly associated with the outcome of the interaction model, but I get an error. I wonder if this test can't be used for gams?

car::Anova(mod1)

Here is the error message:

Error in glm.control(nthreads = 1, ncv.threads = 1, irls.reg = 0, epsilon = 1e-07, : unused arguments (nthreads = 1, ncv.threads = 1, irls.reg = 0, mgcv.tol = 1e-07, mgcv.half = 15, rank.tol = 1.49011611938477e-08, nlm = list(7, 1e-06, 2, 1e-04, 200, FALSE), optim = list(1e+07), newton = list(1e-06, 5, 2, 30, FALSE), outerPIsteps = 0, idLinksBases = TRUE, scalePenalty = TRUE, efs.lspmax = 15, efs.tol = 0.1, keepData = FALSE, scale.est = "fletcher", edge.correct = FALSE)

$\endgroup$
2
  • $\begingroup$ "I'd like to know which predictors is significantly associated with the outcome of the interaction model" the summary table tells you that; your model has two terms and both are estimated such that their effects are unlikely under the null hypothesis of no effect. You'll have to expand on what you want to know if you are to get useful advice. $\endgroup$ Mar 10, 2023 at 19:17
  • $\begingroup$ @GavinSimpson Here is the discussion in the comment section stats.stackexchange.com/questions/603155/… $\endgroup$
    – Ahsk
    Mar 11, 2023 at 1:19

1 Answer 1

2
$\begingroup$

It is my understanding that car::Anova() is a useful function for any type of model where a single predictor is involved in multiple terms (e.g., non-linear terms or interactions).

That's true for many types of models, but a GAM is fit differently from the type of model covered on the page you link. I don't think that car::Anova() can handle your GAM, which uses penalization to trade off the flexibility of the fit against the amount of data available.

You will notice that coefficients aren't reported for the smooths in your GAM model. There is, hiding within the model, effectively a large set of (penalized) coefficients for each smooth, with a Wald test on the entire smooth evaluating the overall significance reported. Within each of your tensor-product smooths, that set of coefficients includes what you might consider all the "main" and "interaction" coefficients involving the included predictors.

Conceptually, the displayed Wald test on each smooth thus accomplishes what a Wald Type II Anova would accomplish in a different type of model: evaluating a combination of multiple coefficient estimates. So there's no need to use something like car::Anova() for this model. You already have the equivalent.

The mgcv package provides an anova.gam() function appropriate to its GAM models. That would be the best choice for evaluating terms in a single model, or for comparing nested GAM models. See its help page for cautions about its use.

$\endgroup$
4
  • $\begingroup$ Thanks. On a side note, how does anova.gam() help choosing a best fit model when none of the models gives a significant variation? I'd think AIC, BIC and R2 should be compared then? For instance, I compared five models using `anova.gam(m1, m2, m3, m4, m5), and got this result. Resid. Df Resid. Dev Df Deviance 1 30.084 -97.266 2 29.759 -97.211 0.3245305 -0.05553 3 29.790 -97.692 -0.0308177 0.48075 4 29.707 -97.370 0.0831208 -0.32170 5 29.713 -97.678 -0.0063971 0.30776 $\endgroup$
    – Ahsk
    Mar 11, 2023 at 21:45
  • $\begingroup$ Hard to decide which one is better based on anova stats.stackexchange.com/questions/608364/… $\endgroup$
    – Ahsk
    Mar 11, 2023 at 21:45
  • $\begingroup$ @Ahsk anova() with multiple models only makes sense if the models are nested. See the help page for anova.glm, which is what anova.gam uses. If models aren't nested, then AIC or BIC could be used. For $R^2$ there are many different types if you aren't doing ordinary least squares, so be cautious. $\endgroup$
    – EdM
    Mar 11, 2023 at 22:48
  • $\begingroup$ thank you very much $\endgroup$
    – Ahsk
    Mar 12, 2023 at 2:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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