Lecture 7 of the Quantitative Genetics series. The multivariate normal distribution, multi-trait models, genetic and environmental correlations, and GxE.
After this lecture you will be able to:
What is covered: the multivariate normal distribution, the two-trait model, genetic and environmental correlation, and genotype-environment interaction.
Lecture 1 explained why that bell shape arises, many small independent gene effects summing under the Central Limit Theorem. This lecture takes the next step: what happens when we consider several traits, or several animals, at once.
A reminder of the one-dimensional case, written \(X\sim N\left(\mu,\sigma^{2}\right)\):
with \(\mathbb{E}(X)=\mu\) and \(\mathbb{V}(X)=\sigma^{2}\). The distribution is completely specified by two numbers, its mean and its variance. That economy is what makes it so useful, and it is exactly what generalises next.
Now consider \(k\) variables together. If all are normally distributed marginally, and jointly they form a multivariate normal distribution, we write \(\boldsymbol{X}\sim N\left(\boldsymbol{\mu},\boldsymbol{\Sigma}\right)\) with density
The two parameters are now a vector of means and a matrix of variances and covariances:
Covariance measures association, but its units are awkward, so we usually standardise it:
which lets the covariance matrix be written entirely in terms of correlations and standard deviations. If two traits are independent they are uncorrelated, so \(\sigma_{XZ}=\rho_{XZ}\sigma_X\sigma_Z=0\) and the corresponding off-diagonal entries vanish. A diagonal \(\boldsymbol{\Sigma}\) therefore means mutually independent traits.
In other words, single-trait theory is the exception in practice. As soon as you have an index over several traits, or you want to predict a breeding value from relatives' records, you are working inside a multivariate normal distribution whether you name it or not.
The single most useful property of the multivariate normal is that its conditional distributions are also normal, with a mean that is a linear function of what you condition on. Partition the vector, mean and covariance matrix as
The conditional mean is the multivariate generalisation of the regression formula from Lecture 2, and the conditional variance is smaller than the unconditional one: observing \(\boldsymbol{X}_2\) removes uncertainty about \(\boldsymbol{X}_1\).
This is where the machinery pays off. Everything below is one application of the conditional formula.
Consider the two variables \(g\) (breeding value) and \(y\) (phenotype) for a single animal:
Suppose \(y=10\), \(\mu=8\), \(\sigma_y^2=4\) and \(h^2=0.25\).
Apply the conditional mean with \(\Sigma_{12}=h^2\sigma_y^2=1\) and \(\Sigma_{22}=\sigma_y^2=4\):
The regression coefficient is exactly \(h^2\), which is the result met in Lecture 6 when deriving the breeder's equation. The conditional variance, \(\sigma_g^2-\left(h^2\sigma_y^2\right)^2/\sigma_y^2=h^2\sigma_y^2(1-h^2)=0.75\), quantifies how much uncertainty remains.
Now stack two related animals. With relationship \(F_{ij}\) between them and independent environments, the joint distribution of both breeding values and both phenotypes is
Conditioning the top block (the breeding values) on the bottom block (the observed phenotypes) predicts both animals' genetic values from both records, weighting each by relationship. Notice that \(F_{ij}\), from the relationship lecture, is what couples the two animals. This is BLUP in miniature: the numerator relationship matrix \(\mathbf{A}\) simply extends this \(2\times 2\) block to the whole pedigree.
Real examples of trait relationships:
Starting from the general model \(y=\mu+\sum_{j=1}^{L}g_j+e\), write it once per trait:
Where does a phenotypic correlation between two traits come from? Decompose it.
There are therefore three correlations between any two traits: phenotypic \(\rho_{12}\), genetic \(\rho_{g_{12}}\) and environmental \(\rho_{e_{12}}\).
The phenotypic correlation you can measure directly is a blend of the genetic and environmental correlations, weighted by the heritabilities. The two components can even have opposite signs, so a near-zero phenotypic correlation does not mean the traits are genetically unrelated. Only the genetic correlation \(\rho_{g_{12}}\) predicts correlated response to selection, which is why the correlated-response formula in Lecture 6 contains \(r_g\) and not \(\rho_{12}\).
The decomposition above rested on the assumption that \(\mathbb{C}(g_1,e_2)=\mathbb{C}(g_2,e_1)=0\). It is worth asking what happens when that fails, and the two ways it can fail are different phenomena.
Genotype-environment correlation occurs when individuals with high genetic values for a trait also experience environments with high values for that trait, that is, when certain genotypes "select" certain environments.
When it is present the cross terms no longer vanish and additional covariance components enter the expression for \(\mathbb{C}(y_1,y_2)\). Practically, it inflates apparent genetic effects: the best-bred animals may also be the best-fed ones, and the phenotype cannot tell the difference.
The distinction matters. A genotype-environment correlation means genotypes are non-randomly distributed across environments; an interaction means their ranking changes between environments. Interaction is the reason a sire that is elite in a temperate intensive system may be mediocre in a tropical extensive one, and it is why ASAP-Bio's concern with locally relevant evaluation is a genetic argument, not merely a political one.
Covariances from correlations. Consider three variables with
1. Since \(\sigma_{ij}=\rho_{ij}\sigma_i\sigma_j\) and all standard deviations are 1, the covariances equal the correlations numerically: \(\sigma_{12}=0.3\), \(\sigma_{13}=-0.9\), \(\sigma_{23}=0.0\). The covariance matrix is
2. In R:
install.packages("mvtnorm"); library(mvtnorm)
mu <- c(-2, 3, 5)
sigma <- matrix(c(1,0.3,-0.9, 0.3,1,0, -0.9,0,1), 3, 3)
x <- rmvnorm(n = 1000, mean = mu, sigma = sigma)
par(mfrow = c(3,3))
for (i in 1:3) for (j in 1:3)
plot(x[,i], x[,j], xlab = paste0("Var", i), ylab = paste0("Var", j))
Var1 against Var3 forms a tight negatively sloping cloud (\(\rho=-0.9\)); Var1 against Var2 is a weak positive cloud; Var2 against Var3 is a circular blob with no tilt (\(\rho=0\)).
Caution. Var2 and Var3 are uncorrelated, yet each is correlated with Var1. Marginal independence of a pair does not survive conditioning: if you condition on Var1, Var2 and Var3 become correlated. This is exactly why multi-trait evaluation cannot be replaced by a series of single-trait analyses.
Predicting genetic values by conditioning. Using the two-block partition, complete the R script below to obtain predicted genetic values from observed phenotypes. Block 1 is the genetic values, block 2 the phenotypes.
mu1 <- c(......)
mu2 <- c(......)
V11 = V12 = V21 <- matrix(c(.......), ncol = 2)
V22 <- matrix(c(.......), ncol = 2)
x2 <- c(......)
x1 <- mu1 + V12 %*% solve(V22) %*% (x2 - mu2)
print(x1)
Take two half sibs (\(F_{ij}=0.25\)) with \(h^2=0.3\), \(\sigma_y^2=100\), \(\mu=50\), and observed phenotypes \(y_1=62\), \(y_2=45\).
With \(h^2\sigma_y^2=30\) and \(F_{ij}h^2\sigma_y^2=0.25\times 30=7.5\):
mu1 <- c(0, 0) # genetic values have mean 0
mu2 <- c(50, 50) # phenotypes have mean mu
V11 = V12 = V21 <- matrix(c(30, 7.5, 7.5, 30), ncol = 2)
V22 <- matrix(c(100, 7.5, 7.5, 100), ncol = 2)
x2 <- c(62, 45)
x1 <- mu1 + V12 %*% solve(V22) %*% (x2 - mu2)
print(x1)
Note that \(V_{12}\) and \(V_{22}\) differ only in their diagonals: the genetic covariance between the two animals is \(7.5\) in both, but the phenotypic variance on the diagonal of \(V_{22}\) is \(100\) whereas the genetic variance in \(V_{12}\) is \(30\).
Evaluating gives approximately \(\hat g_1\approx 3.5\) and \(\hat g_2\approx -0.6\). Two features are worth noting. First, the predictions are regressed towards zero relative to the phenotypic deviations (\(+12\) and \(-5\)), because \(h^2<1\). Second, each animal's prediction is influenced by its half sib's record through the off-diagonal \(7.5\), which is precisely the information-sharing that makes BLUP more accurate than individual phenotypes alone.
Decomposing a phenotypic correlation. Two traits have \(h_1^2=0.25\) and \(h_2^2=0.49\). The genetic correlation is \(\rho_{g_{12}}=-0.50\) and the environmental correlation is \(\rho_{e_{12}}=+0.40\).
1. Apply the decomposition:
2. The phenotypic correlation is small and positive, about \(+0.07\). A breeder seeing only this would conclude the two traits are essentially unrelated, and that selecting on trait 1 is harmless for trait 2.
That conclusion is incorrect, and costly. The genetic correlation is \(-0.50\), strongly antagonistic. Correlated response depends on \(\rho_{g}\), not \(\rho_{12}\), so selection on trait 1 will steadily degrade trait 2. The positive environmental correlation has masked the antagonism at the phenotypic level. This is the quantitative reason that breeding programmes must estimate genetic parameters rather than relying on observed phenotypic associations.
