# Linear mixed-effects modeling with twin study data

Suppose I have some some response variable $y_{ij}$ that was measured from $j$th sibling in $i$th family. In addition, some behavioral data $x_{ij}$ were collected at the same time from each subject. I'm trying to analyze the situation with the following linear mixed-effects model:

$$y_{ij} = \alpha_0 + \alpha_1 x_{ij} + \delta_{1i} x_{ij} + \varepsilon_{ij}$$

where $\alpha_0$ and $\alpha_1$ are the fixed intercept and slope respectively, $\delta_{1i}$ is the random slope, and $\varepsilon_{ij}$ is the residual.

The assumptions for the random effects $\delta_{1i}$ and residual $\varepsilon_{ij}$ are (assuming there are only two siblings within each family)

\begin{align} \delta_{1i} &\stackrel{d}{\sim} N(0, \tau^2) \\[5pt] (\varepsilon_{i1}, \varepsilon_{i2})^T &\stackrel{d}{\sim} N((0, 0)^T, R) \end{align}

where $\tau^2$ is an unknown variance parameter and the variance-covariance structure $R$ is a 2 x 2 symmetric matrix of form

$$\begin{pmatrix} r_1^2&r_{12}^2\\ r_{12}^2&r_2^2 \end{pmatrix}$$

that models the correlation between the two siblings.

1. Is this an appropriate model for such a sibling study?

2. The data are a little bit complicated. Among the 50 families, close to 90% of them are dizygotic (DZ) twins. For the rest families,

1. two have only one sibling;
2. two have one DZ pair plus one sibling; and
3. two have one DZ pair plus two additional siblings.

I believe lme in the R package nlme can easily handle (1) with missing or unbalanced situation. My trouble is, how to deal with (2) and (3)? One possibility I can think of is to break each of those four families in (2) and (3) into two so that each subfamily would have one or two siblings so the above model could be still applied to. Is this fine? Another option would be to simply throw away the data from the extra one or two siblings in (2) and (3), which seems to be a waste. Any better approaches?

3. It seems that lme allows one to fix the $r$ values in the residual variance-covariance matrix $R$, for example $r_{12}^2$ = 0.5. Does it make sense to impose the correlation structure, or should I simply estimate it based on the data?

• What does $x_j$ denote? – Macro May 16 '12 at 18:39
• @Macro: Thanks for spotting that. Just modified the OP to indicate that $x_{ij}$ is an explanatory variable, behavioral measure from each sibling. – bluepole May 16 '12 at 18:49
• Very interesting question and application. I could be missing something but it seems to me this model is over-parameterized. The correlated errors $ϵ_{i1},ϵ_{i2}$ can effectively be factored into a an "unshared" component and a "shared" component, the latter of which has the same function as $δ_{0i}$. You'll have to either delete $\delta_{0i}$, make the $\epsilon$'s iid errors, or impose constraints like $r^2_{12} = .5$ for identifiability - are you doing that on purpose to decouple environmental/genetic components to sibling correlation? – Macro May 16 '12 at 23:56
• @Macro: You're right: $δ_{0i}$ is not necessary in the model. Thanks for pointing this out! Strangely lme does not complain about such redundancy. – bluepole May 17 '12 at 14:04
• Are you still working with this overparameterized model (that part of your question has not been edited)? – Macro May 18 '12 at 16:54

You can include twins and non-twins in a unified model by using a dummy variable and including random slopes in that dummy variable. Since all families have at most one set of twins, this will be relatively simple:

Let $A_{ij} = 1$ if sibling $j$ in family $i$ is a twin, and 0 otherwise. I'm assuming you also want the random slope to differ for twins vs. regular siblings - if not, do not include the $\eta_{i3}$ term in the model below.

Then fit the model:

$$y_{ij} = \alpha_{0} + \alpha_{1} x_{ij} + \eta_{i0} + \eta_{i1} A_{ij} + \eta_{i2} x_{ij} + \eta_{i3} x_{ij} A_{ij} + \varepsilon_{ij}$$

• $\alpha_{0}, \alpha_{1}$ are fixed effect, as in your specifiation

• $\eta_{i0}$ is the 'baseline' sibling random effect and $\eta_{i1}$ is the additional random effect that allows twins to be more similar than regular siblings. The sizes of the corresponding random effect variances quantify how similar siblings are and how much more similar twins are than regular siblings. Note that both twin and non-twin correlations are characterized by this model - twin correlations are calculated by summing random effects appropriately (plug in $A_{ij}=1$).

• $\eta_{i2}$ and $\eta_{i3}$ have analogous roles, only they act as the random slopes of $x_{ij}$

• $\varepsilon_{ij}$ are iid error terms - note that I have written your model slightly differently in terms of random intercepts rather than correlated residual errors.

You can fit the model using the R package lme4. In the code below the dependent variable is y, the dummy variable is A, the predictor is x, the product of the dummy variable and the predictor is Ax and famID is the identifier number for the family. Your data is assumed to be stored in a data frame D, with these variables as columns.

library(lme4)
g <- lmer(y ~ x + (1+A+x+Ax|famID), data=D)


The random effect variables and the fixed effects estimates can be viewed by typing summary(g). Note that this model allows the random effects to be freely correlated with each other.

In many cases, it may make more sense (or be more easily interpretable) to assume independence between the random effects (e.g. this assumption is often made to decompose genetic vs. environmental familial correlation), in which case you'd instead type

g <- lmer(y ~ x + (1|famID) + (A-1|famID) + (x-1|famID) +(Ax-1|famID), data=D)

• This is really a nice solution, and I like it! Will try it out soon, and see it goes... Thanks a lot! – bluepole May 20 '12 at 10:25
• You're welcome. If you've found this solution helpful please consider accepting the answer :) – Macro May 20 '12 at 14:26
• Two issues: 1) Since most subjects are dizygotic twins, your approach seems not modeling the correlation between a DZ twin pair. 2) Only 4 families have extra siblings. I'm worried it would be hard to estimate the random effects for the siblings based on only those 4 families. Because the difference between a DZ twin pair and another sibling is relatively small (mainly environmental, not genetic), maybe I can simply ignore the subtle difference of twin vs. sibling, and treat those few siblings as twins with random effects as in your model or with correlated residuals as in my OP. – bluepole May 20 '12 at 19:40
• This approach does model the correlation between twins. For example, if their predictor values are 0, then the correlation between twins is $$\frac{ \sigma_{0}^{2} + \sigma_{1}^{2} }{ \sigma_{0}^{2} + \sigma_{1}^{2} + \sigma^{2}_{\varepsilon}}$$ where $\sigma_{0}^{2}, \sigma_{1}^{2}$ are the variances of $\eta_{i0}, \eta_{i1}$, respectively and $\sigma^{2}_{\varepsilon}$ is the variance of the error term. When the predictor values are not 0, this expression will also involve the variances of the other two random effects. – Macro May 21 '12 at 0:39
• You're right that, since there are few non-twins, the variances of $\eta_{i0}$ and $\eta_{i2}$ are going to be difficult to estimate. You can leave them out, but you don't lose anything by using the model I've suggested but possibly computational brevity. If you do, you're effectively assuming that non-twin siblings are independent. But you can still be using those observations to estimate the mean parameters (i.e. don't leave them out of the model fitting). Or, as you said, you can just act as though regular siblings are the same as twins, and you wouldn't need to dummy coding at all. – Macro May 21 '12 at 0:41