GAM Analysis – How to Perform Type II Analysis on a GAM Interaction Model in R When car::Anova() Gives an Error

anovabiostatisticsgeneralized-additive-modelmultivariate analysisregression

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)

Best Answer

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.

Related Solutions

Solved – How to determine the type of spline in GAM

In the documentation for the mgcv package, there is a page describing the spline-based smoothers available. Moreover Wood (the package author) offers the following advice:

Broadly speaking the default penalized thin plate regression splines tend to give the best MSE performance, but they are slower to set up than the other bases.

Since in your case you have less than 200 data points, so I don't think you will run into computational issues with the default method. In section 4.1 of Wood's book "Generalized Additive Models: an introduction with R", he has a summary of the major smoothing bases (Thin plate regression splines, Duchon splines, Cubic regression splines, P-splines) available in mgcv along with a discussion of their merits and other practical considerations. I have found the book quite helpful in developing my understanding of GAMs.

Solved – How to choose the type of GAM-parameters

I'm assuming this is better explained in the 2nd edition of Simon's book (which should be out in a couple of days) as he and his students only worked out some of the theory for this years after Simon wrote his book.

What Marra & Wood (2011) showed was that if we want to do selection on a model with smooth terms, then one very good approach is to add an extra penalty to all the smooth terms. This additional penalty works with the smoothness penalty for that term to control both the wiggliness of the term and whether a term should be in the model at all.

So, unless you have any good theory to assume either smooth or linear/parametric forms/effects for the covariates, you could approach the problem as choosing among all models (representable by the additive combination of linear combinations of the basis functions) between one with smooths of each covariate all the way back to a model containing just an intercept.

For example:

library(mgcv)
data(trees)
ct1 <- gam(log(Volume) ~ s(Height) + s(Girth), data=trees, method = "REML", select = TRUE)

> summary(ct1)

Family: gaussian 
Link function: identity 

Formula:
log(Volume) ~ s(Height) + s(Girth)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.27273    0.01492   219.3   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
            edf Ref.df      F  p-value    
s(Height) 0.967      9  3.249 3.51e-06 ***
s(Girth)  2.725      9 75.470  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.975   Deviance explained = 97.8%
-REML = -23.681  Scale est. = 0.0069012  n = 31

Looking at the output (specifically in the section Parametric coefficients), we note that both terms are highly significant. But note the effective degrees of freedom value for the smooth of Height; it is ~1. What these tests are doing is explained in Wood (2013).

This suggests to me that Height should enter the model as a linear parametric term. We can evaluate this by plotting the fitted smooth:

> plot(ct1, select = 1, shade = TRUE, scale = 0, seWithMean = TRUE)

which gives:

This clearly shows that the selected form of the effect of Height is linear.

However, if you didn't know this upfront (and you didn't because otherwise you wouldn't have asked the question) you can't now refit the model to these data using only a linear term for Height. That would cause you real problems with inference down the line. The output in summary() has accounted for the fact that you did this selection. If you refit the model with a linear parametric effect of Height, the output wouldn't know this and you'd get overly optimistic p-values.

As for question 2, as already mentioned in the comments, no, don't exponentiate the coefficients from this model. Also, don't delve into fitted models as the contents of these components is not always what you might expect. Use the extractor functions instead; in this case coef().

Later in the book when Simon gets to GLMs and GAMs, you'll see him model these data via a Gamma GLM:

ct1 <- gam(Volume ~ Height + s(Girth), data=trees, method = "REML",
           family = Gamma(link = "log"))

In that model, because the fitting is being done on the scale of the linear predictor (on the log scale), the coefficients could be exponentiated to get some partial effect, but you are better off using predict(ct1, ...., type = "response") to get back fitted values/predictions on the scale of the response (in m^3).

Marra, G. & Wood, S. N. Practical variable selection for generalized additive models. Comput. Stat. Data Anal. 55, 2372–2387 (2011).

Wood, S. N. On p-values for smooth components of an extended generalized additive model. Biometrika 100, 221–228 (2013).

Best Answer

Related Solutions

Solved – How to determine the type of spline in GAM

Solved – How to choose the type of GAM-parameters

Related Question