[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