[R] mgcv: inclusion of random intercept in model - based on p-value of smooth or anova?
Simon Wood
s.wood at bath.ac.uk
Mon Jun 25 17:31:28 CEST 2012
Martin,
I had a nagging feeling that there must be more to this than my previous
reply suggested, and have been digging a bit further. Basically I would
expect these p-values to not be great if you had nested random effects
(such as main effects + their interaction), but your case looked rather
straightforward...
After some digging it turns out that the p-value for Placename is wrong
in this case, as you suspected. R routine `qr' has a `tol' argument that
by default is set to 1e-7. In the computation of the test statistic for
"re" terms I had called qr without changing this default to the value 0
that is actually appropriate (I should have noticed this, I've been
caught out by it before). With the correct 'tol', Placement is
non-significant. This will be fixed in the next release, but that won't
happen until I've tested the random effects p-value approximations a bit
more thoroughly.
best,
Simon
On 11/06/12 14:56, Martijn Wieling wrote:
> 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
>
--
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603 http://people.bath.ac.uk/sw283
More information about the R-help
mailing list