library(simplePHENOTYPES)
library(dplyr)
library(tibble)14 Individual, family, and multivariate selection
Suppose you have a segregating population organized in families, with a phenotype for each individual, and you want to choose the parents of the next generation. This chapter selects on individual records, on family information, and on several traits at once, and checks the classical prediction of response against what the simulation actually delivers. The simulation knows every individual’s genetic value, so the genetic consequence of each selection decision can be measured directly.
14.1 Worked example
14.1.1 Load the package and build families
The example starts from the bundled maize inbred lines, restricted to three chromosomes to keep the meiosis simulation fast. Ten lines are crossed to ten others in a factorial design, giving 100 F1 hybrids, and each F1 is selfed to produce ten S1 progeny. The result is 1,000 individuals in 100 families, each family being the selfed progeny of one F1. The crossing chapter explains the meiosis model; here only the resulting population matters.
data("SNP55K_maize282_maf04", package = "simplePHENOTYPES")
geno <- as_numeric(SNP55K_maize282_maf04, to_r = TRUE, to_file = FALSE,
verbose = FALSE)
geno <- geno[geno$chr %in% c("1", "2", "3"), ]
parents <- as_population(geno, individuals = 1:20)
plan_f1 <- mating_design(parents[1:10], parents[11:20], design = "factorial")
f1 <- mate(plan_f1, parents, seed = 1, prefix = "F1")
plan_s1 <- data.frame(mother = parentage(f1)$id, father = parentage(f1)$id, n = 10)
s1 <- mate(plan_s1, f1, seed = 2, prefix = "S1")
family <- families(s1, by = "selfed")
s1#> <Population>
#> Individuals: 1000 Markers: 4134 Chromosomes: 3
#> Genetic map: 560 cM total span (168-219 cM per chromosome)
#> Origin: mate(100 rows)
#> Pedigree: 1120 recorded individuals; generation 2
#> IDs: S1_1, S1_2, S1_3, S1_4, S1_5, S1_6, ... (994 more)table(family) |> head()#> family
#> F1_1 F1_2 F1_3 F1_4 F1_5 F1_6
#> 10 10 10 10 10 10mating_design() writes the mother-father plan and mate() executes it. A plan row with the same individual as mother and father is a self. families() labels each S1 individual by the F1 it descends from, which is the grouping the family methods need.
14.1.2 Simulate the phenotype
A single additive trait with 50 causal loci and a genetic share of 0.3 is simulated on the S1 population. The additive chapter explains the model; the selection functions only require a realized simulation object.
sim <- simulate_phenotype(s1, h2 = 0.3, n_qtn = 50, seed = 3)
sim#> <phenotype_sim> (realized · long format)
#> Genotypes: <Population: mate(100 rows)> Traits: 1 Architecture: independent Seed: 3
#> Variance partition (proportions of V_P):
#> additive 0.30 (50 QTNs, geometric)
#> residual 0.70
#> Requested genetic share = 0.30 realized H² = 0.31pheno <- phenotypes_wide(sim)
gv <- genetic_values(sim)[, "Trait_1"]phenotypes_wide() gives one row per individual, and genetic_values() returns the simulated genetic values that selection cannot see but the book can use to judge each method.
14.1.3 Mass selection and the selection intensity
select_ind() ranks individuals on a criterion and keeps the best. The criterion on = "pheno" is the observed phenotype, which is what a breeder without markers can select on. Exactly one of n, prop, or intensity sets how many are kept.
top <- select_ind(sim, prop = 0.1, on = "pheno", method = "mass")
top#> <Population>
#> Individuals: 100 Markers: 4134 Chromosomes: 3
#> Genetic map: 560 cM total span (168-219 cM per chromosome)
#> Origin: mate(100 rows)
#> Pedigree: 175 recorded individuals; generation 2
#> IDs: S1_605, S1_769, S1_803, S1_667, S1_754, S1_829, ... (94 more)tibble(
attribute = c("Individuals kept", "Selection differential",
"Realized intensity", "Method"),
value = c(n_individuals(top), attr(top, "differential"),
attr(top, "intensity"), attr(top, "method"))
)| attribute | value |
|---|---|
| Individuals kept | 100 |
| Selection differential | 1.75086837880987 |
| Realized intensity | 1.76766950479857 |
| Method | mass |
The result is a Population of the selected individuals, ready to be crossed, carrying the selection differential \(S\) (mean criterion of the selected minus mean of all candidates) and the realized intensity \(i = S / s_P\), where \(s_P\) is the sample standard deviation of the criterion. The table compares the realized intensity for several proportions with the large-population value \(i(p) = \phi(z_p)/p\) used when intensity is given instead of prop.
proportions <- c(0.5, 0.2, 0.1, 0.05, 0.02)
intensity <- tibble(
prop = proportions,
kept = sapply(proportions, function(p) n_individuals(select_ind(sim, prop = p))),
realized_i = sapply(proportions, function(p) attr(select_ind(sim, prop = p), "intensity")),
expected_i = dnorm(qnorm(1 - proportions)) / proportions
)
intensity| prop | kept | realized_i | expected_i |
|---|---|---|---|
| 0.50 | 500 | 0.8043168 | 0.7978846 |
| 0.20 | 200 | 1.3984455 | 1.3998096 |
| 0.10 | 100 | 1.7676695 | 1.7549833 |
| 0.05 | 50 | 2.0984916 | 2.0627128 |
| 0.02 | 20 | 2.3520847 | 2.4209068 |
by_intensity <- select_ind(sim, intensity = 1.75)
n_individuals(by_intensity)#> [1] 101The realized intensities lie close to the large-population values. With intensity = 1.75, the function keeps the count whose large-population intensity is nearest to the request.
14.1.4 The breeder’s equation
For mass selection on a purely additive trait, the classical prediction is that the genetic gain equals the heritability times the selection differential, \(R = h^2 S\). The simulation lets us compare that prediction with the genetic selection differential actually achieved: the mean genetic value of the selected individuals minus the mean of all candidates.
genetic_gain <- function(selected) {
mean(gv[attr(selected, "selected")]) - mean(gv)
}
response <- tibble(
prop = proportions,
S = sapply(proportions, function(p) attr(select_ind(sim, prop = p), "differential")),
predicted = 0.3 * S,
simulated = sapply(proportions, function(p) genetic_gain(select_ind(sim, prop = p)))
)
response| prop | S | predicted | simulated |
|---|---|---|---|
| 0.50 | 0.7966721 | 0.2390016 | 0.2455902 |
| 0.20 | 1.3851537 | 0.4155461 | 0.4262686 |
| 0.10 | 1.7508684 | 0.5252605 | 0.5591023 |
| 0.05 | 2.0785461 | 0.6235638 | 0.6624754 |
| 0.02 | 2.3297289 | 0.6989187 | 0.6972038 |
The simulated gains track the prediction, with scatter from the particular residuals drawn. The prediction assumes that the regression of genetic value on phenotype is linear with slope \(h^2\); the theory section below says when that holds.
Selection on the true genetic value (on = "gv") is the idealized upper bound. For this additive model the breeding value (on = "bv") is the same quantity. A direction = "low" selection keeps the smallest values instead.
tibble(
criterion = c("pheno", "gv", "bv", "pheno, direction = low"),
genetic_gain = c(
genetic_gain(select_ind(sim, prop = 0.1, on = "pheno")),
genetic_gain(select_ind(sim, prop = 0.1, on = "gv")),
genetic_gain(select_ind(sim, prop = 0.1, on = "bv")),
genetic_gain(select_ind(sim, prop = 0.1, on = "pheno", direction = "low"))
)
)| criterion | genetic_gain |
|---|---|
| pheno | 0.5591023 |
| gv | 0.8525770 |
| bv | 0.8525770 |
| pheno, direction = low | -0.4334035 |
Any externally computed score, such as the genomic predictions of the prediction chapter, can be passed to on as a named numeric vector or as a function of the simulation object.
14.1.5 The random control
method = "random" draws individuals at random. Its reported differential and intensity are the realized values on the criterion for the individuals drawn, so they scatter around zero. It is the drift control against which a selection method’s gain is judged.
controls <- sapply(1:5, function(k) {
set.seed(k)
genetic_gain(select_ind(sim, prop = 0.1, method = "random"))
})
tibble(replicate = 1:5, genetic_gain = controls)| replicate | genetic_gain |
|---|---|
| 1 | 0.0064095 |
| 2 | -0.0375071 |
| 3 | 0.0100098 |
| 4 | -0.0178430 |
| 5 | 0.0474389 |
14.1.6 Response in the next generation
A selected Population remembers the trait it was selected on, so progeny can be phenotyped on the same scale. The selected individuals are intermated at random, as are the individuals of a random control, and the progeny are scored with the carried trait. population_trait() shows what is carried.
population_trait(top)| trait | h2 | var_e |
|---|---|---|
| Trait_1 | 0.3 | 0.7 |
control <- select_ind(sim, prop = 0.1, method = "random")
progeny_of <- function(selected, seed) {
plan <- mating_design(selected, design = "random", n_crosses = 50,
progeny_per_cross = 10, seed = seed)
mate(plan, selected, seed = seed + 1)
}
next_sel <- simulate_phenotype(progeny_of(top, 10), seed = 12)
next_ctrl <- simulate_phenotype(progeny_of(control, 20), seed = 22)
tibble(
population = c("All candidates", "Selected parents", "Random-control parents",
"Progeny of selected", "Progeny of random control"),
mean_genetic_value = c(mean(gv), mean(gv[attr(top, "selected")]),
mean(gv[attr(control, "selected")]),
mean(genetic_values(next_sel)[, "Trait_1"]),
mean(genetic_values(next_ctrl)[, "Trait_1"])),
mean_phenotype = c(mean(pheno$Trait_1), NA, NA,
mean(phenotypes_wide(next_sel)$Trait_1),
mean(phenotypes_wide(next_ctrl)$Trait_1))
)| population | mean_genetic_value | mean_phenotype |
|---|---|---|
| All candidates | 0.0000000 | 0.0000000 |
| Selected parents | 0.5591023 | NA |
| Random-control parents | 0.0273210 | NA |
| Progeny of selected | 0.5916472 | 0.5916472 |
| Progeny of random control | 0.1371142 | 0.1371142 |
The progeny of the selected parents have a higher mean than the progeny of the random control, and each progeny mean is close to its parents’ mean genetic value, as expected for an additive trait; the remaining gap reflects which parents the random plan happened to use and how often. Because the carried trait keeps the base population’s loci, effects, and residual variance, these means are comparable with the candidates’ means. The progeny phenotype means equal the genetic means because the default fixed residual mode centres the residuals; see the first chapter. The trait chapter explains the carried trait in detail.
14.1.7 Family selection
When individuals are grouped in families, their family mean carries information about their breeding values. Three methods use it. within_family keeps the best individuals inside each family, among_family keeps whole families with the best means, and combined ranks on an index that weights the individual’s own deviation and its family’s deviation. The combined index needs the trait’s heritability and the correlation of breeding values among family members, which for S1 sibs of a non-inbred F1 is two thirds.
within <- select_ind(sim, n = 100, method = "within_family", family = family)
among <- select_ind(sim, n = 100, method = "among_family", family = family)
combined <- select_ind(sim, n = 100, method = "combined", family = family,
h2 = 0.3, family_relationship = 2 / 3)
mass <- select_ind(sim, n = 100, method = "mass")
tibble(
method = c("mass", "within_family", "among_family", "combined"),
kept = c(n_individuals(mass), n_individuals(within),
n_individuals(among), n_individuals(combined)),
families_represented = c(
length(unique(family[attr(mass, "selected")])),
length(unique(family[attr(within, "selected")])),
length(unique(family[attr(among, "selected")])),
length(unique(family[attr(combined, "selected")]))
),
genetic_gain = c(genetic_gain(mass), genetic_gain(within),
genetic_gain(among), genetic_gain(combined))
)| method | kept | families_represented | genetic_gain |
|---|---|---|---|
| mass | 100 | 55 | 0.5591023 |
| within_family | 100 | 100 | 0.1655628 |
| among_family | 100 | 10 | 0.6464029 |
| combined | 100 | 16 | 0.6983419 |
Within-family selection keeps every family represented, which conserves diversity at the cost of gain. Among-family selection keeps whole families, so it returns at least n individuals and concentrates the selection in few families. The combined index uses both sources of information. Which method gains most depends on the heritability and the family structure; at this moderate heritability the methods that use the family mean are competitive with mass selection.
n_per_family keeps a stated number from every family instead of allocating a total.
one_each <- select_ind(sim, method = "within_family", family = family,
n_per_family = 1)
n_individuals(one_each)#> [1] 100head(attr(one_each, "n_per_family"))#> F1_1 F1_10 F1_100 F1_11 F1_12 F1_13
#> 1 1 1 1 1 114.1.8 Selection on several traits
Multi-trait selection needs a multi-trait simulation. Two pleiotropic traits with a requested genetic correlation of 0.3 are simulated on the same population; the multi-trait chapter explains the architecture. The index weights use the simulation’s true breeding-value covariance; the Smith-Hazel index then ranks phenotypes, while the quadratic index ranks the true breeding values.
sim2 <- simulate_phenotype(s1, architecture = "pleiotropy", n_traits = 2,
h2 = c(0.4, 0.4), n_qtn = 50,
cor = matrix(c(1, 0.3, 0.3, 1), 2), seed = 4)
gv2 <- genetic_values(sim2)
gain2 <- function(selected) {
colMeans(gv2[attr(selected, "selected"), , drop = FALSE]) - colMeans(gv2)
}
index <- select_ind(sim2, n = 100, method = "index", weights = c(1, 1))
index_t1 <- select_ind(sim2, n = 100, method = "index", weights = c(1, 0))
linear_bv <- select_ind(sim2, n = 100, method = "quadratic_index",
weights = c(1, 1))
quadratic <- select_ind(sim2, n = 100, method = "quadratic_index",
weights = c(1, 1),
quad_weights = matrix(c(-1, 0, 0, 0), 2))
spread_t1 <- function(selected) sd(gv2[attr(selected, "selected"), "Trait_1"])
bind_rows(
`index, equal weights` = gain2(index),
`index, trait 1 only` = gain2(index_t1),
`quadratic index, linear part only` = gain2(linear_bv),
`quadratic index, penalty on trait 1` = gain2(quadratic),
.id = "method"
) |>
mutate(sd_trait1_selected = c(spread_t1(index), spread_t1(index_t1),
spread_t1(linear_bv), spread_t1(quadratic)))| method | Trait_1 | Trait_2 | sd_trait1_selected |
|---|---|---|---|
| index, equal weights | 0.6577878 | 0.6590633 | 0.5530474 |
| index, trait 1 only | 0.7458763 | 0.5009216 | 0.4943117 |
| quadratic index, linear part only | 0.9881478 | 0.8609306 | 0.3918977 |
| quadratic index, penalty on trait 1 | 0.4863629 | 0.9543652 | 0.3016330 |
method = "index" is the Smith-Hazel index: economic weights, one per trait, are converted to index weights through the phenotypic and genetic covariance matrices, and individuals are ranked on the weighted phenotypes. Putting all weight on the first trait still moves the second because the traits are genetically correlated. method = "quadratic_index" scores the true breeding values directly, as the package documentation states, so its gains are an idealized reference and are not comparable with the phenotype-based index. Its quadratic term, here a penalty on the square of the first trait, pulls the selected set toward intermediate values of that trait: the gain on trait 1 drops and so does the spread of the selected individuals, which is the form for traits with an intermediate optimum.
Independent culling levels keep an individual only if it passes a threshold on every trait. sequential = TRUE culls the traits in order, each among the survivors of the previous ones.
culled <- select_ind(sim2, method = "culling", culling = c(0.4, 0.4))
culled_seq <- select_ind(sim2, method = "culling", culling = c(0.4, 0.4),
sequential = TRUE)
bind_rows(
`simultaneous culling` = c(kept = n_individuals(culled), gain2(culled)),
`sequential culling` = c(kept = n_individuals(culled_seq), gain2(culled_seq)),
.id = "method"
)| method | kept | Trait_1 | Trait_2 |
|---|---|---|---|
| simultaneous culling | 187 | 0.5398436 | 0.5102863 |
| sequential culling | 160 | 0.5703585 | 0.5398739 |
With culling the number kept follows from the proportions, so n is not given. Simultaneous culling at 0.4 on each of two traits keeps fewer than 40 percent, because an individual must pass both thresholds. Sequential culling at the same proportions keeps exactly 0.4 of 0.4.
14.1.10 The method manifest
selection_methods() lists every selection and scoring operator the package exports, with the function that implements it. Front ends build their interfaces from this table.
selection_methods() |>
select(id, fn)| id | fn |
|---|---|
| mass | select_ind(method = “mass”) |
| within_family | select_ind(method = “within_family”) |
| among_family | select_ind(method = “among_family”) |
| combined | select_ind(method = “combined”) |
| index | select_ind(method = “index”) |
| quadratic_index | select_ind(method = “quadratic_index”) |
| bqp | select_ind(method = “bqp”) |
| culling | select_ind(method = “culling”) |
| tandem | pedigree(trait = c(…)); recurrent_selection(trait = c(…)) |
| random | select_ind(method = “random”) |
| ocs | optimum_contribution() |
| usefulness | cross_usefulness() |
| mabc | mabc_select() |
| mas | marker_select() |
| combining_ability | combining_ability() |
| progeny_test | progeny_test() |
| blup | predict_ebv() |
| mt_blup | predict_ebv() |
14.2 Assumptions and theory
The worked example uses an additive trait, a single record per individual, families of one type, and candidates from one generation. There is no dominance, no shared family environment, and no genotype-by-environment interaction.
14.2.1 Differential, intensity, and the count kept
With prop, the number kept is \(k = \max(1, \operatorname{round}(pN))\) for \(N\) candidates, a package definition in .resolve_keep(). The selection differential and realized intensity reported as attributes are
\[ S = \bar y_{\text{sel}} - \bar y, \qquad i = \frac{S}{s_y}, \tag{14.1}\]
with \(s_y\) the sample standard deviation of the criterion over all candidates. With intensity, the count is the one whose large-population intensity \(i(p) = \phi(z_p)/p\) is nearest to the request, where \(z_p\) is the standard normal quantile cutting off the upper fraction \(p\) and \(\phi\) is the standard normal density. The equation map attributes these definitions to Falconer & Mackay 1996 (page unverified). The large-population formula overstates the expected realized intensity in small samples.
14.2.2 Response to mass selection
The package documentation states the classical response \(R = i\,h^2\sigma_P =
h^2 S\) for a purely additive trait, and the more general \(R = i\,\mathrm{Cov}(A, P)/\sigma_P\), the regression of breeding value \(A\) on the criterion \(P\). The equation map attributes the response to Falconer & Mackay 1996 (page unverified) and lists the general index form under Smith 1936 Ann. Eugen. 7:240–250 and Hazel 1943 Genetics 28:476–490 (page unverified). Two conditions are needed for the simple form. The regression of \(A\) on \(P\) must be linear, which holds under joint normality but not for an arbitrary joint distribution. And \(\mathrm{Cov}(A, P)\) must equal the additive variance, which requires every non-additive component of the phenotype to be uncorrelated with the breeding value. With an additive model and independent residuals, as in this chapter, both hold up to sampling. With dominance or epistasis under linkage disequilibrium, or in inbred and selected populations, the covariance need not equal the additive variance; the package documentation for select_ind() gives counterexamples.
The chapter’s check compared the genetic selection differential of the selected parents with \(h^2 S\). The progeny mean then follows from the parents’ breeding values under random mating of the selected set, which the progeny table illustrated.
14.2.3 The combined index
method = "combined" scores each individual on the Lush index of its own deviation and its family’s deviation from the overall mean. With \(h^2\) the heritability, \(r\) the correlation of breeding values among family members, \(n\) the family size, and \(t = r h^2\) the phenotypic intraclass correlation implied by an additive model, the weights are
\[ b_1 = \frac{h^2(1 - r)}{1 - r h^2}, \qquad b_2 = \frac{h^2\, n\, r\,(1 - h^2)}{[1 + (n - 1) r h^2](1 - r h^2)}, \qquad I = b_1 (y - \bar y) + b_2 (\bar y_f - \bar y), \tag{14.2}\]
the solution \(b = V^{-1}c\) of the selection-index equations. The equation map records the reference Lush 1947 81:241–261, 362–379 (verified). The weights are optimal under an additive model in which family members covary only through breeding values; dominance, epistasis, or a shared environment add to the covariance of full-sib or selfed families, so the index is then reasonable but not optimal. The value \(r = 2/3\) for S1 sibs of a non-inbred parent is the package’s derivation from the tabular relationship matrix; the equation map notes that this attribution has no primary page.
14.2.4 The Smith-Hazel index and culling
method = "index" computes \(b = P^{-1} G a\) and ranks on \(I = b'y\), with \(P\) the sample phenotypic covariance matrix, \(G\) the sample covariance matrix of the true breeding values, and \(a\) the economic weights. The equation map cites Smith 1936 7:240–250 and Hazel 1943 28:476–490. When \(P\) is singular a minimum-norm solution is used, a package extension without a genetic source. Hazel’s \(G\) is the covariance of phenotypes with breeding values; the two coincide only when breeding values and dominance deviations are uncorrelated, which holds under random mating but not in general in selfed or selected populations.
method = "culling" is independent culling levels, attributed in the equation map to Hazel & Lush 1942 33(11):393–399 (verified). Under their idealized conditions of equal, uncorrelated traits, the package derivation gives aggregate responses proportional to \(i(p)\sqrt{T}\) for the index, \(T\, i(p^{1/T})\) for simultaneous culling, and \(i(p)\) for tandem selection, so index selection is at least as efficient as culling, which is at least as efficient as tandem selection.
14.3 Implementation and practical limits
select_ind() is an R function. The internal helpers named in the developer index include .resolve_keep() for the count, .combined_score() for the Lush index, .index_weights() for the Smith-Hazel weights, .quadratic_index_score(), .select_culling(), and the breeding-value construction .breeding_value_matrix() that the "bv" criterion and the index methods use.
- The Smith-Hazel index builds its weights from the simulation’s true breeding-value covariance and ranks phenotypes; the quadratic index ranks the true breeding values themselves. Neither uses estimates. To select on predictions, compute the index from the predicted values and pass it through
onwithmethod = "mass". "bv"is refused for models with an epistasis layer, because an epistatic term has no per-locus average effect to reconstruct.among_familyreturns at leastnindividuals, and the reported intensity describes the whole returned set.- The combined index uses one heritability and one family relationship for all families, so it assumes families of one type from comparable parents.
- The random control’s differential is a realized quantity and can take either sign; replicate it over seeds before reading anything into one draw.
- Breeding scheme wrappers such as
recurrent_selection()apply mass selection in each cycle; see the breeding schemes chapter.
TipCarrying the selection forward
Every result of select_ind() on a Population-backed simulation is itself a Population with the trait attached. Crossing it, as the progeny example did, is the first step of a breeding cycle; the contributions chapter shows how to choose parental contributions rather than a fixed number of parents.