[R-sig-ME] Confidence intervals in GAMM4
Joshua Wiley
jwiley.psych at gmail.com
Tue Jun 11 18:42:35 CEST 2013
Hi Rodrigo,
If you have some knowledge about relations among the residuals, rather
than put a structure on the residuals, I would propose you try to it
into the model.
If within your level 2 units, you believe there are residual
correlations over time, this argues for a time effect in the model.
If you already have time in the model, and the residuals are still
correlated, then I would say either:
1) The treatment of time is insufficient --- for example using a
linear effect when you need some other functional form (since you are
using gams anyway, why not put a smooth term on time?)
2) Making the time effect a random effect
This should effectively take an effect that used to have nowhere to go
but the residuals and explicitly put it into the model.
Cheers,
Josh
On Tue, Jun 11, 2013 at 9:08 AM, Gavin Simpson <gavin.simpson at ucl.ac.uk> wrote:
> On Tue, 2013-06-11 at 09:02 -0300, Rodrigo Tardin wrote:
>> Hi Gavin and other member of the list
>>
>> Thanks a lot for your response. You are right about the correlation
>> structure. I was not aware, thank you.
>> One other question that may look like very basic:
>> Without the correlation structure (that as you said, does not exist in
>> GAMM4 or lme4) in the mixed model, does it account for autocorrelation or
>> no, without any specification of correlation structure it does not account
>> for autcorrelation in the residuals. Because my data do have a problem of
>> autocorrelation on the residuals.
>
> There is an induced correlation due to the random effect - all
> observations within a level of RANDOM have the same estimated
> correlation, IIRC. The AR(1) you hoped to specify would have that
> correlation decline by the absolute power of the separation in time
> (i.e. exponentially), so they are different.
>
> If the residuals are correlated then you should do something about it if
> you can as it may change which terms are significant in the model (as
> standard errors are to narrow), and also with GAMMs where the smoothness
> of the splines is determined from the data the procedure may tend to
> under-smooth the data because it assumes that they are less correlated
> than they really are.
>
> HTH
>
> G
>
>> Thanks in advance
>> Rodrigo
>>
>>
>> 2013/6/5 Gavin Simpson <gavin.simpson at ucl.ac.uk>
>>
>> > On Thu, 2013-05-30 at 12:12 -0300, Rodrigo Tardin wrote:
>> > > Hi all,
>> > >
>> > <snip />
>> > >
>> > > I searched in R for confidence intervals in GAMM4 but I did not find it.
>> > I
>> > > could obtain variance and std deviation for the random and fixed effects
>> > > Groups Name Variance Std.Dev.
>> > > RANDOM (Intercept) 1.4973e+01 3.8694319
>> > > Xr.0 s(DISTCOAST) 4.7361e-02 0.2176263
>> > > Xr s(DEPTH) 1.9779e-06 0.0014064
>> > >
>> > >
>> > > Here it is my model.
>> > > n3 <- gamm4(OCC_BIN~s(DEPTH)+s(DISTANCE_TO_COAST)+offset(RT),random = ~
>> > > (1|RANDOM),correlation=corAR1(),method="ML", family=binomial,data=bryde3)
>> >
>> > Unfortunately, this is *totally* wrong. There is no `correlation`
>> > argument in `gamm4()` nor `glmer()`, which is the underlying fitting
>> > function. That this didn't raise an error is due to `gamm4()` etc having
>> > argument `...` which silently mops up any left over, non-used arguments.
>> >
>> > `gamm()` in package *mgcv* does have a `correlation` argument but that
>> > will fit your binomial GLMM via PQL which isn't such a good solution for
>> > such models.
>> >
>> > In nlme:::lme there was a function intervals() which could provide the
>> > CI on the REs - see if there is an equivalent for lme4:::glmer. IIRC
>> > Doug has something on this in his in-prep book on mixed effects models
>> > via lme4, see chapter 1 in
>> > http://lme4.r-forge.r-project.org/lMMwR/lrgprt.pdf
>> >
>> > That presumes that you sort out the AR(1) business; you can't do that in
>> > glmer().
>> >
>> > Options are to move back to mgcv::gamm() but as I said, PQL isn't great
>> > of binomial models. If your REs are simple, then mgcv::gam() can be
>> > used. Again this doesn't have a correlation argument but ?magic (after
>> > loading mgcv) has an example of including the correlation in the fit via
>> > some jiggery-pokery. Alternatively and related to gam() is mgcv::bam()
>> > which can take a known AR(1) parameter into account during the fitting.
>> > This fits the AR(1) given the ordering of the data, which is what
>> > `correlation=corAR1()` would have done - perhaps this just plain won't
>> > work in nlme:::lme if you don't specify any ordering variable? - but I'm
>> > not sure that will be correct given that your data are strictly time
>> > ordered.
>> >
>> > HTH
>> >
>> > G
>> >
>> > > The OFFSET is the boat route (the number of times the boat searched for
>> > > whales in each 2x2km grid)
>> > > RANDOM is the individual whale as done in previous studies (Hazen et al
>> > > 2009 - MEPS - doi:10.3354/meps08108)
>> > >
>> > > Can someone help me, please?
>> > > Sincerely,
>> > > _______________________________________________
>> > > R-sig-mixed-models at r-project.org mailing list
>> > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> >
>> > --
>> > Gavin Simpson, PhD [t] +1 306 337 8863
>> > Adjunct Professor, Department of Biology [f] +1 306 337 2410
>> > Institute of Environmental Change & Society [e] gavin.simpson at uregina.ca
>> > 523 Research and Innovation Centre [tw] @ucfagls
>> > University of Regina
>> > Regina, SK S4S 0A2, Canada
>> >
>> >
>> >
>> >
>>
>>
>> Hi Gavin and other member of the list
>>
>>
>> Thanks a lot for your response. You are right about the correlation
>> structure. I was not aware, thank you.
>> One other question that may look like very basic:
>> Without the correlation structure (that as you said, does not exist in
>> GAMM4 or lme4) in the mixed model, does it account for autocorrelation
>> or no, without any specification of correlation structure it does not
>> account for autcorrelation in the residuals. Because my data do have a
>> problem of autocorrelation on the residuals.
>>
>>
>> Thanks in advance
>> Rodrigo
>>
>>
>> 2013/6/5 Gavin Simpson <gavin.simpson at ucl.ac.uk>
>> On Thu, 2013-05-30 at 12:12 -0300, Rodrigo Tardin wrote:
>> > Hi all,
>> >
>> <snip />
>> >
>> > I searched in R for confidence intervals in GAMM4 but I did
>> not find it. I
>> > could obtain variance and std deviation for the random and
>> fixed effects
>> > Groups Name Variance Std.Dev.
>> > RANDOM (Intercept) 1.4973e+01 3.8694319
>> > Xr.0 s(DISTCOAST) 4.7361e-02 0.2176263
>> > Xr s(DEPTH) 1.9779e-06 0.0014064
>> >
>> >
>> > Here it is my model.
>> > n3 <-
>> gamm4(OCC_BIN~s(DEPTH)+s(DISTANCE_TO_COAST)+offset(RT),random
>> = ~
>> > (1|RANDOM),correlation=corAR1(),method="ML",
>> family=binomial,data=bryde3)
>>
>>
>> Unfortunately, this is *totally* wrong. There is no
>> `correlation`
>> argument in `gamm4()` nor `glmer()`, which is the underlying
>> fitting
>> function. That this didn't raise an error is due to `gamm4()`
>> etc having
>> argument `...` which silently mops up any left over, non-used
>> arguments.
>>
>> `gamm()` in package *mgcv* does have a `correlation` argument
>> but that
>> will fit your binomial GLMM via PQL which isn't such a good
>> solution for
>> such models.
>>
>> In nlme:::lme there was a function intervals() which could
>> provide the
>> CI on the REs - see if there is an equivalent for
>> lme4:::glmer. IIRC
>> Doug has something on this in his in-prep book on mixed
>> effects models
>> via lme4, see chapter 1 in
>> http://lme4.r-forge.r-project.org/lMMwR/lrgprt.pdf
>>
>> That presumes that you sort out the AR(1) business; you can't
>> do that in
>> glmer().
>>
>> Options are to move back to mgcv::gamm() but as I said, PQL
>> isn't great
>> of binomial models. If your REs are simple, then mgcv::gam()
>> can be
>> used. Again this doesn't have a correlation argument
>> but ?magic (after
>> loading mgcv) has an example of including the correlation in
>> the fit via
>> some jiggery-pokery. Alternatively and related to gam() is
>> mgcv::bam()
>> which can take a known AR(1) parameter into account during the
>> fitting.
>> This fits the AR(1) given the ordering of the data, which is
>> what
>> `correlation=corAR1()` would have done - perhaps this just
>> plain won't
>> work in nlme:::lme if you don't specify any ordering variable?
>> - but I'm
>> not sure that will be correct given that your data are
>> strictly time
>> ordered.
>>
>> HTH
>>
>> G
>>
>> > The OFFSET is the boat route (the number of times the boat
>> searched for
>> > whales in each 2x2km grid)
>> > RANDOM is the individual whale as done in previous studies
>> (Hazen et al
>> > 2009 - MEPS - doi:10.3354/meps08108)
>> >
>> > Can someone help me, please?
>> > Sincerely,
>>
>> > _______________________________________________
>> > R-sig-mixed-models at r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>> --
>> Gavin Simpson, PhD [t] +1 306 337
>> 8863
>> Adjunct Professor, Department of Biology [f] +1 306 337
>> 2410
>> Institute of Environmental Change & Society [e]
>> gavin.simpson at uregina.ca
>> 523 Research and Innovation Centre [tw] @ucfagls
>> University of Regina
>> Regina, SK S4S 0A2, Canada
>>
>>
>>
>>
>>
>>
>>
>> --
>> Rodrigo Tardin
>>
>> Doutorando em Ecologia e Conservação - IBRAG - UERJ
>> Mestre em Biologia Animal - PPGBA - UFRRJ
>> Especialista em Docência do Ensino Superior - IAVM
>> Laboratório de Bioacústica e Ecologia de Cetáceos - UFRRJ/ IF/ DCA
>
> --
> Gavin Simpson, PhD [t] +1 306 337 8863
> Adjunct Professor, Department of Biology [f] +1 306 337 2410
> Institute of Environmental Change & Society [e] gavin.simpson at uregina.ca
> 523 Research and Innovation Centre [tw] @ucfagls
> University of Regina
> Regina, SK S4S 0A2, Canada
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://joshuawiley.com/
Senior Analyst - Elkhart Group Ltd.
http://elkhartgroup.com
More information about the R-sig-mixed-models
mailing list