I am a bit stuck with the use of linear mixed models and its random effects in a twin study analysis. What I have is microbiome data for twins that are from MZ twin pairs and DZ twin pairs, so no other subjects than complete twins.

I want to see differences between microbiome signatures in MZ vs DZ twin pairs and I was modelling following model for Shannon diversity:

alpha_divLMM <- lmer(Shannon~ Age+Gender+BMI+Antibiotics+VegetableIntake+FruitIntake+ Zygosity + ind +(1|pair), data=tw_sample_df)

where "ind" refers to either "individual1" or "individual2" in the twin-pair, and "pair" indicating which pair the twin is assigned to as a fixed effect.

However, I am not sure whether this is correct? I plotted "qqnorm" and "qqline" using this model of above, and it is not what it should look like: there is no heteroscedasticity or normality in the residuals. So where am I wrong? I found also in an article: https://www.nature.com/articles/s41598-022-07632-3#data-availability that they also include zygosity as random effect? they model this as follows: y ~ age+BMI+(1|shipmentNumber) + (0 +dummy(Zygosity, "MZ") | IndividualFamilyID) + (0+dummy(Zygosity, "DZ") |IndividualFamilyID) From their github: https://github.com/fnew/New_et_al_2021 , I think this individualFamily ID represents the identifier for the individual itself, but I'm not sure. How do I interpret this setup?

Any help is much appreciated!

thank you!

  • $\begingroup$ Welcome to Cross Validated. Does individual1 and individual2 indicate birth order? If not, why include ind in the model? And what do you mean by "it is not what it should look like"? You don't point to the article, so this is a guess: Identical twins are more similar (genetically) than fraternal twins. So the variance between two MZ twins is smaller than the variance between two DZ twins. This can be modeled as a zygosity random effect. $\endgroup$
    – dipetkov
    May 3, 2022 at 0:52
  • 1
    $\begingroup$ Please visit stats.stackexchange.com/help/merging-accounts to combine your accounts: that will enable you to edit this post without requiring users to vote on it. $\endgroup$
    – whuber
    May 7, 2022 at 13:46

1 Answer 1


I answer your question about the linear mixed model(s) for microbiome gene/species abundances by New et al. described in [1]. You already asked a very similar question here.

In the New et al. paper, IndividualFamilyID is a family ID variable: each family has a unique ID and each pair of twins share the same family ID. New et al. use IndividualFamilyID to capture the relatedness of twins.

Here is the New et al. mixed-effects model as a formula in lme4 syntax:

y ~ 1
  + Age.at.metagenomics.sample
  + IndividualBMI
  + (0 + dummy(twinstatus, "MZ") | IndividualFamilyID)
  + (0 + dummy(twinstatus, "DZ") | IndividualFamilyID)

where dummy(twinstatus, "MZ") is an indicator variable equal to 1 if the twins are monozygotic and 0 otherwise. Similary, dummy(twinstatus, "DZ") is an indicator variable equal to 1 if the twins are dizygotic. In this model, each family has a random intercept; the random intercepts for MZ families and for DZ families come from two different distributions.

To understand this better, here is the formal specification of the New et al. mixed-effects model:

\begin{aligned} Y_i &\sim N \left(\alpha_{\operatorname{Z}j[i]} + X\beta, \sigma^2 \right) \\ \alpha_{\operatorname{MZ}j} &\sim N \left(0, \sigma_{\operatorname{MZ}}^2 \right) \text{, for MZ pair j = 1,} \dots \text{,J}_{\operatorname{MZ}} \\ \alpha_{\operatorname{DZ}j} &\sim N \left(0, \sigma_{\operatorname{DZ}}^2 \right) \text{, for DZ pair j = 1,} \dots \text{,J}_{\operatorname{DZ}} \end{aligned}

Again, the twins in the same family are related because they share the same random intercept $\alpha_{\operatorname{Z}j[i]}$. Furthermore, monozygotic and dyzygotic twins have different degree of relatedness represented by $\sigma_{\operatorname{MZ}}^2$ and $\sigma_{\operatorname{DZ}}^2$ respectively. Naturally, we expect that identical twins are more similar than fraternal twins, ie, we expect $\sigma_{MZ}^2$ is higher than $\sigma_{DZ}^2$.

[1] New, F.N., Baer, B.R., Clark, A.G. et al. Collective effects of human genomic variation on microbiome function. Sci Rep 12, 3839 (2022). https://doi.org/10.1038/s41598-022-07632-3

Related CV posts about analysis of twins studies:

Linear mixed-effects modeling with twin study data
Linear mixed models nested random effects: can you nest within a factor that has just 2 observations?
Choice of referent twin in twin difference model


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.