R Package Conflict Between Gam and Mgcv

R Package conflict between gam and mgcv?

The problem, as you ANTICIPATED, is that both gam and mgcv packages install S3 methods for "gam" objects.
But, as stated in the ?detach documentation:

registered S3 methods from the namespace will not be removed.

So in your case it is easy to see that is the cause of your problems:

library(gam)
# installs gam::print.summary.gam
identical(getS3method('print', 'summary.gam'), gam:::print.summary.gam)
[1] TRUE

library(mgcv)
# installs mgcv::print.summary.gam
identical(getS3method('print', 'summary.gam'), mgcv:::print.summary.gam)
[1] TRUE

# save a pointer before unloading namespaces
mgcv_psgam <- mgcv:::print.summary.gam

detach('package:mgcv',unload = TRUE, character.only = TRUE)
# after the detach, the method from mgcv is still installed !!!
identical(getS3method('print', 'summary.gam'), mgcv_psgam)
[1] TRUE

Conclusion: make sure you never load mgcv

AIC() and model$aic give different result in mgcv::gam()

The discrepancy is because what is stored in $aic takes as the degrees of freedom for the complexity correction in AIC is the effective degrees of freedom (EDF) of the model. This has been shown to be too liberal or conservative and can result in AIC always selecting the more complex model or the simpler model depending on whether a marginal or conditional AIC is used.

There are approaches to correct for this behaviour and mgcv implements the one of Wood et al (2016), which applies a correction to the degrees of freedom. This is done via the logLik.gam() function, which is called by AIC.gam(). This also explains the difference as $aic is the standard AIC without the correction applied and IIRC is a component of the GAM object that significantly predates the work of Wood et al (2016).

As to why you couldn't replicate this with the simple example, that's because the correction requires the use of components of the fit that are only available when the method used to fit is "REML" or "ML" (including "fREML" for bam() and also not when the Extended Fellner Schall or BFGS optimizers are used:

> library('mgcv')
Loading required package: nlme
This is mgcv 1.8-34. For overview type 'help("mgcv-package")'.
> set.seed(2)
> dat <- gamSim(1,n=400,dist="normal",scale=2)
Gu & Wahba 4 term additive model
> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat, method = 'REML')
> AIC(b)
[1] 1698.504
> b$aic
[1] 1696.894

where the default for method is to use GCV.

References

Wood, S.N., Pya, N., Säfken, B., 2016. Smoothing Parameter and Model Selection for General Smooth Models. J. Am. Stat. Assoc. 111, 1548–1563. https://doi.org/10.1080/01621459.2016.1180986

Are there known compatibility issues with R package mgcv? Are there general rules for compatibility?

Loading mgcv as the first package solved my problem ... strange but true.

mgcv: Difference between s(x, by=cat) and s(cat, bs='re')

The factor by model, gam0, contains a separate smooth of x1 for each level of cat, but doesn't include anything specifically for the means of y in each group[*] because it is miss-specified. Compare this with gam1, which has a single smooth of x1 plus group means for the levels of cat.

Even though you generated random data without any smooth or group level effects, the gam0 model is potentially much more complex and flexible a model as it contains 4 separate smooths, each using potentially 9 degrees of freedom. Your gam1 has a single smooth of x1 which uses up to 9 degrees of freedom, plus something between 4 and 0 degrees of freedom for the random effect smooth. gam0 is simply exploiting random variation in the data that can be explained a little bit by those extra potential degrees of freedom. You can see this in the adjusted R-sq.(adj), which is lower for gam0 despite it explaining ~ twice the deviance as does gam1 (not that either is a good amount of deviance explained).

r$> library("gratia")                                                           

r$> smooths(gam0)
[1] "s(x1):cat1" "s(x1):cat2" "s(x1):cat3" "s(x1):cat4"

r$> smooths(gam1)
[1] "s(x1)" "s(cat)"

[*] Note that your by model should be

gam0 <- gam(y ~ cat + s(x1, by=cat), data=gam.df)

because the smooths created by s(x1, by=cat) are subject to an identifiability constraint (as there's a constant term — the intercept — in the model). This constraint is a sum-to-zero constraint which means that the individual smooths do not contain the group means. This forces the smooths to not only model the way Y changes as a function of x1 in each group but also model the magnitude of Y in the respective groups, but without functions in the span of the basis that could model such constant (magnitude) effects.



Related Topics



Leave a reply



Submit