Random Effect in twin study of both MZ and DZ twin pairs

Question

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!

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