Epistatic Relationship Matrices For Genome-Wide Models

Why Build an Epistatic Relationship Matrix?

A genome-wide model can represent additive genetic resemblance with a relationship matrix whose entries summarize how similar two individuals are across markers. To represent additive-by-additive epistasis, the model needs a different kernel: one that summarizes similarity in products of additive marker scores. The matrix is a way to encode that resemblance across all individuals without estimating every pairwise interaction as a separate fixed effect.

The construction has an important detail. Squaring an additive relationship matrix entry by entry produces a useful all-pairs product kernel, but it includes products involving the same marker as well as products involving distinct markers. When the scientific target is strictly interaction between two different loci, the corresponding matrix can be built directly from distinct-locus score products.

The scale of \(K_A\) depends on how marker scores are coded and scaled. A common choice under Hardy–Weinberg equilibrium is to center a genotype count by its expectation \(2p_j\) and divide by \(\sqrt{2p_j(1-p_j)}\), where \(p_j\) is the allele frequency used for marker \(j\). With that coding, the score has variance one under the assumed genotype distribution. In practice, allele frequencies, missing-data handling, and any further matrix scaling should be specified, because they affect the numerical scale of the variance component.

From Marker Relationships to an Additive-by-Additive Kernel

For each pair of individuals \(i\) and \(l\), the squared additive relationship is

This expansion shows why the Hadamard square represents product-score resemblance. It compares the product of an individual's scores at markers \(j\) and \(k\) with the corresponding product for another individual. The sum includes all ordered pairs \((j,k)\), including \(j=k\). For \(j\ne k\), reversing the order gives the same product, so each distinct unordered pair appears twice.

A conventional additive-by-additive kernel is therefore \(K_{AA}=K_A\circ K_A\). Its definition is precise and convenient, but its interpretation is an all-pairs product kernel, not exclusively a sum over distinct two-locus interactions. If the target is only interactions between different loci, construct the matrix from those pair features directly.

These two constructions are related, but they need not be numerically identical. In particular, \(K_A\circ K_A\) contains the terms \(z_{ij}^{2}z_{lj}^{2}\) from \(j=k\), whereas \(K_{AA}^{\mathrm{distinct}}\) omits them. A model should state which kernel it uses rather than treating the labels “additive-by-additive” and “Hadamard square” as interchangeable in every interpretation.

Why the Kernel Is a Valid Covariance Structure

A relationship matrix used as a covariance matrix must be positive semidefinite: for every vector \(a\), \(a^\mathsf{T}Ka\ge 0\). The additive matrix has this property because it is a scaled Gram matrix. Specifically, \(a^\mathsf{T}K_Aa=(Z^\mathsf{T}a)^\mathsf{T}(Z^\mathsf{T}a)/m\ge 0\).

The Hadamard square is also positive semidefinite. One justification is to write each row of \(Z\) as a vector \(z_i\). The entry \(K_{A,il}\) is the inner product \(z_i^\mathsf{T}z_l/m\), so its square is the inner product between the tensor-product feature vectors \(z_i\otimes z_i\) and \(z_l\otimes z_l\), divided by \(m^2\). Thus \(K_A\circ K_A\) is itself a Gram matrix. Equivalently, the Schur product theorem states that the element-by-element product of two positive semidefinite matrices is positive semidefinite. The distinct-locus matrix is also a Gram matrix, since it is proportional to \(VV^\mathsf{T}\).

Variance-Component Interpretation

A genome-wide additive-by-additive relationship matrix is commonly used in a linear mixed model. For a quantitative phenotype vector \(y\), one possible model is

Here \(X\) is the fixed-effect design matrix, \(\beta\) is its coefficient vector, \(I\) is the identity matrix, and the random terms are assumed mutually independent. \(K_{AA}\) must be defined explicitly as either the all-pairs Hadamard-square kernel or a distinct-locus kernel. Under the stated model, the phenotype covariance is

The parameter \(\sigma_{AA}^2\) is the variance scale attached to the chosen kernel in the modeled population. It is not automatically a scale-free percentage of genetic variance: its numerical value depends on the kernel's normalization, marker coding, and model specification. A proportion such as \(\sigma_{AA}^2/(\sigma_A^2+\sigma_{AA}^2+\sigma_e^2)\) has an interpretation as a variance fraction only when the component scales and the target population make those terms commensurate.

One way to see the connection with interaction effects is to let \(q_{jk}\) be the vector across individuals whose \(i\)th element is \(z_{ij}z_{ik}\). If distinct-pair effects are independent with mean zero and common variance \(\tau^2\), then \(g_{AA}=\sum_{j<k}q_{jk}\gamma_{jk}\) has covariance \(\tau^2\sum_{j<k}q_{jk}q_{jk}^\mathsf{T}\). With \(K_{AA}^{\mathrm{distinct}}=VV^\mathsf{T}/\binom{m}{2}\), this is \(\sigma_{AA}^2K_{AA}^{\mathrm{distinct}}\) when \(\sigma_{AA}^2=\binom{m}{2}\tau^2\). The relationship-matrix variance component therefore aggregates the assumed variance of many interaction effects.

Constructing a Distinct-Locus Matrix

The next example makes the distinction between all-pairs and distinct-locus kernels concrete. For three markers there are \(\binom{3}{2}=3\) unordered distinct pairs. Compute one product feature per pair and individual, then take the normalized Gram matrix.

Assumptions for Genome-Wide Estimation

The matrix construction is algebraic; interpreting an estimated component as genome-wide epistatic variance requires additional assumptions. The model should state how genotype scores and allele frequencies were defined, how missing genotypes were handled, and how each kernel was normalized. Marker scores should represent the intended population and be coded consistently across individuals. Poorly handled missingness or frequency differences can change the relationship structure rather than merely reduce precision.

The mixed-model interpretation also depends on the random-effects assumptions. Effects represented by a kernel are modeled as mean-zero random quantities with a covariance proportional to that kernel. The residual covariance must be specified, and fixed effects in \(X\) should account for relevant design variables and population structure. The components are most straightforward to distinguish when their covariance patterns are sufficiently different in the sample. If kernels are strongly correlated, the data may provide limited information for allocating covariance uniquely among them. Linkage disequilibrium can affect that separation and the interpretation of component estimates; the kernel alone does not remove that issue.

Genome-wide estimates are estimates for the modeled marker set and population under these assumptions. They do not imply that every marker pair has an effect, that all interactions share the same biological mechanism, or that the estimated component will transport unchanged to another population. A variance-component estimate is a summary of covariance captured by a specified model, not a list of identified causal interactions.

Common Mistakes

• Using a wrong dot product. The entry \(K_{A,il}\) is the dot product of rows \(i\) and \(l\), divided by \(m\). For score rows \((1,0)\) and \((0,1)\), that dot product is zero, not one.
• Confusing elementwise and ordinary matrix squaring. \(K_A\circ K_A\) squares each entry. It is not the matrix product \(K_AK_A\).
• Calling every Hadamard square a strictly distinct-locus kernel. Its expansion includes \(j=k\) terms. Use a product-feature matrix indexed by \(j<k\) when distinct loci are required.
• Ignoring normalization. Rescaling a relationship matrix changes the numerical variance component attached to it. Compare variance estimates only after checking the kernel definition and scale.
• Interpreting covariance as causation. A positive kernel component describes covariance under a statistical model; it does not establish a biological mechanism or identify causal marker pairs.