[R-sig-ME] subject-specific interpetation two random effects

[ben pelzer]

Dear Dimitris and others,

I would like to get back to the problem posted a while ago: different
estimates for correlated vs. uncorrelated intercepts of two groups, males
and females. The problem can be put more simply than I did before.

Let's assume a LINEAR random intercept null-model for the Y value of male i
in school j: Yij = b0 + uj + eij, with b0 being the grand mean (for males),
uj being the deviation of school j from the grand mean (for males), and eij
the individual error of male pupil i in school j. For females the model is
similar: Yij = a0 + vj + eij with grand mean a0, and vj  the deviation of
school j from the grand mean (for females).

I first estimated both models separately (with lme), leading to particular
estimates of the fixed effects a0 and b0.

mmale    <- lme(y ~ 1, random = ~ 1|school, data=damale,
control=lmeControl(opt="optim"))
mfemale  <- lme(y ~ 1, random = ~ 1|school, data=dafemale,
control=lmeControl(opt="optim"))

I then estimated both models simultaneously, in two ways: the first model
below has correlated intercepts (uj and vj), the second model has
independent intercepts:

mboth_corr <- lme(y ~ 0+male+female, random = ~ 0+male+female|school,
                      weights = varIdent(form= ~ 1 | male), data=da,
                      control=lmeControl(opt="optim"))

mboth_indep <- lme(y ~ 0+male+female, random = list(school = pdDiag(~
0+male+female)),
                       weights = varIdent(form= ~ 1 | male), data=da)

"male" and "female" are 0/1 indicators of gender.

For model "mboth_indep" the fixed estimates of "male" and "female" are
equal to the fixed intercept estimates of the separate models "mmale" and
"mfemale".
However, for model "mboth_corr, the fixed estimates  of "male" and "female"
are (in general) NOT equal to the fixed intercept estimates of the separate
models "mmale" and "mfemale". I find this hard to understand... Apparently,
allowing correlation between uj and vj leads to a better model (lower
deviance) but also to a different grand mean estimate for males as well as
for females.

What is more: if I simulate data in which all schools have the same (even)
number of pupils, half of which are male, half of which female, then using
"mboth_corr" with correlated intercepts does NOT lead to estimates
different from the ones in the separate models. Why is this the case?

I would really appreciate some explanation of these complex (for me)
results. Thanks in advance for any help!

Ben.





On Fri, 30 Dec 2022 at 15:38, ben pelzer <benpelzer using gmail.com> wrote:

> Dear Dimitris,
>
> Thanks for your response, that was really helpful and the estimates are
> now as good as equal!
>
> I thought the difference was caused by the subject-specific
> interpretation, which for the interaction model "mab" is different than for
> the two separate models, if I'm right ...
>
> After reading your comment, I also tried linear models (using lme) for
> males/females separately and for the two groups combined using interaction
> (for linear models the subject-specific interpretation would not make a
> difference). Again it appeared that allowing for correlation of the
> intercepts for males and females leads to slightly different estimates of
> the fixed effects, whereas no correlation leads to exactly the same fixed
> effects as in the separate analyses.
>
> Thanks for your help, and, of course, a happy and healthy 2023,
>
> Ben.
>
>
> On Thu, 29 Dec 2022 at 07:37, Dimitris Rizopoulos <
> d.rizopoulos using erasmusmc.nl> wrote:
>
>> The last model says that the random effects for males and females are
>> correlated-you would need to assume independent random effects.
>>
>> Best,
>> Dimitris
>>
>>
>> ——
>> Dimitris Rizopoulos
>> Professor of Biostatistics
>> Erasmus University Medical Center
>> The Netherlands
>> ------------------------------
>> *From:* R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on
>> behalf of ben pelzer <benpelzer using gmail.com>
>> *Sent:* Wednesday, December 28, 2022 3:13:26 PM
>> *To:* r-sig-mixed-models <r-sig-mixed-models using r-project.org>
>> *Subject:* [R-sig-ME] subject-specific interpetation two random effects
>>
>>
>>
>> Waarschuwing: Deze e-mail is afkomstig van buiten de organisatie. Klik
>> niet op links en open geen bijlagen, tenzij u de afzender herkent en weet
>> dat de inhoud veilig is.
>> Caution: This email originated from outside of the organization. Do not
>> click links or open attachments unless you recognize the sender and know
>> the content is safe.
>>
>>
>>
>> Dear list,
>>
>> Sorry for cross-posting this question; I got no response on Cross
>> Validated. Hope to get one here.
>>
>>
>> I have two groups of pupils, males and females in a number of schools. So,
>> two-level data, pupils nested in schools. For each gender group, a
>> logistic
>> model was estimated with a school-level predictor "type" of school. Next,
>> a
>> "complete" model was estimated for both gender groups simultaneously. Here
>> is the syntax I used:
>>
>>
>> library(GLMMadaptive)
>>
>>
>> ma  <- mixed_model(y ~ 1+type, random=~ 1|school, family=binomial,
>> data=maledata)
>>
>>
>>
>> mb  <- mixed_model(y ~ 1+type, random=~ 1|school, family=binomial,
>> data=femaledata)
>>
>>
>>
>> mab <- mixed_model(y ~ 0+male+female+type:male+type:female,
>>
>>                    random=~ 0+male+female|school, family=binomial,
>> data=completedata)
>>
>>
>> The final model mab uses interaction terms and is specified without
>> intercepts, so that its estimates can directly be compared with those of
>> the separate group models ma and mb.
>>
>>
>> I expected that the final "complete" model mab would give the same
>> estimates for the effects of "type" as the separate models do. More in
>> particular, that the regression coefficients of type:male and type:female
>> in model mab would be equal to the coefficients of "type" in ma and mb,
>> respectively.
>>
>>
>> However, this is not (exactly) the case: the coefficients are relatively
>> close, but clearly not equal. I first thought that these dissimilarities
>> might be due to the number of quadrature points used by mixed_model.
>> Hence,
>> I chose larger values for nAGQ but this did not help: the dissimilarities
>> persist, and all estimates hardly change. Also, no convergence problems
>> exist, everything seems all right.
>>
>>
>> Now my guess is that these differences are caused by the fact that the
>> regression coefficients are subject-specific. That is, for model ma, the
>> effect of "type" expresses the influence of schooltype for two schools
>> with
>> the same random school effect across male pupils. In contrast, the effect
>> of "type" in model mab expresses the effect of schooltype for two schools
>> with the same random school effect across male pupils and also across
>> female pupils. So in mab, there are two random school effects, one for
>> males and the other for females, which BOTH have to be equal. It toke me
>> quite a while to realise this, and still I'm not completely sure. I would
>> really appreciate it if someone could confirm my suspicion. Thanks!!!
>>
>>
>> Ben.
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>>
>> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=05%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C09a5f446bc1443dd190608dae8ddbf2a%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C638078336334964127%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=FAr0RrkoPQSniGE4JOjKQzQfjkohGTucLO0VC3R93oE%3D&reserved=0
>>
>

	[[alternative HTML version deleted]]