Estimating genetic parameters

Practical 3: Pedigree relationships, variance components and heritability

Learning objectives

Construct and interpret additive pedigree relationships; align records by identity; fit a linear mixed model; and distinguish a model-based heritability estimate from a raw phenotype variance ratio.

Continue from Practical 2. Read alongside variance partitioning. You need base R and the qgg package (install.packages("qgg") if needed).

1. Prepare the mouse records and pedigree

Body weight is measured at eight weeks in this F2 teaching population. The phenotype file has 1,177 records; the pedigree includes ancestors. Recorded trait units are not supplied in file metadata. Shared helpers download pinned public data, check their checksums and validate the identifier mapping. Download the helper script with the tutorial, retaining the scripts/ and tutorials/ folders.

Show R code 
source("../scripts/mouse_practical_helpers.R")
data <- mouse_pedigree_data()
mouse <- data$mouse
pedigree <- data$pedigree
ids <- data$ids
c(records = nrow(mouse), pedigree_individuals = nrow(pedigree))
             records pedigree_individuals 
                1177                 1267
Show R code 
head(pedigree)
  id sire  dam family    sex generation
1  1 <NA> <NA>    0/0   Male         M6
2  2 <NA> <NA>    0/0 Female         IC
3  3 <NA> <NA>    0/0   Male         IC
4  4 <NA> <NA>    0/0 Female         M6
5  5 <NA> <NA>    0/0   Male         M6
6  6 <NA> <NA>    0/0 Female         IC
Show R code 
table(pedigree$generation)

  F1   F2   IC   M6 
  66 1177   12   12

Question 1. Why does the pedigree contain more individuals than the phenotype file?

Ancestors establish relationships even when they have no body-weight record in the analysis. The historical pedigree uses 0 for unknown parents; the helper converts it to NA. Source phenotype row labels are checked against explicit pedigree IDs and parent labels, rather than treated as row positions.

2. Calculate pedigree relationships

The additive relationship matrix \(A\) is twice the kinship matrix. Under unrelated, noninbred founder assumptions, founder diagonals are 1. With known parents \(s,d\), the tabular recursion is

\[A_{ij}=\tfrac12(A_{sj}+A_{dj}),\qquad A_{ii}=1+\tfrac12A_{sd}.\]

An unknown parent’s contribution to the off-diagonal recursion is zero. Parents must precede offspring. Here is a small family: two unrelated founders and two full siblings.

Show R code 
family <- data.frame(id = c("sire", "dam", "child1", "child2"),
                     sire = c(NA, NA, "sire", "sire"),
                     dam = c(NA, NA, "dam", "dam"))
A_small <- mouse_relationship(family)
A_small
       sire dam child1 child2
sire    1.0 0.0    0.5    0.5
dam     0.0 1.0    0.5    0.5
child1  0.5 0.5    1.0    0.5
child2  0.5 0.5    0.5    1.0
Show R code 
expected <- matrix(c(1,0,.5,.5, 0,1,.5,.5,
                     .5,.5,1,.5, .5,.5,.5,1), 4, 4)
stopifnot(max(abs(A_small-expected)) < 1e-12)
# Verify the teaching recursion against an installed public interface.
family_qgg <- family
family_qgg[is.na(family_qgg)] <- "0"
stopifnot(max(abs(A_small-qgg::grm(pedigree=family_qgg))) < 1e-12)
A <- mouse_relationship(pedigree)
K <- A[ids, ids, drop = FALSE]
stopifnot(identical(rownames(K), ids), identical(colnames(K), ids))
K[1:5, 1:5]
    91  92  93  94  95
91 1.0 0.5 0.5 0.5 0.5
92 0.5 1.0 0.5 0.5 0.5
93 0.5 0.5 1.0 0.5 0.5
94 0.5 0.5 0.5 1.0 0.5
95 0.5 0.5 0.5 0.5 1.0
Show R code 
range(diag(K))
[1] 1 1

Question 2. Why do the siblings have relationship 0.5, not 1? What does a diagonal above 1 mean?

Full siblings share an expected half of their additive inheritance from unrelated parents; their realized genomes need not be identical. A diagonal is \(1+F_i\), where \(F_i\) is pedigree inbreeding under the declared founder base. Founder relatedness is an assumption, not established by missing parent records.

Show R code 
image(seq_len(40), seq_len(40), K[1:40,1:40],
      xlab="Record index", ylab="Record index", col=hcl.colors(20,"Blues"))
Figure 1: Pedigree relationships among the first 40 phenotype records, in source order.

The full relationship matrix is retained for this small teaching dataset. Plotting a labelled subset is easier to interpret than drawing every matrix entry. Dense whole-population matrices are not the default for large genomic workflows.

3. Specify a variance-component model

For the recorded individuals,

\[y=Xb+u+e,\quad u\sim N(0,A\sigma_A^2),\quad e\sim N(0,I\sigma_E^2),\]

with independent \(u,e\). Fixed effects \(b\) describe sex and replicate means. The pedigree base, additive model and independent residual assumptions determine interpretation. Maternal/shared-litter effects and dominance are not included; family resemblance can therefore be attributed incorrectly if these omissions matter.

Show R code 
X <- model.matrix(~ sex + reps, data=mouse)
y <- mouse$BW
names(y) <- ids
fitA <- mouse_fit(y, X, K, "A")
fitA$theta
        A         E 
 6.569611 13.384147
Show R code 
c(iterations=fitA$niter, qgg_version=as.character(packageVersion("qgg")))
 iterations qgg_version 
        "9"     "1.1.6"
Show R code 
VA <- unname(fitA$theta["A"])
VE <- unname(fitA$theta["E"])
h2 <- VA/(VA+VE)
h2
[1] 0.3292418

We use the installed qgg REML interface. The helper stops if the iteration limit is reached or fitted components are invalid; reaching convergence supports a numerical result, not the biological adequacy of the model.

Question 3. Is \(\sigma_A^2+\sigma_E^2\) the raw variance of body weight?

No. Sex and replicate variation are represented by fixed effects. The model’s genetic and residual scales describe conditional variation; the raw phenotype variance also includes fixed group differences. For individual \(i\), conditional variance is \(A_{ii}\sigma_A^2+\sigma_E^2\). The reported \(h^2\) uses the noninbred base convention, not an individual-specific diagonal adjustment.

4. Interpret heritability and uncertainty

Under that base convention, narrow-sense heritability is

\[h^2=\frac{\sigma_A^2}{\sigma_A^2+\sigma_E^2}.\]

The fitted value is 0.329. It describes the model and population; it is not the fraction of an individual’s phenotype caused by genes, nor a universal property of body weight.

The installed qgg 1.1.6 implementation returns the estimated covariance matrix of the variance components in asd. Apply the delta method to propagate its uncertainty to the ratio. Check the version before relying on this representation.

Show R code 
stopifnot(packageVersion("qgg") == "1.1.6",
          is.matrix(fitA$asd), all(dim(fitA$asd)==c(2,2)))
gradient <- c(VE, -VA)/(VA+VE)^2
se_h2 <- sqrt(drop(t(gradient) %*% fitA$asd %*% gradient))
data.frame(h2=h2, approximate_SE=se_h2,
           Wald_lower=h2-1.96*se_h2, Wald_upper=h2+1.96*se_h2)
         h2 approximate_SE Wald_lower Wald_upper
1 0.3292418     0.09276235  0.1474276   0.511056

This is a local asymptotic approximation, not a bounded confidence interval procedure. Do not clip its endpoints and call that a calibrated interval. Limited founder families, omitted shared environment and relationship assumptions can matter beyond this reported uncertainty.

Regression on relatives: a complementary model

Regression on relatives connects family resemblance to the variance model. In an illustrative noninbred population, suppose \(V_P=V_A+V_E\), mates are unrelated, and environmental deviations are independent across individuals and generations. Additive inheritance gives covariance \(V_A/2\) between a parent’s phenotype and its offspring’s phenotype. The offspring-on-one-parent regression slope is therefore

\[b_{\mathrm{one\ parent}}=\frac{V_A/2}{V_P}=\frac12h^2.\]

For the midparent phenotype \(M=(P_s+P_d)/2\), unrelated mates give \(\operatorname{Var}(M)=V_P/2\) and \(\operatorname{Cov}(P_o,M)=V_A/2\). Thus

\[b_{\mathrm{midparent}}=\frac{V_A/2}{V_P/2}=h^2.\]

At a hypothetical \(h^2=0.4\), these slopes are 0.2 and 0.4. They are expected slopes under the stated population model, rather than regressions fitted to the mouse data. Assortative mating, shared environment, maternal effects or a different sampling population can change their interpretation. Family resemblance alone does not separate all these sources.

Interpretation question. Why would doubling every regression slope on relatives fail to give a general heritability estimator?

The factor depends on which relatives and observation units enter the regression and on their covariance model. Doubling the one-parent slope gives \(h^2\) under the stated assumptions; doubling the midparent slope does not. Environmental covariance can contribute to either slope. The mouse mixed model uses its pedigree and declared fixed and random effects rather than replacing them with a generic family regression.

Exercises

  1. Calculate phenotype variance and compare it with the sum of model components.
  2. Add a hypothetical offspring of child1 and child2 to the small family. Calculate its pedigree inbreeding.
  3. Explain why changing environment or population can change heritability even if gene effects do not change.
Show R code 
c(raw_variance=var(y), model_component_sum=VA+VE)
       raw_variance model_component_sum 
           37.84458            19.95376
Show R code 
inbred_family <- rbind(family,
  data.frame(id="offspring", sire="child1", dam="child2"))
A_inbred <- mouse_relationship(inbred_family)
F_offspring <- A_inbred["offspring","offspring"]-1
stopifnot(abs(F_offspring-0.25)<1e-12)
F_offspring
[1] 0.25

The offspring of the full siblings has \(F=0.25\) under this founder base. Heritability is a ratio: altered allele frequencies, genetic base, environmental variance or measurement design can change it.

Interactive exploration

Open Heritability Estimation Using Different Types of Relatives to explore regression on relatives. Compare the relatives and observation units used in the app with the one-parent and midparent models above, and explain the assumptions behind a heritability interpretation.

Next practical

Practical 4 uses the fitted relationship model to predict breeding values and explains shrinkage, prediction uncertainty and ranking.