In your solution, the model takes care only if the samples belong to the same family. I would like to know how we can add the level of relatedness such as MZ twin (100%) and DZ twin (50%). Currently, I have a variable zygosity with 2 levels (MZ and DZ).
The only real way to do this in a generalized linear model is to add a column to the model matrix denoting the twin type. In cases where you have only a single set of twins per family, then there's no gain from doing this (in fact, the model matrix will be rank deficient). This would only work for cases where you have more individuals per family than simply the twins. If this approach isn't acceptable, then the only alternative is to use a different kind of model, where you can input a sample correlation or relatedness structure of some sort. Limma provides some methods in that regard, though I don't know enough about them to say how much more useful they may be in this case.