Your first evaluation, starting from your own files
Source:vignettes/getting-started.Rmd
getting-started.RmdThis guide goes from two files on your disk to breeding values you can act on. It assumes you know what a breeding value is and does not assume you know anything about this package. Everything here runs: the files are written first so the reading code is real, and you can replace those two lines with your own paths.
What you need
Two text files.
A phenotype file, one row per record. It needs the animal’s identifier, the trait you are evaluating, and the fixed effects you want to correct for. Contemporary group is almost always one of them.
A pedigree file, one row per animal, with three
columns in this order: animal, sire, dam. An unknown parent is
0. Every animal that appears as a parent must have its own
row, and the package will stop rather than guess if one is missing. A
parent quietly turned into an unknown changes the Mendelian variance of
its offspring and the relationship of everything downstream, so it is an
error and not a warning. An animal of the data with no row at all is a
different case: it enters as a founder, with unknown parents, and the
fit warns with how many and which.
The identifiers must be the same strings in both files. This is the single most common reason a first fit returns nothing useful.
Let me write two files so the rest of the guide can read them. Yours already exist, so skip this chunk and point the next one at your paths.
Read them in
dados <- read.csv(file.path(pasta, "phenotypes.csv"), colClasses = "character")
ped <- read.csv(file.path(pasta, "pedigree.csv"), colClasses = "character")
dados$y <- as.numeric(dados$y) # the trait is the only numeric column
str(dados)
#> 'data.frame': 420 obs. of 3 variables:
#> $ id: chr "F001" "F002" "F003" "F004" ...
#> $ cg: chr "g3" "g3" "g2" "g1" ...
#> $ y : num 9.29 10 9.57 10.46 9.79 ...
head(ped, 3)
#> id sire dam
#> 1 F001 0 0
#> 2 F002 0 0
#> 3 F003 0 0Reading everything as character first and converting only the trait
is deliberate. Identifiers like 0012345 lose their leading
zero when R decides they are numbers, and then they stop matching the
pedigree. Convert the trait, and any covariate, by hand.
Look at the data before you model it
describe() reports the shape of what you have, and
suggest_model() reads that shape and names the terms it
calls for. Neither fits anything.
describe(dados, trait = "y", classes = "cg", pedigree = ped)
#> Trait 'y': 420 record(s)
#> mean 9.9266, sd 1.0651
#> quartiles 6.633 / 9.170 / 9.962 / 10.686 / 13.027
#> 3 record(s) outside the Tukey fences [6.896, 12.96]
#> class cg: 3 level(s), smallest 139, median 140, 0 with a single record
#> pedigree: 420 animals, 60 founder(s), mean F 0.005655, maximum 0.25
suggest_model(dados, trait = "y", animal = "id")
#> the shape of the data asks for nothing beyond the base modelqc_phenotypes() flags problems. It marks the missing
code, turns outliers beyond a Tukey fence into missing values, and
reports class levels too small to estimate. It flags rather than
deletes: dropping a row would reshape contemporary groups
without saying so, and you would never see it happen.
q <- qc_phenotypes(dados, "y", classes = "cg")
q[c("n_missing", "n_outliers", "small_classes")]
#> $n_missing
#> [1] 0
#>
#> $n_outliers
#> [1] 0
#>
#> $small_classes
#> $small_classes$cg
#> [1] 0
dados <- q$dataThe fit
The model is a formula. An unmarked term is a fixed class effect,
cov(x) is a covariate, and animal(id) is the
additive genetic effect that the pedigree gives structure to.
fit <- model(y ~ cg + animal(id), dados, ped)
fit
#> AI-REML fit of 'y'
#> converged in 7 iteration(s), relDelta 4.88e-09, Newton decrement 6.56e-18, 0.006 s
#> -2logL 432.5132793
#> 420 record(s), 423 column(s) in the equations
#> fixed column(s) removed for linear dependence: cg=g1
#>
#> component estimate std_error share correlation
#> var(animal) 0.483684 0.1135616 0.4322
#> var(residual) 0.635468 0.0806489 0.5678
#>
#> fixed effects (implicit intercept; dropped columns are zero; compare by contrast):
#> intercept cg=g3 cg=g2
#> 9.908111 0.164347 -0.171770Read that from the top. The variance components come from AI-REML.
converged says whether it got there, and it is not a
formality: the package refuses to call a fit converged unless the Newton
decrement of the free components passes tolerance as well, so a
FALSE there is a real statement and not a rounding
complaint.
Reading the answer
Heritability is the additive share of the phenotypic
variance. h2() picks the numerator by name and puts the
covariances in the denominator, where a direct-maternal covariance
belongs (Willham, 1972); it is the same number as the share
of var(animal) in the printed fit:
The file was simulated with a heritability of 0.35, so this is the estimate landing near a truth you happen to know. On your own data you will not have that check, which is why the last section of this guide is about what to look at instead.
Breeding values and their accuracy:
e <- ebv(fit)
acc <- accuracy(fit, ped)
head(round(e, 3), 5)
#> F001 F002 F003 F004 F005
#> -0.192 -0.058 0.282 0.456 0.196
summary(acc)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.7051 0.7244 0.7437 0.7480 0.7654 0.8355ebv() returns every animal of the data, every genotyped
animal and all their ancestors in the pedigree, including parents
without records, because a parent gets an estimate from its offspring.
The fit prunes the animals with no record, no genotype and no recorded
descendant, and a message says how many lines it left out; to rank such
a young candidate, give it a row of the data with the trait missing.
Accuracy carries the (1 + F) of the pedigree in its
denominator, which is the honest one: an inbred animal has less
Mendelian variance left to predict.
Ranking, which is usually the actual question.
solutions() puts the breeding value, its standard error and
its accuracy in one table, sorted:
head(solutions(fit, ped), 5)
#> id ebv se acc
#> 141 G1_081 1.412237 0.4288105 0.7872976
#> 137 G1_077 1.303333 0.4641909 0.7446586
#> 329 G3_029 1.267224 0.4840481 0.7180441
#> 11 F011 1.235056 0.3939246 0.8241224
#> 7 F007 1.229029 0.4026998 0.8153067Read the accuracy column before you act on the ranking. A high breeding value with low accuracy is a young animal with little information, not a proven one.
The complication you will actually hit: repeated records
If an animal has more than one record, the additive term alone will
absorb what those records share and heritability will come out too high.
The permanent environment term is what catches it. model()
notices the gap: the first fit below warns that every animal has more
than one record and the model has no pe(id), and says what
to add.
rep_dados <- rbind(dados, transform(dados, y = dados$y + rnorm(nrow(dados), 0, 0.5)))
sem_pe <- model(y ~ cg + animal(id), rep_dados, ped)
#> Warning: animal(id): 420 of the 420 animals with records of y have more than
#> one record (840 records between them) and the model has no pe(id): the
#> permanent-environment covariance between records of the same animal goes into
#> var(animal), which inflates h2 by more the larger the share of repeated
#> records; add pe(id), the repeatability model
com_pe <- model(y ~ cg + animal(id) + pe(id), rep_dados, ped)
c(sem_pe = round(h2(sem_pe), 3), com_pe = round(h2(com_pe), 3))
#> sem_pe com_pe
#> 0.940 0.369The first number is inflated and the second is the heritability. The difference between them is exactly the within-animal covariance that had nowhere else to go. Repeatability is the two shares added together.
When it does not converge
converged = FALSE has a small number of causes, and the
message says which.
Hitting the iteration cap means it was still moving. Raise
maxiter=. Starting values far from the answer make this
common, and start= from a simpler fit usually fixes it. The
exception is a cap reached with the Newton decrement already under
tolerance: the fit then reports converged = TRUE and says
the likelihood is flat along some direction, because more iterations
would not choose a point on that ridge.
A message that says SINGULAR names components the data cannot
separate. Their standard errors are NaN and their values depend on
start=; only combinations of them are estimable. A
pe(id) on data with one record per animal is the common
case: there it is the residual under another name.
A component walking to zero means the data does not support that term. This is information, not a failure: a permanent environment that goes to zero is telling you the repeated records share nothing beyond the genetics.
A correlation stuck at plus or minus one, in a model with a
covariance group such as direct plus maternal, is a real boundary and
not a numerical accident. The fit says so in its message. A standard
error there is meaningless, and the interval has to come from
profile_theta() instead.
fit$converged
#> [1] TRUE
fit$message
#> [1] ""Adding genotypes
Genotypes are an argument, not a different program. The matrix is animals by markers, coded 0, 1, 2, with a vector of identifiers that match the pedigree.
g <- qc_genotypes(s$genotypes, min_maf = 0.01)
g[c("n_animal", "n_call", "n_maf", "passes")]
#> $n_animal
#> [1] 0
#>
#> $n_call
#> [1] 0
#>
#> $n_maf
#> [1] 0
#>
#> $passes
#> [1] 1
fit_g <- model(y ~ cg + animal(id), dados, ped, genotypes = g)
round(c(pedigree = h2(fit), single_step = h2(fit_g)), 3)
#> pedigree single_step
#> 0.432 0.294
fit_g$dense_block
#> dense columns
#> 421 423That is single-step: genotyped and non-genotyped animals in one
evaluation, with the relationship matrix built from both sources. In
this fit accuracy() divides a genotyped animal by its own
diagonal of G* rather than 1 + F.
fit_g$dense_block says how large the dense corner of the
factor is, the k of the k^3 each factorization pays: without
apy_core= it is at least the number of genotyped
animals.
The two heritabilities are not the same number and there is no reason to expect them to be. The denominator changed. Pedigree relationships are expectations, while genomic ones are what actually happened at the markers, so full sibs stop being 0.5 apiece and spread around it. Which of the two you report depends on which base population you mean, and the honest habit is to say which one you fitted rather than to pick the number you like.
If your genotyped set is large enough that the dense route hurts,
apy_core= and vecchia_k= are the two
approximations, and snp_blup() solves the marker equations
without forming the genomic relationship matrix at all.
apy_core_select() chooses the APY core by the eigenvalues
of G that explain 98% of its trace (Pocrnic et al., 2016); choose it
once and pass it to every fit.
Before you believe any of it
Four checks, in the order they catch things.
Do the identifiers match?
length(intersect(dados$id, ped$id)) should be the number of
animals with records. A near-zero intersection is the leading-zero
problem from the first section.
Did quality control remove something that mattered?
n_missing, n_outliers and
small_classes say what it flagged.
Is the contemporary group doing what you think? A group with very many levels can absorb the very signal you are trying to estimate, and a group with one animal in it contributes nothing but a lost degree of freedom.
Does the heritability sit where the literature puts it? If it does not, the reason is usually a missing term in the model rather than a discovery.