Random Effect in twin study of both MZ and DZ twin pairs 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!

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 
