[R] mgcv: inclusion of random intercept in model - based on p-value of smooth or anova?
Martijn Wieling
wieling at gmail.com
Mon Jun 11 15:56:32 CEST 2012
Dear Simon,
I ran an additional analysis using bam (mgcv 1.7-17) with three random
intercepts and no non-linearities, and compared these to the results
of lmer (lme4). Using bam results in a significant random intercept
(even though it has a very low edf-value), while the lmer results show
no variance associated to the random intercept of Placename. Should I
drop the random intercept of Placename and if so, how is this apparent
from the results of bam?
Summaries of both models are shown below.
With kind regards,
Martijn
#### l1 = lmer(RefPMIdistMeanLog.c ~ geogamfit + (1|Word) + (1|Key) +
(1|Placename), data=wrddst); print(l1,cor=F)
Linear mixed model fit by REML
Formula: RefPMIdistMeanLog.c ~ geogamfit + (1 | Word) + (1 | Key) + (1
| Placename)
Data: wrddst
AIC BIC logLik deviance REMLdev
-44985 -44927 22498 -45009 -44997
Random effects:
Groups Name Variance Std.Dev.
Word (Intercept) 0.0800944 0.283009
Key (Intercept) 0.0013641 0.036933
Placename (Intercept) 0.0000000 0.000000
Residual 0.0381774 0.195390
Number of obs: 112608, groups: Word, 357; Key, 320; Placename, 40
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.00342 0.01513 -0.23
geogamfit 0.99249 0.02612 37.99
#### m1 = bam(RefPMIdistMeanLog.c ~ geogamfit + s(Word,bs="re") +
s(Key,bs="re") + s(Placename,bs="re"), data=wrddst,method="REML");
summary(m1,freq=F)
Family: gaussian
Link function: identity
Formula:
RefPMIdistMeanLog.c ~ geogamfit + s(Word, bs = "re") + s(Key,
bs = "re") + s(Placename, bs = "re")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.00342 0.01513 -0.226 0.821
geogamfit 0.99249 0.02612 37.991 <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(Word) 3.554e+02 347 634.716 <2e-16 ***
s(Key) 2.946e+02 316 23.054 <2e-16 ***
s(Placename) 1.489e-04 38 7.282 <2e-16 ***
---
Signif. codes: 0 â***â 0.001 â**â 0.01 â*â 0.05 â.â 0.1 â â 1
R-sq.(adj) = 0.693 Deviance explained = 69.4%
REML score = -22498 Scale est. = 0.038177 n = 112608
On Wed, May 23, 2012 at 11:30 AM, Simon Wood-4 [via R]
<ml-node+s789695n4631060h52 at n4.nabble.com> wrote:
> Having looked at this further, I've made some changes in mgcv_1.7-17 to
> the p-value computations for terms that can be penalized to zero during
> fitting (e.g. s(x,bs="re"), s(x,m=1) etc).
>
> The Wald statistic based p-values from summary.gam and anova.gam (i.e.
> what you get from e.g. anova(a) where a is a fitted gam object) are
> quite well founded for smooth terms that are non-zero under full
> penalization (e.g. a cubic spline is a straight line under full
> penalization). For such smooths, an extension of Nychka's (1988) result
> on CI's for splines gives a well founded distributional result on which
> to base a Wald statistic. However, the Nychka result requires the
> smoothing bias to be substantially less than the smoothing estimator
> variance, and this will often not be the case if smoothing can actually
> penalize a term to zero (to understand why, see argument in appendix of
> Marra & Wood, 2012, Scandinavian Journal of Statistics, 39,53-74).
>
> Simulation testing shows that this theoretical concern has serious
> practical consequences. So for terms that can be penalized to zero,
> alternative approximations have to be used, and these are now
> implemented in mgcv_1.7-17 (see ?summary.gam).
>
> The approximate test performed by anova(a,b) (a and b are fitted "gam"
> objects) is less well founded. It is a reasonable approximation when
> each smooth term in the models could in principle be well approximated
> by an unpenalized term of rank approximately equal to the edf of the
> smooth term, but otherwise the p-values produced are likely to be much
> too small. In particular simulation testing suggests that the test is
> not to be trusted with s(...,bs="re") terms, and can be poor if the
> models being compared involve any terms that can be penalized to zero
> during fitting. (Although the mechanisms are a little different, this is
> similar to the problem we would have if the models were viewed as
> regular mixed models and we tried to use a GLRT to test variance
> components for equality to zero).
>
> These issues are now documented in ?anova.gam and ?summary.gam...
>
> Simon
>
> On 08/05/12 15:01, Martijn Wieling wrote:
>
>> Dear useRs,
>>
>> I am using mgcv version 1.7-16. When I create a model with a few
>> non-linear terms and a random intercept for (in my case) country using
>> s(Country,bs="re"), the representative line in my model (i.e.
>> approximate significance of smooth terms) for the random intercept
>> reads:
>> edf Ref.df F p-value
>> s(Country) 36.127 58.551 0.644 0.982
>>
>> Can I interpret this as there being no support for a random intercept
>> for country? However, when I compare the simpler model to the model
>> including the random intercept, the latter appears to be a significant
>> improvement.
>>
>>> anova(gam1,gam2,test="F")
>> Model 1: ....
>> Model 2: .... + s(BirthNation, bs="re")
>> Resid. Df Resid. Dev Df Deviance F Pr(>F)
>> 1 789.44 416.54
>> 2 753.15 373.54 36.292 43.003 2.3891 1.225e-05 ***
>>
>> I hope somebody could help me in how I should proceed in these
>> situations. Do I include the random intercept or not?
>>
>> I also have a related question. When I used to create a mixed-effects
>> regression model using lmer and included e.g., an interaction in the
>> fixed-effects structure, I would test if the inclusion of this
>> interaction was warranted using anova(lmer1,lmer2). It then would show
>> me that I invested 1 additional df and the resulting (possibly
>> significant) improvement in fit of my model.
>>
>> This approach does not seem to work when using gam. In this case an
>> apparent investment of 1 degree of freedom for the interaction, might
>> result in an actual decrease of the degrees of freedom invested by the
>> total model (caused by a decrease of the edf's of splines in the model
>> with the interaction). In this case, how would I proceed in
>> determining if the model including the interaction term is better?
>>
>> With kind regards,
>> Martijn Wieling
>>
>> --
>> *******************************************
>> Martijn Wieling
>> http://www.martijnwieling.nl
>> [hidden email]
>> +31(0)614108622
>> *******************************************
>> University of Groningen
>> http://www.rug.nl/staff/m.b.wieling
>>
>> ______________________________________________
>> [hidden email] 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
>
> ______________________________________________
> [hidden email] 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.
>
>
> ________________________________
> If you reply to this email, your message will be added to the discussion
> below:
> http://r.789695.n4.nabble.com/mgcv-inclusion-of-random-intercept-in-model-based-on-p-value-of-smooth-or-anova-tp4617585p4631060.html
> To unsubscribe from mgcv: inclusion of random intercept in model - based on
> p-value of smooth or anova?, click here.
> NAML
More information about the R-help
mailing list