[R-sig-ME] Testing for differences in random effects between two groups
Stephen T
stwebvanuatu at yahoo.com.au
Tue May 28 15:06:29 CEST 2013
Yes, I think this is what I want.
Here's an example:
> # homoscedastic model
> # I checked and these results match lmer
> model1=lme(MH11~SEX, random=~1|BIRD/NIGHT, data=calls)
> model1
Linear mixed-effects model fit by REML
Data: calls
Log-restricted-likelihood: -1825.803
Fixed: MH11 ~ SEX
(Intercept) SEXM
445.05914 94.77928
Random effects:
Formula: ~1 | BIRD
(Intercept)
StdDev: 32.91933
Formula: ~1 | NIGHT %in% BIRD
(Intercept) Residual
StdDev: 25.54875 31.08874
Number of Observations: 356
Number of Groups:
BIRD NIGHT %in% BIRD
57 138
> # heteroscedastic model
> model2=lme(MH11~SEX, random=~1|BIRD/NIGHT, data=calls, weights=varIdent(form=~1|SEX))
> model2
Linear mixed-effects model fit by REML
Data: calls
Log-restricted-likelihood: -1811.211
Fixed: MH11 ~ SEX
(Intercept) SEXM
445.06240 94.67614
Random effects:
Formula: ~1 | BIRD
(Intercept)
StdDev: 32.77675
Formula: ~1 | NIGHT %in% BIRD
(Intercept) Residual
StdDev: 22.54168 36.82237
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | SEX
Parameter estimates:
M F
1.0000000 0.5704416
Number of Observations: 356
Number of Groups:
BIRD NIGHT %in% BIRD
57 138
# likelihood ratio test
# I have read that one should actually use method="ML" in the model1 and model2 for this test
> anova(model1, model2)
Model df AIC BIC logLik Test L.Ratio p-value
model1 1 5 3661.606 3680.953 -1825.803
model2 2 6 3634.421 3657.637 -1811.210 1 vs 2 29.18513 <.0001
I will now have to read up on interpretation of nlme results.
The variance function results always sets SEXM = 1, which is unexpected because F is the first factor.
Anyhow, model2 indicated that SEXF sd was about 0.57 times SEXM, correct me if I'm wrong? I calculated a mean sd ratio = 0.58 directly from the data, ignoring NIGHT effects. Indeed, boxplots suggest that SEXM measurements were more variable than SEXF measurements.
More about what I am doing: model2 seems appropriate for my study. Looking at the box plots, I can see a shift between the groups SEXF and SEXM. Secondly, there is variation between- and within- individual BIRDs. Between-BIRD variation suggests individual vocal signatures. Some of the within-BIRD variation might be because different stimuli (call playbacks) were used on different NIGHTs to elicit responses, i.e. producing different excitedness in responses from the BIRDs. Thirdly, the groups SEXF and SEXM are heteroscedastic as noted above.
I hope to make some good progress now and maximise the value of my hard-earned data.
Thank you very much, Ewart.
----- Original Message -----
From: Ewart Thomas <ethomas at stanford.edu>
To: Stephen T <stwebvanuatu at yahoo.com.au>
Cc:
Sent: Tuesday, 28 May 2013 4:26 AM
Subject: Re: [R-sig-ME] Testing for differences in random effects between two groups
stephen, see if this site is useful: http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/7108
your code might look like:
library(nlme)
model1 = lme(MH11 ~ SEX, random = ~ 1 | BIRD / NIGHT, data = calls)
model2 = lme(MH11 ~ SEX, random = ~ 1 | BIRD / NIGHT, data = calls, weights = varIdent(form = ~ 1 | SEX))
# Test if equal var model is sig worse than unequal var model
anova(model1, model2)
good luck!
ewart
On May 27, 2013, at 6:48 AM, Stephen T wrote:
> Thanks,
>
> I hope that when I hit reply this message also goes to the mailing list. Please correct me if it doesn't.
>
>
> I am weak on model definition and don't understand the output from lmer here.
>
>
> The first model you suggested is equivalent to:
>> lmer(MH11~SEX+(1|SEX/BIRD/NIGHT), data=calls)
> Linear mixed model fit by REML
> Formula: MH11 ~ SEX + (1 | SEX/BIRD/NIGHT)
> Data: calls
> AIC BIC logLik deviance REMLdev
> 3664 3687 -1826 3664 3652
> Random effects:
> Groups Name Variance Std.Dev.
> NIGHT:(BIRD:SEX) (Intercept) 652.770 25.5494
> BIRD:SEX (Intercept) 1083.670 32.9191
> SEX (Intercept) 14.479 3.8052
> Residual 966.500 31.0886
> Number of obs: 356, groups: NIGHT:(BIRD:SEX), 138; BIRD:SEX, 57; SEX, 2
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 445.059 9.052 49.17
> SEXM 94.779 11.878 7.98
>
> Correlation of Fixed Effects:
> (Intr)
> SEXM -0.762
>
> The above results show the correct nesting structure. Your code actually returned NIGHT and not NIGHT:(BIRD:SEX) but the numbers are the same.
>
> The SEX random effect above is clearly not important but how to interpret the results?
> I was hoping to see something like this:
> NIGHT
> BIRD:SEXF
> BIRD:SEXM
> Residual
> Then I could do a likelihood ratio test of the above versus the simpler model. Or is the small SEX random effect variance itself the result I am seeking???
>
>
> The next model below has split the BIRD:SEX effect in two. It looks like a bit of a mess:
>> lmer(MH11~SEX+(1|SEX/BIRD)+(1|SEX/BIRD/NIGHT), data=calls)
> Linear mixed model fit by REML
> Formula: MH11 ~ SEX + (1 | SEX/BIRD) + (1 | SEX/BIRD/NIGHT)
> Data: calls
> AIC BIC logLik deviance REMLdev
> 3668 3699 -1826 3664 3652
> Random effects:
> Groups Name Variance Std.Dev.
> NIGHT:(BIRD:SEX) (Intercept) 652.762 25.5492
> BIRD:SEX (Intercept) 541.811 23.2768
> BIRD:SEX (Intercept) 541.867 23.2780
> SEX (Intercept) 14.473 3.8043
> SEX (Intercept) 14.475 3.8046
> Residual 966.502 31.0886
> Number of obs: 356, groups: NIGHT:(BIRD:SEX), 138; BIRD:SEX, 57; SEX, 2
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 445.059 9.803 45.40
> SEXM 94.779 13.016 7.28
>
> Correlation of Fixed Effects:
> (Intr)
> SEXM -0.753
>
>
> ----- Original Message -----
> From: Ewart Thomas <ethomas at stanford.edu>
> To: Stephen T <stwebvanuatu at yahoo.com.au>
> Cc:
> Sent: Monday, 27 May 2013 10:34 PM
> Subject: Re: [R-sig-ME] Testing for differences in random effects between two groups
>
> stephen, try inserting the nesting structure into the model:
>
> model1 = lmer(MH11 ~ SEX + (1 | SEX / BIRD) + (1 | NIGHT), data = calls),
>
> or even
>
> model2 = lmer(MH11 ~ SEX + (1 | SEX / BIRD) + (1 | SEX / BIRD / NIGHT), data = calls),
>
> i think model2 is the correct model. hope that works.
> ewart
>
> On May 27, 2013, at 5:13 AM, Stephen T wrote:
>
>> Hello,
>>
>> I am trying to apply linear models to some observational data.
>> The data are not balanced.
>>
>> Basically, I am wish to describe sexual and individual signatures in calls of a bird species. I have many call recordings and measured several acoustic variables from each recording.
>>
>> Here's an analysis for one variable:
>>
>>> lmer(MH11~SEX+(1|BIRD)+(1|NIGHT), data=calls)
>> Linear mixed model fit by REML
>> Formula: MH11 ~ SEX + (1 | BIRD) + (1 | NIGHT)
>> Data: calls
>> AIC BIC logLik deviance REMLdev
>> 3662 3681 -1826 3663 3652
>> Random effects:
>> Groups Name Variance Std.Dev.
>> NIGHT (Intercept) 652.77 25.549
>> BIRD (Intercept) 1083.67 32.919
>> Residual 966.50 31.089
>> Number of obs: 356, groups: NIGHT, 138; BIRD, 57
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 445.059 8.212 54.19
>> SEXM 94.779 10.588 8.95
>>
>> Correlation of Fixed Effects:
>> (Intr)
>> SEXM -0.776
>>
>> SEX is a fixed effect (Female/Male). BIRD (individuals) and NIGHT are random effects.
>> The nesting is SEX/BIRD/NIGHT/CALL.
>>
>> There is a clear difference between-SEXes (fixed effect), i.e. SEXF and SEX
>> M have different voices.
>>
>> Secondly, there is some interesting variation between-BIRDs. Individual BIRDs have different voices.
>>
>>
>> What if the BIRD random effect was not the same in both SEXF and SEXM groups? That would confound the model. Since there is a fixed effect (sexual signature), perhaps the random effect (individual signature) is different between groups SEXF and SEXM as well (unexpected, but should be tested).
>>
>> I'm not sure how to test for differences between groups in the random effects. One possibility is to split SEXF and SEXM and run separate lmer models. How would I the compare the results? Another is to compare the distributions of random effects between SEXF and SEXM (a two-sample test?).
>>
>> I would appreciate some advice on how to proceed?
>>
>> Stephen from Australia.
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list