[R-sig-ME] Understanding variance components

René bimonosom at gmail.com
Thu Jan 19 16:55:00 CET 2017


Hey Gang, hey Henrik (!) :)

very insightful query, hope it is okay that I briefly join it, on the last
question.
Variance is a result of division of squared deviations by sample size. In
the complete sample, this means divided by 240; and in the gender samples
by 120 each.
Have a look on the reversed equations for obtaining the overall variance by
combining gender groups, weighted by sample size:

(120 * 0.9944589   +   4.137316 * 120)/240 = overall variance (should be
equal to)  = (0.9944589 + 4.137316)/2   acording to your equation gang
left hand can be written as:
120 * (0.9944589   +   4.137316) / 240
which is:
(0.9944589   +   4.137316)/2

so nothing to worry, I think :)

Best wishes,
René



2017-01-19 16:14 GMT+01:00 Chen, Gang (NIH/NIMH) [C] <gangchen at mail.nih.gov>
:

> Nice, Henrik!
>
> One thing we need to resolve, though, is this.
>
> ====================
> The variances for model m1:
>
> (1) intercept
>
> summary(m1)$varcor$id[1]
> [1] 2.565887
>
> (2) residuals
>
> attr(summary(m1)$varcor, "sc")^2
> [1] 2.196212
>
> ====================
> And the variances for model m2:
>
> (1) female
>
> summary(m2)$varcor$id[1]
> [1] 0.9944589
>
> (2) male
>
> summary(m2)$varcor$id.1[1]
> [1] 4.137316
>
> (3) residuals
>
> attr(summary(m2)$varcor, "sc")^2
> [1] 2.196212
> ====================
>
> It’s great to see that the residual variance matches between the two
> models m1 and m2. However, the intercept variance in m1, 2.565887, is not
> equal to the sum of female and male variances, 0.9944589 + 4.137316 =
> 5.131775. However, if we divide the total variance of female and male by 2,
> we have (0.9944589 + 4.137316)/2 = 2.565887. Why is that?
>
> If we code the two groups as
>
> obk.long$gender_F <- sqrt(2)*as.numeric(obk.long$gender == "F")
> obk.long$gender_M <- sqrt(2)*as.numeric(obk.long$gender == "M”)
>
> then we have the desired result,
>
> m2 <- lmer(value ~ gender*phase+(0+gender_F|id)+(0+gender_M|id),
> data=obk.long)
> summary(m2)$varcor$id[1]+summary(m2)$varcor$id.1[1]
> [1] 2.565888
>
> Even though the variance part is reconciled, I cannot come up with a good
> explanation as to why this coding strategy is required. Any thought?
>
> Thanks,
> Gang
>
>
> On Jan 19, 2017, at 6:13 AM, Henrik Singmann <singmann at psychologie.uzh.ch<
> mailto:singmann at psychologie.uzh.ch>> wrote:
>
> Hi Gang,
>
> I have an idea which is based on the last example given on ?lmer:
> ## Fit sex-specific variances by constructing numeric dummy variables
> ...
>
> I am not sure if this is entirely correct, but it looks good to me. If
> not, hopefully someone more knowledgeable will jump in.
>
> ## Original model without gender specific variance:
> data(obk.long, package = "afex")
>
> m1 <- lmer(value ~ gender*phase+(1|id), data=obk.long)
> summary(m1)$varcor
> ## Groups   Name        Std.Dev.
> ## id       (Intercept) 1.6018
> ## Residual             1.4820
> REMLcrit(m1)
> ## [1] 911.1599
>
> ## to get gender specific vari8ances, we construct two dummy variables:
> obk.long$gender_F <- as.numeric(obk.long$gender == "F")
> obk.long$gender_M <- as.numeric(obk.long$gender == "M")
> m2 <- lmer(value ~ gender*phase+(0+gender_F|id)+(0+gender_M|id),
> data=obk.long)
> summary(m2)$varcor
> ## Groups   Name     Std.Dev.
> ## id       gender_F 0.99723
> ## id.1     gender_M 2.03404
> ## Residual          1.48196
> REMLcrit(m2)
> ## [1] 908.297
>
> So far, looks reasonably close. Same for the conditional modes (thx
> Phillip). Left two columns is separate, right is joint variance.
> cbind(ranef(m2)$id, rep(NA, 16), ranef(m1)$id)
> ##      gender_F   gender_M rep(NA, 16) (Intercept)
> ## 1   0.0000000 -3.0986753          NA  -3.0351426
> ## 2   0.0000000 -2.1328544          NA  -2.0891242
> ## 3   0.0000000  0.1207276          NA   0.1182523
> ## 4  -0.9806229  0.0000000          NA  -1.0642708
> ## 5  -0.3995130  0.0000000          NA  -0.4335918
> ## 6   0.0000000  2.6962499          NA   2.6409683
> ## 7   0.0000000  1.4084888          NA   1.3796103
> ## 8  -0.3995130  0.0000000          NA  -0.4335918
> ## 9  -0.6900680  0.0000000          NA  -0.7489313
> ## 10  0.0000000  0.4426679          NA   0.4335918
> ## 11  0.0000000 -1.1670336          NA  -1.1431057
> ## 12  0.0000000  1.7304291          NA   1.6949498
> ## 13  1.3438166  0.0000000          NA   1.4584452
> ## 14 -0.6900680  0.0000000          NA  -0.7489313
> ## 15  0.4721518  0.0000000          NA   0.5124267
> ## 16  1.3438166  0.0000000          NA   1.4584452
>
> And, finally, the same for the fixed effects:
> fixef(m1)
> ##       (Intercept)           genderM         phasepost
> ##      6.000000e+00      7.500000e-01     -6.250000e-01
> ##          phasepre genderM:phasepost  genderM:phasepre
> ##     -2.000000e+00     -1.243450e-15     -1.687539e-15
> fixef(m2)
> ##      (Intercept)           genderM         phasepost
> ##     6.000000e+00      7.500000e-01     -6.250000e-01
> ##         phasepre genderM:phasepost  genderM:phasepre
> ##    -2.000000e+00      1.829648e-14      1.820766e-14
>
> Hope that helps,
> Henrik
>
> On Jan 18, 2017, at 5:18 PM, Chen, Gang (NIH/NIMH) [C] <
> gangchen at mail.nih.gov<mailto:gangchen at mail.nih.gov>> wrote:
>
> Happy New Year, Henrik! Thanks for explaining the details. A couple of
> days after I posted the question, I realized that my question was silly!
> Once I laid out the LME model equation, my original confusion was resolved.
>
> Actually I meant to ask a slightly different question. Let me use the
> dataset embedded in your ‘afex’ package as an example:
>
> data(obk.long, package = "afex”)
>
> Suppose that my base model is
>
> lmer(value ~ gender*phase+(1|id), data=obk.long)
>
> Is there a way to specify a different variance for each gender in one
> model?
>
> Thanks,
> Gang
>
>
> On Jan 13, 2017, at 12:06 PM, Henrik Singmann <singmann at psychologie.uzh.ch
> <mailto:singmann at psychologie.uzh.ch>> wrote:
>
> Hi Gang,
>
> Sorry that I so am late to the party, but in case you are still interested
> I will reply (and, of course, for the archive).
>
> The answer is basically given in the old faq:
> http://glmm.wikidot.com/faq#toc27
>
> (1|site/block) = (1|site)+(1|site:block)
>
> Which is exactly what is given in your output. A random intercept for
> Worker and a random intercept for each worker:Machine interaction.
>
> To answer your questions. The random intercepts do not have base or
> reference levels. They are increments or decrements to the overall
> intercept for each level of Worker or the Machine:Worker combination. The
> reported variance is the estimated variance of these increments, which is
> most likely unequal to the actual variance you would obtain by calculating
> it from the estimated increments, which are sometimes called BLUPs (I
> wonder if a better term for those exist).
>
> Hope that helps,
> Henrik
>
> PS: Belated Happy New Year to everyone.
>
>
> Am 05.01.2017 um 17:28 schrieb Chen, Gang (NIH/NIMH) [C]:
> Suppose that I have the following dataset in R:
>
> library(lme4)
> data(Machines,package="nlme")
> mydata <- Machines[Machines$Machine!='C’,]
>
> With the following model:
>
> print(lmer(score ~ 1 + (1|Worker/Machine), data=mydata), ranef.comp="Var")
>
> I have the variance components as shown below:
>
> Random effects:
> Groups         Name        Variance
> Machine:Worker (Intercept) 46.00
> Worker         (Intercept) 13.84
> Residual                    1.16
>
> I have trouble understanding exactly what the first two components are:
> Machine:Worker and Worker? Specifically,
>
> 1) What is the variance for Worker: corresponding to the base (or
> reference) level of the factor Machine? If so, what is the base level: the
> first level in the dataset or alphabetically the first level (it happens to
> be the same in this particular dataset)?
>
> 2) What is the variance for Machine:Worker? Is it the variance for the
> second level of the factor Machine, or the extra variance relative to the
> variance for Worker?
>
> Furthermore, for the model:
>
> print(lmer(score ~ 1 + (1|Worker/Machine), data=Machines),
> ranef.comp="Var")
>
> what is the variance for Machine:Worker in the following result since
> there are 3 levels involved in the factor Machine?
>
> Random effects:
> Groups         Name        Variance
> Machine:Worker (Intercept) 60.2972
> Worker         (Intercept)  7.3959
> Residual                    0.9246
>
> Thanks,
> Gang
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> _______________________________________________
> R-sig-mixed-models at r-project.org<mailto:R-sig-mixed-models at r-project.org>
> mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]



More information about the R-sig-mixed-models mailing list