Multiple Random-Effect Components in Genetic Models

Why Use More Than One Random-Effect Component?

Individuals can resemble one another for more than one reason. Genetic similarity may induce covariance between their phenotypes, while sharing a household may contribute additional covariance. A model with only one genetic covariance component can attribute too much of the observed resemblance to genetics if household structure is relevant. A model with multiple components represents these sources separately, provided the available covariance patterns contain enough information to distinguish them.

The key idea is to model the covariance between every pair of observations as the sum of contributions from specified sources. A variance parameter controls the size of a contribution; a covariance kernel specifies how that contribution is distributed across individuals. The model does not automatically establish that a component has a causal interpretation. That interpretation depends on the design, the kernels, and the assumptions used to construct them.

Model and Covariance Structure

Definition: A multiple-component linear mixed model for $n$ individuals is $$ y=X\beta+g+h+e, $$ where $y$ is the $n$-vector of phenotypes, $X$ is a known fixed-effect design matrix, $\beta$ is a vector of fixed-effect coefficients, $g$ is a genetic random-effect vector, $h$ is a household random-effect vector, and $e$ is the residual vector. A standard covariance specification is $$ g\sim N(0,\sigma_g^2K_g),\qquad h\sim N(0,\sigma_h^2K_h),\qquad e\sim N(0,\sigma_e^2I), $$ with $g,h,e$ mutually independent. Here $K_g$ and $K_h$ are specified positive semidefinite covariance kernels, and $I$ is the $n\times n$ identity matrix. Consequently, $$ y\sim N(X\beta,V),\qquad V=\sigma_g^2K_g+\sigma_h^2K_h+\sigma_e^2I. $$

A kernel describes a pattern, not the amount of variance by itself. For example, entries of $K_g$ encode the relative genetic similarity of pairs of individuals, while entries of $K_h$ encode household sharing. If both kernels have diagonal entries equal to 1, then $\sigma_g^2$ and $\sigma_h^2$ are the corresponding marginal variance contributions for each individual. If a kernel is scaled differently, its variance parameter changes scale accordingly. Therefore, component estimates should be interpreted together with the definitions and scaling of their kernels.

The diagonal entries of $V$ are marginal variances, and its off-diagonal entries are covariances between distinct individuals. Under the model above, individuals $i$ and $j$ have covariance $$ \operatorname{Cov}(y_i,y_j)=\sigma_g^2(K_g)_{ij}+\sigma_h^2(K_h)_{ij} $$ for $i\ne j$; the independent residual term contributes only to the diagonal. If the random effects are not mutually independent, cross-covariance terms may also be required, so the displayed sum is not appropriate without the independence assumption.

Worked Example: Combining Genetic and Household Kernels

Suppose three hypothetical individuals have the following genetic and household kernels:

K_g= \begin{pmatrix} 1&0.4&0.1\\ 0.4&1&0.1\\ 0.1&0.1&1 \end{pmatrix}, \qquad K_h= \begin{pmatrix} 1&1&0\\ 1&1&0\\ 0&0&1 \end{pmatrix}.

For illustration, take $\sigma_g^2=2$, $\sigma_h^2=1$, and $\sigma_e^2=3$, in squared phenotype units. The marginal covariance matrix is obtained by adding the three scaled components:

V=2K_g+K_h+3I = \begin{pmatrix} 2&0.8&0.2\\ 0.8&2&0.2\\ 0.2&0.2&2 \end{pmatrix} + \begin{pmatrix} 1&1&0\\ 1&1&0\\ 0&0&1 \end{pmatrix} + \begin{pmatrix} 3&0&0\\ 0&3&0\\ 0&0&3 \end{pmatrix} = \begin{pmatrix} 6&1.8&0.2\\ 1.8&6&0.2\\ 0.2&0.2&6 \end{pmatrix}.

All three diagonal entries are $6$ squared phenotype units. Individuals 1 and 2 have covariance $2(0.4)+1(1)=1.8$: both genetic similarity and shared household contribute. Individuals 1 and 3 have covariance $2(0.1)+1(0)=0.2$, entirely from the genetic kernel in this specification. These are model-implied covariances, not evidence by themselves that either source causes the observed resemblance.

Joint Estimation of the Variance Parameters

The variance parameters are estimated together because they all determine the same marginal covariance matrix $V$. For instance, a large covariance between two household members could be explained partly by their genetic similarity and partly by their shared household. The likelihood compares the overall covariance pattern predicted by all components with the pattern in the observed phenotypes; it does not generally estimate each source in isolation.

For REML, assume the model is Gaussian, $X$ has full column rank $p$, and $V$ is positive definite for the parameter values considered. Define

P=V^{-1}-V^{-1}X(X'V^{-1}X)^{-1}X'V^{-1}.

The restricted log likelihood, up to the displayed normalizing constant, is

\ell_R=-\frac12\left[(n-p)\log(2\pi)+\log|V|+\log|X'V^{-1}X|+y'Py\right].

Here $P$ accounts for estimating the fixed effects, so REML evaluates the covariance information in contrasts that remove the fixed-effect space. With $\theta_g=\sigma_g^2$, $\theta_h=\sigma_h^2$, and $\theta_e=\sigma_e^2$, the model has $V=\theta_gK_g+\theta_hK_h+\theta_eI$. REML chooses these parameters jointly to maximize $\ell_R$, subject to the variance constraints $\theta_g,\theta_h,\theta_e\ge 0$. Maximum likelihood is another option; it uses the full likelihood rather than the restricted likelihood.

Why the REML Score Combines the Components

A short derivative shows explicitly why the estimates depend on the kernels together. If $\theta_j$ is one variance parameter and $K_j$ its kernel, then $\partial V/\partial\theta_j=K_j$. Differentiating the restricted log likelihood gives

\frac{\partial\ell_R}{\partial\theta_j} =\frac12\left[y'PK_jPy-\operatorname{tr}(PK_j)\right].

To see the form of this expression, the differential of the log determinant terms contributes $-\tfrac12\operatorname{tr}(P\,dV)$, while differentiating the quadratic term contributes $\tfrac12 y'P(dV)Py$. Substituting $dV=K_j\,d\theta_j$ and dividing by $d\theta_j$ yields the score above. Each score depends on $P$, which is constructed from the full $V$, not from one component alone. Thus changing one variance parameter changes the fitted covariance structure used to assess the others. At a boundary, such as $\theta_j=0$, a maximum can occur without the unconstrained score being zero.

Worked Example: Illustrating Covariance-Based Separation

Consider a hypothetical, centered cohort with unit-diagonal kernels. Suppose the average cross-product for pairs who are genetically related but do not share a household is $0.30$ phenotype-units squared, and their genetic-kernel entry is $0.5$. Suppose the average cross-product for unrelated pairs sharing a household is $0.45$ phenotype-units squared. For this illustration, take their household-kernel entry as $1$; the other kernel entry is $0$ for each pair type. Finally, suppose the average marginal variance is $2.00$ phenotype-units squared.

Under the model, the expected cross-product for the first pair type is $0.5\sigma_g^2$. Equating it to the stated average gives the illustrative moment estimate

\widehat{\sigma_g^2}=\frac{0.30}{0.5}=0.600.

For the unrelated household-sharing pairs, the expected cross-product is $\sigma_h^2$, so

\widehat{\sigma_h^2}=\frac{0.45}{1}=0.450.

Since both kernels have unit diagonal, the model's marginal variance is $\sigma_g^2+\sigma_h^2+\sigma_e^2$. The residual component estimate from the stated marginal variance is therefore

\widehat{\sigma_e^2}=2.00-0.600-0.450=0.950.

These arithmetic steps illustrate how distinct pair types can supply information about distinct covariance patterns. They are simple method-of-moments calculations, not a claim that REML uses only these averages. In an actual analysis, the estimates and their uncertainty depend jointly on all observations, fixed effects, and kernel entries. The example's values are internally consistent and nonnegative, but finite-sample moment estimates need not always be.

Identifiability: Can the Components Be Distinguished?

Separate estimates require covariance structures that are sufficiently distinct in the observed design. A simple failure occurs when two kernels are identical. If $K_g=K_h=K$, then

V=(\sigma_g^2+\sigma_h^2)K+\sigma_e^2I.

The likelihood can then contain information about the sum $\sigma_g^2+\sigma_h^2$, but not about how that sum is divided between the genetic and household components. More generally, identical or nearly proportional kernels make separate components unidentifiable or weakly identified. Having two named sources in a model does not guarantee that the data can estimate their contributions separately.

Worked Example: Identical Kernels Do Not Identify Two Components

Suppose a hypothetical analysis uses $K_g=K_h=K$, where $K$ is a specified kernel, and fixes the residual variance at $0.5$ squared phenotype units. Compare the parameter pairs $(\sigma_g^2,\sigma_h^2)=(1.2,0.8)$ and $(1.5,0.5)$. Both pairs have the same sum:

1.2+0.8=2.0,\qquad 1.5+0.5=2.0.

For either pair, the marginal covariance is

$$ V=2.0K+0.5I. $$

The two choices therefore assign different values to the genetic and household parameters but give exactly the same $V$. Under this model they produce the same likelihood for every possible $y$ and $X$. No estimation procedure using this covariance model can determine which division is correct. The remedy is not simply to collect more observations under the same indistinguishable covariance design; separation requires informative differences between the kernels or other justified design information.

Interpretation and Practical Checks

A fitted variance parameter describes the scale of its model component, conditional on the chosen kernels, fixed effects, and independence assumptions. With unit-diagonal kernels, for example, $\widehat{\sigma_g^2}$ is the fitted diagonal contribution associated with the genetic kernel. Its off-diagonal contribution for a pair is $\widehat{\sigma_g^2}(K_g)_{ij}$, so the parameter alone does not give the covariance for every pair. Changing kernel scaling changes the numerical value of the corresponding parameter while leaving the intended covariance model potentially unchanged.

Before interpreting component estimates, check that each kernel is defined on the same individuals and that its scaling and entries have the intended meaning. Verify that the covariance matrix is valid, that the random-effect independence assumption is defensible, and that the design provides contrasts capable of distinguishing the kernels. The nonnegative constraints matter: a component may be estimated at zero, which indicates a boundary estimate under the fitted model, not proof that the underlying biological source is absent.

A common mistake is to read each estimate as though it were calculated from a separate set of observations: “the household estimate is the covariance among household members, and the genetic estimate is the covariance among relatives.” This is generally wrong because pair types can share both kernels, and every component affects the joint likelihood through $V$. The correction is to interpret each estimate as part of a jointly fitted covariance decomposition and examine which pairwise covariance patterns the full model implies.

Key takeaway: Multiple random-effect components add scaled kernels to form one marginal covariance matrix. Their variance parameters are estimated jointly, and separate interpretation requires the kernels to encode distinguishable patterns in the data.

Check Your Understanding

Use the model and assumptions above to answer these questions.

If the diagonal entries of $K_g$ equal 1, what is the genetic component's contribution to an individual's marginal variance?
For two individuals with genetic-kernel entry $0.25$ and household-kernel entry $1$, write their model-implied covariance in terms of $\sigma_g^2$ and $\sigma_h^2$.
Why are the variance parameters generally estimated jointly rather than one at a time?
If $K_g=K_h$, which aspect of the two variance parameters can the covariance model identify?
What does a nonnegative constraint on a variance parameter imply about an estimate at zero, and what does it not establish biologically?