[R] GAM Chi-Square Difference Test

[Simon Wood]

Hi Will,

Your edf interpretation is not quite right. The smooths are subject to a 
centring constraint for identifiability reasons, and this removes a 
degree of freedom, so EDF=1 corresponds to a straight line fit.

On your second point and third points. anova(gamb1.1,gamb1.2, 
test="Chisq") conducts the crude approximate GLRT described in the 
?anova.gam. It will not produce sensible results for comparing two 
models that are essentially identical (nor, although it's not the issue 
here, does it do well for testing random effects for equality to zero). 
mgcv generally tries not to second guess the user, so it won't stop you 
trying to compare two essentially identical models, but that doesn't 
mean that the results will be meaningful.

Basically, the results of anova(gamb1.1,gamb1.2, test="Chisq") will give 
reasonable ball-park p-values when both models can be reasonably well 
approximated by some un-penalized model (approximately matching the 
effective degrees of freedom of the penalized version), and the null 
model at least has a lower EDF than the alternative (better still, one 
would condition on the larger model smoothing parameter estimates, to 
force proper nesting). The test is in any case conditional on the 
smoothing parameter estimates. This can lead to the p-values being 
somewhat too low when smoothing psrameters are highly uncertain (as they 
may be at small sample sizes). Note that ML or REML smoothing parameter 
estimation (see gam argument 'method') tends to lead to best p-value 
performance in simulations.

Generally speaking, the p-values reported by single model calls to 
`anova.gam' or `summary.gam' are better justified than the approximate 
GLRT (but prior to 1.7-19, only for terms which cannot be penalized away 
altogether). Even then, what is reported is conditional on the smoothing 
parameter estimates, so can be too low when these are not well identified.

best,
Simon

On 14/07/12 17:41, wshadish wrote:
> We are using GAM in mgcv (Wood), relatively new users, and wonder if anyone
> can advise us on a problem we are encountering as we analyze many short time
> series datasets. For each dataset, we have four models, each with intercept,
> predictor x (trend), z (treatment), and int (interaction between x and z).
> Our models are
>
> Model 1: gama1.1 <- gam(y~x+z+int, family=quasipoisson) ##no smooths
> Model 2: gama1.2 <- gam(y~x+z+s(int, bs="cr"), family=quasipoisson) ##smooth
> the interaction
> Model 3: gama1.3 <- gam(y~s(x, bs="cr")+z+int, family=quasipoisson) ##smooth
> the trend
> Model 4: gama1.4 <- gam(y~s(x, bs="cr")+z+s(int, bs="cr"),
> family=quasipoisson) ##smooth trend and interaction
>
> We have three questions. One question is simple. We occasionally obtain edf
> =1 and Ref.df=1 for some smoothed predictors (x, int). Because Wood says
> that edf can be interpreted roughly as functional form (quadratic, cubic
> etc) + 1, this would imply x^0 functional form for the predictor, and that
> doesn't make a lot of sense. Does such a result for edf and rdf indicate a
> problem (e.g., collinearity) or any particular interpretation?
>
> The other two questions concern which model fits the data best. We do look
> at the usual various fit statistics (R^2, Dev, etc), but our question
> concerns using the anova function to do model comparisons, e.g.,
>
> anova(gama2.1,gama2.2, test="Chisq").
>
> 1. Is there research on the power of the model comparison test? Anecdotally,
> the test seems to reject the null even in cases that would appear to have
> only small differences. These are not hugely long time series, ranging from
> about 17 to about 49, so we would not have thought them to yield large
> power.
>
> 2. More important, in a few cases, we are getting a result that looks like
> this:
>
> anova(gamb1.1,gamb1.2, test="Chisq")
> Analysis of Deviance Table
>
> Model 1: y ~ x + z + int
> Model 2: y ~ x + z + s(int, bs = "cr")
>    Resid. Df Resid. Dev         Df   Deviance P(>|Chi|)
> 1        30     36.713
> 2        30     36.713 1.1469e-05 1.0301e-05 6.767e-05 ***
>
> We are inclined to think that the significance p value here is simply a
> result of rounding error in the computation of the df difference and
> deviance difference, and that we should treat this as indicating the models
> are not different from each other. Has anyone experienced this before? Is
> our interpretation reasonable?
>
> Thanks to anyone who is able to offer advice.
>
> Will Shadish
>
> --
> View this message in context: http://r.789695.n4.nabble.com/GAM-Chi-Square-Difference-Test-tp4636523.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603               http://people.bath.ac.uk/sw283