

Dear Rhelpers,
I'd like to understand how to test the statistical significance of a
random effect in gamm. I am using gamm because I want to test a model
with an AR(1) error structure, and it is my understanding neither gam
nor gamm4 will do the latter.
The data set includes nine short interrupted time series (single case
designs in education, sometimes called Nof1 trials in medicine) from
one study. They report a proportion as outcome (y), computed from a
behavior either observed or not out of 10 trials per time point. Hence I
use binomial (I believe quasibinomial is not available in gamm). Each
of the nine series has an average of 30 observations give or take (total
264 observations), some under treatment (z) and some not. xc is centered
session number, int is the z*xc interaction. Based on prior work, xc is
also smoothed
Consider, for example, two models, both with AR(1) but one allowing a
random effect on xc:
g1 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial,
correlation=corAR1())
g2 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random =
list(xc=~1),correlation=corAR1())
I include the output for g1 and g2 below, but the question is how to
test the significance of the random effect on xc. I considered a test
comparing the LogLikelihoods, but have no idea what the degrees of
freedom would be given that s(xc) is smoothed. I also tried:
anova(g1$gam, g2$gam)
that did not seem to return anything useful for this question.
A related question is how to test the significance of adding a second
random effect to a model that already has a random effect, such as:
g3 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
list(Case=~1, z=~1),correlation=corAR1())
g4 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
list(Case=~1, z=~1, int=~1),correlation=corAR1())
Any help would be appreciated.
Thanks.
Will Shadish
********************************************
g1
$lme
Linear mixedeffects model fit by maximum likelihood
Data: data
Loglikelihood: 437.696
Fixed: fixed
X(Intercept) Xz Xint Xs(xc)Fx1
0.6738466 2.5688317 0.0137415 0.1801294
Random effects:
Formula: ~Xr  1  g
Structure: pdIdnot
Xr1 Xr2 Xr3 Xr4
Xr5 Xr6 Xr7 Xr8 Residual
StdDev: 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781
0.0004377781 0.0004377781 0.0004377781 1.693177
Correlation Structure: AR(1)
Formula: ~1  g
Parameter estimate(s):
Phi
0.3110725
Variance function:
Structure: fixed weights
Formula: ~invwt
Number of Observations: 264
Number of Groups: 1
$gam
Family: binomial
Link function: logit
Formula:
y ~ s(xc) + z + int
Estimated degrees of freedom:
1 total = 4
attr(,"class")
[1] "gamm" "list"
****************************
> g2
$lme
Linear mixedeffects model fit by maximum likelihood
Data: data
Loglikelihood: 443.9495
Fixed: fixed
X(Intercept) Xz Xint Xs(xc)Fx1
0.720018143 2.562155820 0.003457463 0.045821030
Random effects:
Formula: ~Xr  1  g
Structure: pdIdnot
Xr1 Xr2 Xr3 Xr4
Xr5 Xr6 Xr7 Xr8
StdDev: 7.056078e06 7.056078e06 7.056078e06 7.056078e06 7.056078e06
7.056078e06 7.056078e06 7.056078e06
Formula: ~1  xc %in% g
(Intercept) Residual
StdDev: 6.277279e05 1.683007
Correlation Structure: AR(1)
Formula: ~1  g/xc
Parameter estimate(s):
Phi
0.1809409
Variance function:
Structure: fixed weights
Formula: ~invwt
Number of Observations: 264
Number of Groups:
g xc %in% g
1 34
$gam
Family: binomial
Link function: logit
Formula:
y ~ s(xc) + z + int
Estimated degrees of freedom:
1 total = 4
attr(,"class")
[1] "gamm" "list"

William R. Shadish
Distinguished Professor
Founding Faculty
Mailing Address:
William R. Shadish
University of California
School of Social Sciences, Humanities and Arts
5200 North Lake Rd
Merced CA 95343
Physical/Delivery Address:
University of California Merced
ATTN: William Shadish
School of Social Sciences, Humanities and Arts
Facilities Services Building A
5200 North Lake Rd.
Merced, CA 95343
2092284372 voice
2092284007 fax (communal fax: be sure to include cover sheet)
[hidden email]
http://faculty.ucmerced.edu/wshadish/index.htmhttp://psychology.ucmerced.edu [[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


On Fri, 20130607 at 13:12 0700, William Shadish wrote:
> Dear Rhelpers,
>
> I'd like to understand how to test the statistical significance of a
> random effect in gamm. I am using gamm because I want to test a model
> with an AR(1) error structure, and it is my understanding neither gam
> nor gamm4 will do the latter.
gamm4() can't yes and out of the box mgcv::gam can't either but
see ?magic for an example of correlated errors and how the fits can be
manipulated to take the AR(1) (or any structure really as far as I can
tell) into account.
You might like to look at mgcv::bam() which allows an known AR(1) term
but do check that it does what you think; with a random effect spline
I'm not at all certain that it will nest the AR(1) in the random effect
level.
<snip />
> Consider, for example, two models, both with AR(1) but one allowing a
> random effect on xc:
>
> g1 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial,
> correlation=corAR1())
> g2 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random =
> list(xc=~1),correlation=corAR1())
Shouldn't you specify how the AR(1) is nested in the hierarchy here,
i.e. AR(1) within xc? maybe I'm not following your data structure
correctly.
> I include the output for g1 and g2 below, but the question is how to
> test the significance of the random effect on xc. I considered a test
> comparing the LogLikelihoods, but have no idea what the degrees of
> freedom would be given that s(xc) is smoothed. I also tried:
>
> anova(g1$gam, g2$gam)
gamm() fits via the lme() function of package nlme. To do what you want,
you need the anova() method for objects of class "lme", e.g.
anova(g1$lme, g2$lme)
Then I think you should check if the fits were done via REML and also be
aware of the issue of testing wether a variance term is 0.
> that did not seem to return anything useful for this question.
>
> A related question is how to test the significance of adding a second
> random effect to a model that already has a random effect, such as:
>
> g3 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
> list(Case=~1, z=~1),correlation=corAR1())
> g4 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
> list(Case=~1, z=~1, int=~1),correlation=corAR1())
Again, I think you need anova() on the $lme components.
HTH
G
> Any help would be appreciated.
>
> Thanks.
>
> Will Shadish
> ********************************************
> g1
> $lme
> Linear mixedeffects model fit by maximum likelihood
> Data: data
> Loglikelihood: 437.696
> Fixed: fixed
> X(Intercept) Xz Xint Xs(xc)Fx1
> 0.6738466 2.5688317 0.0137415 0.1801294
>
> Random effects:
> Formula: ~Xr  1  g
> Structure: pdIdnot
> Xr1 Xr2 Xr3 Xr4
> Xr5 Xr6 Xr7 Xr8 Residual
> StdDev: 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781
> 0.0004377781 0.0004377781 0.0004377781 1.693177
>
> Correlation Structure: AR(1)
> Formula: ~1  g
> Parameter estimate(s):
> Phi
> 0.3110725
> Variance function:
> Structure: fixed weights
> Formula: ~invwt
> Number of Observations: 264
> Number of Groups: 1
>
> $gam
>
> Family: binomial
> Link function: logit
>
> Formula:
> y ~ s(xc) + z + int
>
> Estimated degrees of freedom:
> 1 total = 4
>
> attr(,"class")
> [1] "gamm" "list"
> ****************************
> > g2
> $lme
> Linear mixedeffects model fit by maximum likelihood
> Data: data
> Loglikelihood: 443.9495
> Fixed: fixed
> X(Intercept) Xz Xint Xs(xc)Fx1
> 0.720018143 2.562155820 0.003457463 0.045821030
>
> Random effects:
> Formula: ~Xr  1  g
> Structure: pdIdnot
> Xr1 Xr2 Xr3 Xr4
> Xr5 Xr6 Xr7 Xr8
> StdDev: 7.056078e06 7.056078e06 7.056078e06 7.056078e06 7.056078e06
> 7.056078e06 7.056078e06 7.056078e06
>
> Formula: ~1  xc %in% g
> (Intercept) Residual
> StdDev: 6.277279e05 1.683007
>
> Correlation Structure: AR(1)
> Formula: ~1  g/xc
> Parameter estimate(s):
> Phi
> 0.1809409
> Variance function:
> Structure: fixed weights
> Formula: ~invwt
> Number of Observations: 264
> Number of Groups:
> g xc %in% g
> 1 34
>
> $gam
>
> Family: binomial
> Link function: logit
>
> Formula:
> y ~ s(xc) + z + int
>
> Estimated degrees of freedom:
> 1 total = 4
>
> attr(,"class")
> [1] "gamm" "list"
>
>

Gavin Simpson, PhD [t] +1 306 337 8863
Adjunct Professor, Department of Biology [f] +1 306 337 2410
Institute of Environmental Change & Society [e] [hidden email]
523 Research and Innovation Centre [tw] @ucfagls
University of Regina
Regina, SK S4S 0A2, Canada
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Gavin et al.,
Thanks so much for the help. Unfortunately, the command
> anova(g1$lme, g2$lme)
gives "Error in eval(expr, envir, enclos) : object 'fixed' not found
and for bam (which is the one that can use a known ar1 term), the error is
> AR1 parameter rho unused with generalized model
Apparently it cannot run for binomial distributions, and presumably also
Poisson.
I did find a Frequently Asked Questions for package mgcv page that said
"How can I compare gamm models? In the identity link normal errors case,
then AIC and hypotheis testing based methods are fine. Otherwise it is
best to work out a strategy based on the summary.gam"
So putting all this together, I take it that my binomial example will
not support a direct model comparison to test the significance of the
random effects. I'm guessing the best strategy based on the summary.gam
is probably just to compare fit indices like Log Likelihoods.
If anyone has any other suggestions, though, please do let me know.
Thanks so much.
Will Shadish
On 6/7/2013 3:02 PM, Gavin Simpson wrote:
> On Fri, 20130607 at 13:12 0700, William Shadish wrote:
>> Dear Rhelpers,
>>
>> I'd like to understand how to test the statistical significance of a
>> random effect in gamm. I am using gamm because I want to test a model
>> with an AR(1) error structure, and it is my understanding neither gam
>> nor gamm4 will do the latter.
>
> gamm4() can't yes and out of the box mgcv::gam can't either but
> see ?magic for an example of correlated errors and how the fits can be
> manipulated to take the AR(1) (or any structure really as far as I can
> tell) into account.
>
> You might like to look at mgcv::bam() which allows an known AR(1) term
> but do check that it does what you think; with a random effect spline
> I'm not at all certain that it will nest the AR(1) in the random effect
> level.
>
> <snip />
>> Consider, for example, two models, both with AR(1) but one allowing a
>> random effect on xc:
>>
>> g1 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial,
>> correlation=corAR1())
>> g2 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random =
>> list(xc=~1),correlation=corAR1())
>
> Shouldn't you specify how the AR(1) is nested in the hierarchy here,
> i.e. AR(1) within xc? maybe I'm not following your data structure
> correctly.
>
>> I include the output for g1 and g2 below, but the question is how to
>> test the significance of the random effect on xc. I considered a test
>> comparing the LogLikelihoods, but have no idea what the degrees of
>> freedom would be given that s(xc) is smoothed. I also tried:
>>
>> anova(g1$gam, g2$gam)
>
> gamm() fits via the lme() function of package nlme. To do what you want,
> you need the anova() method for objects of class "lme", e.g.
>
> anova(g1$lme, g2$lme)
>
> Then I think you should check if the fits were done via REML and also be
> aware of the issue of testing wether a variance term is 0.
>
>> that did not seem to return anything useful for this question.
>>
>> A related question is how to test the significance of adding a second
>> random effect to a model that already has a random effect, such as:
>>
>> g3 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
>> list(Case=~1, z=~1),correlation=corAR1())
>> g4 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
>> list(Case=~1, z=~1, int=~1),correlation=corAR1())
>
> Again, I think you need anova() on the $lme components.
>
> HTH
>
> G
>
>> Any help would be appreciated.
>>
>> Thanks.
>>
>> Will Shadish
>> ********************************************
>> g1
>> $lme
>> Linear mixedeffects model fit by maximum likelihood
>> Data: data
>> Loglikelihood: 437.696
>> Fixed: fixed
>> X(Intercept) Xz Xint Xs(xc)Fx1
>> 0.6738466 2.5688317 0.0137415 0.1801294
>>
>> Random effects:
>> Formula: ~Xr  1  g
>> Structure: pdIdnot
>> Xr1 Xr2 Xr3 Xr4
>> Xr5 Xr6 Xr7 Xr8 Residual
>> StdDev: 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781
>> 0.0004377781 0.0004377781 0.0004377781 1.693177
>>
>> Correlation Structure: AR(1)
>> Formula: ~1  g
>> Parameter estimate(s):
>> Phi
>> 0.3110725
>> Variance function:
>> Structure: fixed weights
>> Formula: ~invwt
>> Number of Observations: 264
>> Number of Groups: 1
>>
>> $gam
>>
>> Family: binomial
>> Link function: logit
>>
>> Formula:
>> y ~ s(xc) + z + int
>>
>> Estimated degrees of freedom:
>> 1 total = 4
>>
>> attr(,"class")
>> [1] "gamm" "list"
>> ****************************
>> > g2
>> $lme
>> Linear mixedeffects model fit by maximum likelihood
>> Data: data
>> Loglikelihood: 443.9495
>> Fixed: fixed
>> X(Intercept) Xz Xint Xs(xc)Fx1
>> 0.720018143 2.562155820 0.003457463 0.045821030
>>
>> Random effects:
>> Formula: ~Xr  1  g
>> Structure: pdIdnot
>> Xr1 Xr2 Xr3 Xr4
>> Xr5 Xr6 Xr7 Xr8
>> StdDev: 7.056078e06 7.056078e06 7.056078e06 7.056078e06 7.056078e06
>> 7.056078e06 7.056078e06 7.056078e06
>>
>> Formula: ~1  xc %in% g
>> (Intercept) Residual
>> StdDev: 6.277279e05 1.683007
>>
>> Correlation Structure: AR(1)
>> Formula: ~1  g/xc
>> Parameter estimate(s):
>> Phi
>> 0.1809409
>> Variance function:
>> Structure: fixed weights
>> Formula: ~invwt
>> Number of Observations: 264
>> Number of Groups:
>> g xc %in% g
>> 1 34
>>
>> $gam
>>
>> Family: binomial
>> Link function: logit
>>
>> Formula:
>> y ~ s(xc) + z + int
>>
>> Estimated degrees of freedom:
>> 1 total = 4
>>
>> attr(,"class")
>> [1] "gamm" "list"
>>
>>
>

William R. Shadish
Distinguished Professor
Founding Faculty
Mailing Address:
William R. Shadish
University of California
School of Social Sciences, Humanities and Arts
5200 North Lake Rd
Merced CA 95343
Physical/Delivery Address:
University of California Merced
ATTN: William Shadish
School of Social Sciences, Humanities and Arts
Facilities Services Building A
5200 North Lake Rd.
Merced, CA 95343
2092284372 voice
2092284007 fax (communal fax: be sure to include cover sheet)
[hidden email]
http://faculty.ucmerced.edu/wshadish/index.htmhttp://psychology.ucmerced.edu______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


On Tue, 20130611 at 10:08 0700, William Shadish wrote:
> Gavin et al.,
>
> Thanks so much for the help. Unfortunately, the command
>
> > anova(g1$lme, g2$lme)
>
> gives "Error in eval(expr, envir, enclos) : object 'fixed' not found
This is with mgcv:::gamm yes? Strange  did you load nlme first? I think
mgcv now imports from the nlme package so to use its functions/methods
explicitly you need to load nlme  before loading mgcv also loaded nlme,
but it no longer does so.
That *should* work.
HTH
G
> and for bam (which is the one that can use a known ar1 term), the error is
>
> > AR1 parameter rho unused with generalized model
>
> Apparently it cannot run for binomial distributions, and presumably also
> Poisson.
>
> I did find a Frequently Asked Questions for package mgcv page that said
>
> "How can I compare gamm models? In the identity link normal errors case,
> then AIC and hypotheis testing based methods are fine. Otherwise it is
> best to work out a strategy based on the summary.gam"
>
> So putting all this together, I take it that my binomial example will
> not support a direct model comparison to test the significance of the
> random effects. I'm guessing the best strategy based on the summary.gam
> is probably just to compare fit indices like Log Likelihoods.
>
> If anyone has any other suggestions, though, please do let me know.
>
> Thanks so much.
>
> Will Shadish
>
> On 6/7/2013 3:02 PM, Gavin Simpson wrote:
> > On Fri, 20130607 at 13:12 0700, William Shadish wrote:
> >> Dear Rhelpers,
> >>
> >> I'd like to understand how to test the statistical significance of a
> >> random effect in gamm. I am using gamm because I want to test a model
> >> with an AR(1) error structure, and it is my understanding neither gam
> >> nor gamm4 will do the latter.
> >
> > gamm4() can't yes and out of the box mgcv::gam can't either but
> > see ?magic for an example of correlated errors and how the fits can be
> > manipulated to take the AR(1) (or any structure really as far as I can
> > tell) into account.
> >
> > You might like to look at mgcv::bam() which allows an known AR(1) term
> > but do check that it does what you think; with a random effect spline
> > I'm not at all certain that it will nest the AR(1) in the random effect
> > level.
> >
> > <snip />
> >> Consider, for example, two models, both with AR(1) but one allowing a
> >> random effect on xc:
> >>
> >> g1 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial,
> >> correlation=corAR1())
> >> g2 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random =
> >> list(xc=~1),correlation=corAR1())
> >
> > Shouldn't you specify how the AR(1) is nested in the hierarchy here,
> > i.e. AR(1) within xc? maybe I'm not following your data structure
> > correctly.
> >
> >> I include the output for g1 and g2 below, but the question is how to
> >> test the significance of the random effect on xc. I considered a test
> >> comparing the LogLikelihoods, but have no idea what the degrees of
> >> freedom would be given that s(xc) is smoothed. I also tried:
> >>
> >> anova(g1$gam, g2$gam)
> >
> > gamm() fits via the lme() function of package nlme. To do what you want,
> > you need the anova() method for objects of class "lme", e.g.
> >
> > anova(g1$lme, g2$lme)
> >
> > Then I think you should check if the fits were done via REML and also be
> > aware of the issue of testing wether a variance term is 0.
> >
> >> that did not seem to return anything useful for this question.
> >>
> >> A related question is how to test the significance of adding a second
> >> random effect to a model that already has a random effect, such as:
> >>
> >> g3 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
> >> list(Case=~1, z=~1),correlation=corAR1())
> >> g4 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
> >> list(Case=~1, z=~1, int=~1),correlation=corAR1())
> >
> > Again, I think you need anova() on the $lme components.
> >
> > HTH
> >
> > G
> >
> >> Any help would be appreciated.
> >>
> >> Thanks.
> >>
> >> Will Shadish
> >> ********************************************
> >> g1
> >> $lme
> >> Linear mixedeffects model fit by maximum likelihood
> >> Data: data
> >> Loglikelihood: 437.696
> >> Fixed: fixed
> >> X(Intercept) Xz Xint Xs(xc)Fx1
> >> 0.6738466 2.5688317 0.0137415 0.1801294
> >>
> >> Random effects:
> >> Formula: ~Xr  1  g
> >> Structure: pdIdnot
> >> Xr1 Xr2 Xr3 Xr4
> >> Xr5 Xr6 Xr7 Xr8 Residual
> >> StdDev: 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781
> >> 0.0004377781 0.0004377781 0.0004377781 1.693177
> >>
> >> Correlation Structure: AR(1)
> >> Formula: ~1  g
> >> Parameter estimate(s):
> >> Phi
> >> 0.3110725
> >> Variance function:
> >> Structure: fixed weights
> >> Formula: ~invwt
> >> Number of Observations: 264
> >> Number of Groups: 1
> >>
> >> $gam
> >>
> >> Family: binomial
> >> Link function: logit
> >>
> >> Formula:
> >> y ~ s(xc) + z + int
> >>
> >> Estimated degrees of freedom:
> >> 1 total = 4
> >>
> >> attr(,"class")
> >> [1] "gamm" "list"
> >> ****************************
> >> > g2
> >> $lme
> >> Linear mixedeffects model fit by maximum likelihood
> >> Data: data
> >> Loglikelihood: 443.9495
> >> Fixed: fixed
> >> X(Intercept) Xz Xint Xs(xc)Fx1
> >> 0.720018143 2.562155820 0.003457463 0.045821030
> >>
> >> Random effects:
> >> Formula: ~Xr  1  g
> >> Structure: pdIdnot
> >> Xr1 Xr2 Xr3 Xr4
> >> Xr5 Xr6 Xr7 Xr8
> >> StdDev: 7.056078e06 7.056078e06 7.056078e06 7.056078e06 7.056078e06
> >> 7.056078e06 7.056078e06 7.056078e06
> >>
> >> Formula: ~1  xc %in% g
> >> (Intercept) Residual
> >> StdDev: 6.277279e05 1.683007
> >>
> >> Correlation Structure: AR(1)
> >> Formula: ~1  g/xc
> >> Parameter estimate(s):
> >> Phi
> >> 0.1809409
> >> Variance function:
> >> Structure: fixed weights
> >> Formula: ~invwt
> >> Number of Observations: 264
> >> Number of Groups:
> >> g xc %in% g
> >> 1 34
> >>
> >> $gam
> >>
> >> Family: binomial
> >> Link function: logit
> >>
> >> Formula:
> >> y ~ s(xc) + z + int
> >>
> >> Estimated degrees of freedom:
> >> 1 total = 4
> >>
> >> attr(,"class")
> >> [1] "gamm" "list"
> >>
> >>
> >
>

Gavin Simpson, PhD [t] +1 306 337 8863
Adjunct Professor, Department of Biology [f] +1 306 337 2410
Institute of Environmental Change & Society [e] [hidden email]
523 Research and Innovation Centre [tw] @ucfagls
University of Regina
Regina, SK S4S 0A2, Canada
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


I would be very nervous about relying on an anova call here. It will
attempt a generalized likelihood ratio test, but gamm is using penalized
quasi likelihood and there is really no likelihood here (even without
the problem that if there was a likelihood the null hypothesis would
still be on the edge of the feasible parameter space making the GLRT
problematic). The best hope might be to model the random effect of xc
using a term s(xc,bs="re") in the model formula (xc will need to be a
factor for this), and then use summary on the gam part of the fitted
model object to assess significance. If you do this you'll need to
include the grouping factor explicitly in corAR1 (at present it's picked
up from the random effect, so is nested in xc).
i.e
g2 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random =
list(xc=~1),correlation=corAR1())
becomes something like...
xf < factor(xc)
g2 < gamm(y ~ s(xc) +z+ int + s(xf,bs="re"),family=binomial,
weights=trial,
correlation=corAR1(form=~1xf))
summary(g2$gam)
... I'm also a bit nervous about xc entering as an iid random effect and
the argument of a smooth, however  does that model structure really
make sense?
best,
Simon
On 11/06/13 18:08, William Shadish wrote:
> Gavin et al.,
>
> Thanks so much for the help. Unfortunately, the command
>
> > anova(g1$lme, g2$lme)
>
> gives "Error in eval(expr, envir, enclos) : object 'fixed' not found
>
> and for bam (which is the one that can use a known ar1 term), the
> error is
>
> > AR1 parameter rho unused with generalized model
>
> Apparently it cannot run for binomial distributions, and presumably
> also Poisson.
>
> I did find a Frequently Asked Questions for package mgcv page that said
>
> "How can I compare gamm models? In the identity link normal errors
> case, then AIC and hypotheis testing based methods are fine. Otherwise
> it is best to work out a strategy based on the summary.gam"
>
> So putting all this together, I take it that my binomial example will
> not support a direct model comparison to test the significance of the
> random effects. I'm guessing the best strategy based on the
> summary.gam is probably just to compare fit indices like Log Likelihoods.
>
> If anyone has any other suggestions, though, please do let me know.
>
> Thanks so much.
>
> Will Shadish
>
> On 6/7/2013 3:02 PM, Gavin Simpson wrote:
>> On Fri, 20130607 at 13:12 0700, William Shadish wrote:
>>> Dear Rhelpers,
>>>
>>> I'd like to understand how to test the statistical significance of a
>>> random effect in gamm. I am using gamm because I want to test a model
>>> with an AR(1) error structure, and it is my understanding neither gam
>>> nor gamm4 will do the latter.
>>
>> gamm4() can't yes and out of the box mgcv::gam can't either but
>> see ?magic for an example of correlated errors and how the fits can be
>> manipulated to take the AR(1) (or any structure really as far as I can
>> tell) into account.
>>
>> You might like to look at mgcv::bam() which allows an known AR(1) term
>> but do check that it does what you think; with a random effect spline
>> I'm not at all certain that it will nest the AR(1) in the random effect
>> level.
>>
>> <snip />
>>> Consider, for example, two models, both with AR(1) but one allowing a
>>> random effect on xc:
>>>
>>> g1 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial,
>>> correlation=corAR1())
>>> g2 < gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random =
>>> list(xc=~1),correlation=corAR1())
>>
>> Shouldn't you specify how the AR(1) is nested in the hierarchy here,
>> i.e. AR(1) within xc? maybe I'm not following your data structure
>> correctly.
>>
>>> I include the output for g1 and g2 below, but the question is how to
>>> test the significance of the random effect on xc. I considered a test
>>> comparing the LogLikelihoods, but have no idea what the degrees of
>>> freedom would be given that s(xc) is smoothed. I also tried:
>>>
>>> anova(g1$gam, g2$gam)
>>
>> gamm() fits via the lme() function of package nlme. To do what you want,
>> you need the anova() method for objects of class "lme", e.g.
>>
>> anova(g1$lme, g2$lme)
>>
>> Then I think you should check if the fits were done via REML and also be
>> aware of the issue of testing wether a variance term is 0.
>>
>>> that did not seem to return anything useful for this question.
>>>
>>> A related question is how to test the significance of adding a second
>>> random effect to a model that already has a random effect, such as:
>>>
>>> g3 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
>>> list(Case=~1, z=~1),correlation=corAR1())
>>> g4 < gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random =
>>> list(Case=~1, z=~1, int=~1),correlation=corAR1())
>>
>> Again, I think you need anova() on the $lme components.
>>
>> HTH
>>
>> G
>>
>>> Any help would be appreciated.
>>>
>>> Thanks.
>>>
>>> Will Shadish
>>> ********************************************
>>> g1
>>> $lme
>>> Linear mixedeffects model fit by maximum likelihood
>>> Data: data
>>> Loglikelihood: 437.696
>>> Fixed: fixed
>>> X(Intercept) Xz Xint Xs(xc)Fx1
>>> 0.6738466 2.5688317 0.0137415 0.1801294
>>>
>>> Random effects:
>>> Formula: ~Xr  1  g
>>> Structure: pdIdnot
>>> Xr1 Xr2 Xr3 Xr4
>>> Xr5 Xr6 Xr7 Xr8 Residual
>>> StdDev: 0.0004377781 0.0004377781 0.0004377781 0.0004377781
>>> 0.0004377781
>>> 0.0004377781 0.0004377781 0.0004377781 1.693177
>>>
>>> Correlation Structure: AR(1)
>>> Formula: ~1  g
>>> Parameter estimate(s):
>>> Phi
>>> 0.3110725
>>> Variance function:
>>> Structure: fixed weights
>>> Formula: ~invwt
>>> Number of Observations: 264
>>> Number of Groups: 1
>>>
>>> $gam
>>>
>>> Family: binomial
>>> Link function: logit
>>>
>>> Formula:
>>> y ~ s(xc) + z + int
>>>
>>> Estimated degrees of freedom:
>>> 1 total = 4
>>>
>>> attr(,"class")
>>> [1] "gamm" "list"
>>> ****************************
>>> > g2
>>> $lme
>>> Linear mixedeffects model fit by maximum likelihood
>>> Data: data
>>> Loglikelihood: 443.9495
>>> Fixed: fixed
>>> X(Intercept) Xz Xint Xs(xc)Fx1
>>> 0.720018143 2.562155820 0.003457463 0.045821030
>>>
>>> Random effects:
>>> Formula: ~Xr  1  g
>>> Structure: pdIdnot
>>> Xr1 Xr2 Xr3 Xr4
>>> Xr5 Xr6 Xr7 Xr8
>>> StdDev: 7.056078e06 7.056078e06 7.056078e06 7.056078e06
>>> 7.056078e06
>>> 7.056078e06 7.056078e06 7.056078e06
>>>
>>> Formula: ~1  xc %in% g
>>> (Intercept) Residual
>>> StdDev: 6.277279e05 1.683007
>>>
>>> Correlation Structure: AR(1)
>>> Formula: ~1  g/xc
>>> Parameter estimate(s):
>>> Phi
>>> 0.1809409
>>> Variance function:
>>> Structure: fixed weights
>>> Formula: ~invwt
>>> Number of Observations: 264
>>> Number of Groups:
>>> g xc %in% g
>>> 1 34
>>>
>>> $gam
>>>
>>> Family: binomial
>>> Link function: logit
>>>
>>> Formula:
>>> y ~ s(xc) + z + int
>>>
>>> Estimated degrees of freedom:
>>> 1 total = 4
>>>
>>> attr(,"class")
>>> [1] "gamm" "list"
>>>
>>>
>>
>
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.

