Gibbs sampler for a single-trait mixed model

The structure checks of model() (see its section "Structure checks") apply here: a design that one of them finds unidentified stops the chain unless theta_fixed = gives the components (a random term inside a fixed effect leaves the chain walking the prior of that variance), and the warnings go to message. Those checks read the structure of the design; a ridge that only the likelihood shows is not among them, and the chain has no SINGULAR verdict where model() reports the information matrix SINGULAR. Measured on 30 data sets of 1600 records each (validation/indirect_gibbs_recovery.R): with pens all of 4 animals, two full sibs from each of two families (indirect()), model() was SINGULAR in 30 of 30 and the message of the chain empty in 29 (in the 30th a chain got the boundary warning described in prior, not a SINGULAR verdict); with records on unrelated founders only (var(animal) against var(residual)), SINGULAR in 30 of 30 and the message empty in 30 of 30. On such a ridge the chain wanders (R-hat up to 2.5 over the separate components, while the combinations the data identify are recovered), and under the default prior it can end on that boundary. Fit model() with the same formula first and read its SINGULAR verdict. The h2 row of summary() and of mcmc_diagnostics() follows the rule of h2(): the direct genetic term, 4 var(sire) in a sire model, the components as estimated (the base of Gamma with metafounders), a kernel(fixed =) or an associative_matrix() at its mean diagonal over the records and the grandsire part of sire(mgs =) in the denominator; with only pe() or random() terms there is no h2 row.

Usage

gibbs(
  formula,
  data,
  pedigree = NULL,
  genotypes = NULL,
  blend = 0.05,
  apy_core = NULL,
  missing_code = NULL,
  vecchia_k = NULL,
  n_iter = 20000L,
  burnin = 2000L,
  thin = 10L,
  theta_fixed = NULL,
  metafounders = NULL,
  gamma = NULL,
  prior = "jeffreys",
  family = c("gaussian", "probit"),
  chains = 1L,
  cores = 1L,
  verbose = interactive(),
  residual_draws = 20L
)

Arguments

formula

as in model() (all the same markers, including rn() and indirect())

data

data.frame

pedigree

data.frame animal, sire, dam. The ids are matched to the data by the rule of model() (see pedigree there): an absent animal enters as a founder, an unknown 0 is no effect, the same id written two ways stops the fit, and the pedigree is pruned to the data, the genotyped and their ancestors

genotypes

single step as in model()

blend

as in model()

apy_core

as in model()

missing_code

missing-value code for the trait

vecchia_k

neighbors per animal in the Vecchia inverse of G*: each genotyped animal conditions on its k strongest previous relationships instead of a global core (the generalization of APY; Henderson's A^-1 is the pedigree case, with the parents as the conditioning set). k >= n - 1 reproduces the exact inverse; the result message says it is an approximation and with which k. Mutually exclusive with apy_core

n_iter

total iterations of the chain

burnin

iterations discarded from the start

thin

keep one sample every thin iterations

theta_fixed

fix the variance components at this vector and sample only the locations (the mode the gates use: the sample mean must reproduce the BLUP and the sample variance the PEV). A kernel() term with fixed = holds that one component and samples the others

metafounders

labels of unknown-parent groups; a parent with one of these labels needs no line of its own (any OTHER cited-without-line parent is still a declared error). Metafounders enter as virtual base rows of A(Gamma) after Legarra et al. (2015)

gamma

the base relationship matrix of the metafounders. Either a vector of length n_metafounders, read as the DIAGONAL, or a full symmetric n_metafounders x n_metafounders matrix. The off-diagonal gamma_jk is the ancestral relationship BETWEEN two base populations, which is what a multibreed analysis turns on; it may be NEGATIVE, for bases pulled apart by selection in opposite directions. Admissibility is positive SEMI-definiteness, tested by a spectral decomposition and not by a Cholesky, plus a diagonal below 2 so that a metafounder's offspring keeps a positive Mendelian variance. A SINGULAR gamma is accepted, through the Moore-Penrose pseudo-inverse: that covers gamma = 0, the unknown-parent-group limit, where the pseudo-inverse reproduces the A-inverse of unknown parent groups exactly, and two metafounders standing for one population, whose rows are identical. What is still refused is an INDEFINITE gamma, a negative eigenvalue, which does not generate a covariance matrix at all

prior

prior of the variance components, the same for every covariance group and the residual. "jeffreys" (the default) is p(s2) ~ 1/s2, and |C|^-(d+1)/2 for a d x d group: the full conditional is (u'K^-1 u) / chisq(nl). It is improper, and so, in principle, is the posterior of every random effect whose likelihood stays positive at a variance of zero (Hobert and Casella 1996). For a scalar variance the divergence at zero is logarithmic: with a well identified component the spurious mass near zero is negligible, with a weak one it is not (the residual too: on 1600 records and 2000 animals with indirect(), 12.6% of the kept draws of one chain had var(residual) below 0.1, against a REML estimate of 0.597, and the smallest bulk ESS was 82 after 200000 iterations). For a covariance group it is a power, |C|^-(d+1)/2 as C nears a singular matrix (a correlation of +-1), and the posterior has infinite mass on that boundary even when the data identify every component: a chain that gets there stays. Measured with animal() and indirect() in one group on pens of 2 to 8 (validation/indirect_gibbs_recovery.R): 4 of 30 data sets of 1600 records had a chain held at a correlation of +-1 (the smallest eigenvalue of the correlation matrix 1.2e-8 to 2.0e-8, against REML correlations of 0.53 to 0.71, none SINGULAR), still there at 200000 iterations; on 200 records, 15 of 20 data sets within 10000 iterations, 12 of them with a REML correlation below 0.99. With the flat prior no chain was held in any of these 50 data sets (on the 20 small ones 1 - |cor| stayed above 9e-3 over the same 10000 iterations). gibbs() warns, adds the reason to message, and mcmc_diagnostics() says "not mixed", when more than 1% of the kept draws of a chain sit within 1e-6 of a singular correlation matrix. With a covariance group on a Gaussian trait prefer "flat" (with enough levels for the group) or a proper prior; with family = "probit" only a proper c(df =, scale =), since the binary likelihood stays positive as the variances grow and the flat prior leaves that posterior improper too (on one data set of validation/indirect_gibbs_probit_recovery.R the flat chains took var(animal) past 30). "flat" is uniform on the variance (chisq(nl - 2); for a group, uniform on C, InvWishart(nl - d - 1)). "uniform_sd" is uniform on the standard deviation (chisq(nl - 1)), the choice Gelman (2006) recommends for a variance with few levels; scalar groups only. c(df = nu0, scale = S0) is the proper scaled inverse chi-square of the BLUPF90 Gibbs programs (OPTION prior), (u'K^-1 u + nu0 S0) / chisq(nl + nu0), and InvWishart(nl + nu0, U'K^-1 U + nu0 S0 I) for a group: S0 is the prior guess of a variance, nu0 its degree of belief. A ridge of 1e-10 times the mean diagonal of the scale matrix keeps the draw off an exactly singular matrix at the start of the chain; it scales with K, so the chain is equivariant in the scale of a declared covariance

family

"gaussian" (default) or "probit", the threshold model for a binary trait by data augmentation (Albert and Chib 1993; Sorensen, Andersen, Gianola and Korsgaard 1995): each record gets a liability drawn from a normal truncated at 0 on the side of its category, the residual variance is FIXED at 1 (it sets the scale of the liability), and the implicit intercept stands for the threshold. The two values are ordered by the rule of model_threshold(), and the second is the upper category: numbers (or logical) by value, a factor by its levels, text that reads as numbers by value; text that is not a number is refused, with the order sorting would give (sorted, "healthy" < "sick" happens to put the disease on top and "diseased" < "healthy" does not), and a factor whose levels are in the alphabetical order factor() gives text runs with a warning that names the upper category. A level of a fixed class effect with every record in the same category is refused, naming the level: with the flat prior of the fixed effects its posterior is improper (Hoeschele and Tier 1995) and the chain drifts instead of converging; merge the level, drop its records, or fit the effect as random. A missing observation (NA, NaN, Inf, -Inf, blank text, the text "NaN" or missing_code) drops its record, as in model_threshold(); it is never one of the two values. Everything else is the chain of the Gaussian case: kernel(), indirect(), genotypes, APY and prior =. It does not use the Laplace approximation of model_threshold(), but it does not rescue a design the binary records leave weak: with indirect() and the components estimated, on pens of 2 to 8 with 3200 records and 4000 animals, no chain of 30 data sets reached a bulk ESS of 400 (ESS 2.5 to 23.3 on the variances after 2 x 6000 iterations, the mean still drifting), where REML on the liability itself recovers the components (validation/indirect_gibbs_probit_recovery.R)

chains

number of independent chains, each from its own seed. The seeds are drawn from R's generator, so set.seed() governs every chain, and the result is the same run in series or in parallel. With more than one, samples, mean, sd, ebv, ebv_sd, b and b_sd pool the chains that ran to the end: the mean and the sd (denominator N - 1) of all their draws stacked. Each chain stays in chain_samples, and the per-chain means and sds of the locations in chain_ebv, chain_ebv_sd, chain_b and chain_b_sd (levels x chains), which mcmc_diagnostics() reads. A chain that stops early (an inadmissible theta or a failed Cholesky inside the chain) is left out of the pooled results and of the diagnostics, and the verdict says so; when no chain ran to the end, the longest one is pooled alone. chain_n holds the draws each chain kept, chains_used the chains pooled and chain_messages what each chain reported. The per-chain location summaries add 2 q M numbers to the fit for q levels and M chains (about 64 MB for a million levels and 4 chains). All chains start from the same point; one chain (the default) leaves the between-chain checks inactive, and 4 is the usual choice for them

cores

R processes that run the chains at the same time (a PSOCK cluster of the parallel package, which also works on Windows); at most chains

verbose

print the chain as it runs, the iteration count every so often, so a long chain is a progress report instead of silence. Defaults to interactive(), live in a session and quiet in scripts and checks. Every fitter also honors Ctrl+C now

residual_draws

how many kept iterations, evenly spaced, store their realized residual e_s / sqrt(s2e_s), with e_s = y - W loc_s (for family = "probit", the drawn liability minus W loc_s, with variance 1). Each such draw is iid N(0, 1) under the model, which is what residual_diagnostics() uses for its pivotal discrepancy checks (Johnson 2007; Yuan and Johnson 2012). They cost 8 * residual_draws * chains bytes per record in the fit (pdm_residuals, one column per draw, the pooled chains side by side), which with cores > 1 also cross the link from each worker: 96 MB per chain for 600,000 records at the default 20. The chains themselves do not change (no extra random draws). 0 stores none

Value

samples matrix (kept iterations x parameters), posterior mean and sd, the posterior mean/sd of every random effect (ebv, ebv_sd) and of every fixed effect (b, b_sd, named term=level as in model(), whose parametrization note applies: dropped columns are in dropped_x and only contrasts compare against a reference-level convention), dense_block and k_level_type as in model(), start, the components the chain starts from (the first chain's, with several), and the convergence diagnostics: diagnostics is mcmc_diagnostics() of the fit, and ess (the multi-chain bulk-ESS), rhat (rank-normalized split R-hat, also for one chain, from its two halves) and geweke (the Geweke z combined over the chains, with the spectral density at zero, on the scale given by the start_scale column of the diagnostics) are its columns for the components; they are NA and the message says why if the diagnostics fail. n_iter, burnin, thin and n_expected (the draws each chain keeps when it runs to the end) record the run, and fixed_components marks the components the model holds (theta_fixed, kernel(fixed =), the probit residual), which the diagnostics report as "fixed". residual is the posterior mean of the realized residual, one value per row of data (NA where the record did not enter), which by linearity is y - X E[b] - Z E[u] (the liability for the probit), fitted and fitted_fixed the posterior means X E[b] + Z E[u] and X E[b] (on the liability scale for the probit), and pdm_residuals the stored draws described in residual_draws; with several chains all four pool the chains pooled for b and ebv. residuals.breeding_fit() and residual_diagnostics() read them, with missing_code as given. b_standardized (internal) is E[b] on the engine's standardized scale (each covariate centred and scaled), pooled the same way, which residual_diagnostics() uses to check the chain against the data

References

Hobert, J.P. & Casella, G. (1996). The effect of improper priors on Gibbs sampling in hierarchical linear mixed models. Journal of the American Statistical Association 91:1461-1473.

Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis 1:515-534.

Hoeschele, I. & Tier, B. (1995). Estimation of variance components of threshold characters by marginal posterior modes and means via Gibbs sampling. Genetics Selection Evolution 27:519-540.

Geman, S. & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE TPAMI 6:721-741.

Wang, C.S., Rutledge, J.J. & Gianola, D. (1993). Genetics Selection Evolution 25:41-62; (1994) 26:91-115.