Kernel Models for Predicting Phenotypes across Multi-Environment Conditions

Fits Bayesian linear mixed models for multi-environment trials (MET) using the kernels produced by get_kernel(). The model is $$y = 1\mu + X\beta + \sum_j u_j + e$$ with \(u_j \sim N(0, K_j \sigma^2_j)\) and \(e \sim N(0, I\sigma^2_e)\). Sampling runs in the spectral basis of each kernel, so every kernel is used through its eigendecomposition rather than as an n x n matrix.

The function returns genotype predictions, per-kernel variance components, broad- and narrow-sense heritabilities, and – with keep_effects = TRUE – posterior effect chains and breeder-facing report tables (per-environment rankings, stability, selection candidates).

Set verbose = TRUE to trace the six execution steps with per-step and total timings; verbose = FALSE (default) is entirely silent.

Usage

kernel_model(
  y,
  data = NULL,
  random = NULL,
  fixed = NULL,
  env,
  gid,
  verbose = FALSE,
  iterations = 1000,
  burnin = 200,
  thining = 10,
  tol = 1e-10,
  R2 = 0.5,
  digits = 4,
  seed = NULL,
  probs = c(0.025, 0.975),
  variance = c("realized", "scale"),
  scale_kernels = TRUE,
  keep_effects = FALSE,
  compute_CI = TRUE,
  R2_split = c("per_kernel", "shared"),
  prior_scale = c("eigen", "diag"),
  varcomp_method = c("postmean", "gibbs")
)

Arguments

y

character. Name of the phenotype column in data.

data

data.frame with the environment, genotype and phenotype columns. Rows must be grouped by environment if any kernel is block-diagonal.

random

list of kernels from get_kernel. Each element is a list carrying $decomp (the spectral decomposition) and $Type ("D" or "BD"); $Kernel is optional and is not required by the sampler. Factored kernels (no $Kernel) and dense kernels may be mixed in one list. Names should follow the KG_ / KE_ / KGE_ convention, which drives the variance-component roles.

fixed

matrix (n x p) of fixed effects, or NULL.

env, gid

character. Environment and genotype column names in data.

verbose

logical or 0/1. TRUE traces the six execution steps with per-step and total timings; FALSE (default) is entirely silent – no cat(), no message(). Results are identical either way.

iterations

integer. Total Gibbs iterations. Default 1000.

burnin

integer. Iterations discarded before collecting draws. Default 200.

thining

integer. Keep every thining-th post-burn-in draw. Default 10.

tol

numeric. Eigenvalue cutoff. Must match get_kernel(tol =) or the cached decomposition is rejected and recomputed. Default 1e-10.

R2

numeric in (0, 1). Prior proportion of phenotypic variance explained by the random terms. Default 0.5.

digits

integer. Rounding for the variance-component tables.

seed

integer or NULL. Sets the RNG state for reproducibility.

probs

length-2 numeric. Quantiles for posterior prediction intervals. Default c(0.025, 0.975).

variance

"realized" (default) or "scale". See Details.

scale_kernels

logical. Normalise each kernel to mean(diag(K)) = 1. Must match get_kernel(scale_kernels =) for the cache to be reused. Default TRUE.

keep_effects

logical. Retain per-draw effect chains. Required for varcomp_method = "postmean" and for the breeder-facing report tables ($predictions, $genotype_effects, $env_summary). Costs O(k * draws * n) memory. Default FALSE.

compute_CI

logical. Compute posterior prediction intervals for yHat. Requires storing effect chains. Default TRUE.

R2_split

"per_kernel" (default) or "shared". Controls how the prior variance R2 is distributed across kernels. See Prior parameterisation.

prior_scale

"eigen" (default) or "diag". Whether the prior scale is divided by the mean eigenvalue of each kernel, making kernels of differing rank comparable. See Prior parameterisation.

varcomp_method

"postmean" (default) or "gibbs". The estimator used for $varcomp. "postmean" needs keep_effects = TRUE; when unavailable it falls back to "gibbs" (reported in the trace when verbose = TRUE).

Value

An object of class "kernel_model", a list with:

yHat

numeric vector of fitted/predicted values, length n.

yHat.CI

matrix of posterior prediction intervals (NULL unless compute_CI = TRUE).

varE

posterior mean residual variance.

random

per-kernel posterior effects and variances.

fit

the full sampler object: chains, kernel names, MCMC settings and – with keep_effects = TRUE – Uchain. (Named $BGGE in versions before v1.2.3.)

meta

the aligned env / gid / y data frame.

VarComp

variance-component table with confidence intervals.

varcomp

validated components and heritabilities, class "km_varcomp"; see varcomp_summary.

varcomp_gibbs

the same computed by the Gibbs-chain estimator, for comparison.

runtime

sampler wall time in seconds.

With keep_effects = TRUE it additionally carries predictions, genetic_cor_env, genotype_effects, env_summary and variance_proportions.

Details

Fits $$y = 1\mu + X\beta + \sum_j u_j + e$$ with \(u_j \sim N(0, K_j \sigma^2_{u_j})\) and \(e \sim N(0, I\sigma^2_e)\), by Gibbs sampling in the spectral basis of each kernel. Writing \(K_j = U_j S_j U_j'\), the effects are sampled as \(u_j = U_j b_j\) with \(b_j \sim N(0, S_j \sigma^2_{u_j})\), so the cost per iteration scales with the retained rank of each kernel, not with \(n^2\). Trimming the eigenvalue tail via get_kernel(keep_var =) is therefore the most direct way to speed up sampling.

variance = "scale" reports the raw kernel-scale parameter sigb[j] (original BGGE behaviour, inflated for low-rank incidence kernels such as KE_). variance = "realized" reports var(u_j) averaged over draws – what most users mean by a variance component. The raw scale is always returned as $sigb.

Original Version

Costa-Neto et al (2021)

EnvRtype v.0.1.3, Sep 2026

Performance changes vs. the original

• Cached eigendecomposition. If random carries $decomp (attached by get_kernel(decompose = TRUE)) and the cache matches tol / scale_kernels / n, eigen() is skipped entirely. Cross-validation folds and multiple chains reuse one decomposition.

• No stored transposes. The original kept both U and t(U) per kernel (doubling eigenvector memory) and computed crossprod(tU, b). Now only U is stored and the back-transform is U %*% b.

• Optional chain storage. Full per-draw effect matrices (k x nCum x n doubles) are allocated only when needed (keep_effects = TRUE or yHat.CI). Otherwise running sums reduce memory from O(k*nCum*n) to O(k*n).

• Vectorised post-processing. Row-wise apply() over draw matrices replaced by rowSums/rowMeans identities; sweep() replaced by recycled arithmetic.

• Scalar chains in a matrix rather than a nested list, removing per-iteration copy-on-modify of the chain container.

Prior parameterisation

Two defects in the original parameterisation are addressed by arguments:

R2_split. The original gave every kernel the full R2 share, so with nk kernels the priors collectively assert nk * R2 of the phenotypic variance – with RNMDs (4 kernels, R2 = 0.5) that is twice var(y). "per_kernel" (default) uses R2 / nk so total prior mass stays at R2 regardless of model size, keeping model comparisons on equal footing. "shared" restores the original behaviour.

prior_scale. Kernels normalised to mean(diag(K)) = 1 are equal in scale but not in effective dimension: a rank-5 and a rank-750 kernel with the same mean diagonal carry very different total variance. "eigen" (default) divides the prior scale by mean(s) (= trace / rank), making priors comparable across kernels of differing rank.

References

Costa-Neto G, Galli G, Carvalho HF, Crossa J, Fritsche-Neto R (2021). EnvRtype: a software to interplay enviromics and quantitative genomics in agriculture. G3, 11(4), jkab040.

Granato I, Cuevas J, Luna-Vazquez F, Crossa J, Montesinos-Lopez O, Burgueno J, Fritsche-Neto R (2018). BGGE: a new package for genomic-enabled prediction incorporating genotype x environment interaction models. G3, 8(9), 3039-3047.

See also

get_kernel, varcomp_summary, kernel_model_mc, kernel_cv

Author

Germano Costa Neto. Gibbs sampler adapted from Granato et al. (2018), BGGE package. Refactored version.

Examples

# \donttest{
data("maizeYield"); data("maizeG"); data("maizeWTH")

ECs <- W_matrix(env.data = maizeWTH[maizeWTH$daysFromStart < 100, ],
                var.id = c("FRUE", "PETP", "SRAD", "T2M_MAX"),
                statistic = "mean")
#> ---------------------------------------------------------------
#> W_matrix -- builds the environmental covariable (W) matrix
#> Last updated at 09/27/2026
#> ---------------------------------------------------------------
#>   - summarising weather data into environmental covariables
#>   - centring, scaling and quality-controlling W
K <- get_kernel(K_G = list(G = maizeG),
                K_E = list(W = env_kernel(env.data = ECs)[[2]]),
                data = maizeYield, model = "RNMM",
                env = "env", gid = "gid", y = "value")
#> ---------------------------------------------------------------
#> env_kernel -- builds environmental relatedness kernels
#> Last updated at 09/27/2026
#> ---------------------------------------------------------------
#>   - computing a single environmental kernel

## silent fit
fit <- kernel_model(y = "value", data = maizeYield, random = K,
                    env = "env", gid = "gid",
                    iterations = 5000, burnin = 1000)
fit
#> <kernel_model> Bayesian G + GxE fit (BGGE sampler, Granato et al. 2018)
#>   observations : 750 
#>   kernels      : KE_W, KG_G, KGE_GW 
#>   iterations   : 5000 (burn-in 1000 , thinning 10 )
#>   residual var : 0.2967 
#>   runtime (s)  : 2.61 
#>   variance components:
#>             Type        K    Var CI_lower CI_upper SD.var
#>  Environment (E)        W 0.5952   0.5477   0.6495 0.0331
#>     Genotype (G)        G 0.0735   0.0670   0.0810 0.0114
#>              GxE       GW 0.0471   0.0429   0.0520 0.0081
#>         Residual Residual 0.2967   0.2704   0.3273 0.0172
#>   (re-fit with keep_effects = TRUE to include the report tables)
#>   validated varcomp (method = 'gibbs'):
#>     h2 entry-mean 0.0730 | h2 plot 0.0730 | rg NA
#>     ! q < 20 environments: sigma2_GE only approximate
fit$varcomp
#> ------------------------------------------------------------
#> Variance components  (method = 'gibbs')
#>   genotypes p = 150   environments q = 5
#> ------------------------------------------------------------
#>        Component Variance  Prop
#>     Genotype (G)   0.0735 0.073
#>  Environment (E)       NA    NA
#>              GxE       NA    NA
#>         Residual   0.9333 0.927
#>   Total (var(y))   1.0068 1.000
#> ------------------------------------------------------------
#>   h2 (entry-mean) : 0.0730
#>   h2 (plot-level) : 0.0730
#>   rg  vG/(vG+vGE) : NA
#> ------------------------------------------------------------
#>   ! q < 20 environments: sigma2_GE is only approximate.
#>     Validated GxE bias: +0.90 at q=8 vs +0.09 at q=20.

## traced fit, with effect chains for the report tables
fit2 <- kernel_model(y = "value", data = maizeYield, random = K,
                     env = "env", gid = "gid",
                     iterations = 5000, burnin = 1000,
                     keep_effects = TRUE, verbose = TRUE)
#> 
#> =========================================================
#>  kernel_model: Bayesian genomic prediction for multi-
#>    environment trials. Fits y = 1mu + Xb + sum_j u_j + e
#>    with u_j ~ N(0, K_j sigma2_j) by Gibbs sampling in the
#>    spectral basis of each kernel, returning predictions,
#>    variance components and heritabilities.
#>  Original Version: Costa-Neto et al (2021)
#>  EnvRtype v.0.1.3, Sep 2026
#> =========================================================
#>  Gibbs sampler after Granato et al (2018) G3
#> ---------------------------------------------------------
#>   [1/6] Validating inputs ...                             0.00s  750 obs, 3 kernels
#>   [2/6] Assembling response and design ...                0.00s  5 env, 150 gid
#>   [3/6] Spectral decomposition ...                        0.00s  all kernels pre-factored; eigen skipped
#>   [4/6] Computing prior scales ...                        0.00s  R2 = 0.50 split 'per_kernel', prior_scale 'eigen'
#>   [5/6] Gibbs sampling (5000 iter) ...                    2.56s  400 draws kept
#>   [6/6] Variance components and heritability ...          0.00s  method 'postmean', h2 = 0.121
#> ---------------------------------------------------------
#>  Total running time: 2.62s  (sampler 2.56s, 98%)
#>  Residual variance  : 0.2984
#> =========================================================
#> 
head(fit2$genotype_effects)
#>      gid mean_effect min_effect max_effect     range stability_sd   post_sd
#> 92  G092   0.5700298  0.2636566  0.7998444 0.5361878    0.1951406 0.1756836
#> 132 G132   0.5269830  0.3201177  0.6947423 0.3746245    0.1390777 0.1728714
#> 71  G071   0.4861233  0.2374689  0.8138733 0.5764044    0.2185658 0.1661643
#> 65  G065   0.4861134  0.2374571  0.8138807 0.5764236    0.2185805 0.1661877
#> 14  G014   0.4671437  0.3349233  0.7488440 0.4139206    0.1641661 0.1802682
#> 103 G103   0.4626073  0.3660940  0.7462059 0.3801119    0.1594609 0.1742715

## leave-one-environment-out: kernels do not depend on y, so build once
for (e in unique(maizeYield$env)) {
  tr <- maizeYield; tr$value[tr$env == e] <- NA
  f <- kernel_model(y = "value", data = tr, random = K,
                    env = "env", gid = "gid", iterations = 2000, burnin = 500)
}
# }