2) Example of GMA-SMA models for varietal mixtures in microplots
T. Flutre
14/10/2026
Abstract
This document aims at simulating data with GMA-SMA models, fitting them with various softwares, and checking their estimates.1 Preamble
Dependencies:
suppressPackageStartupMessages(library(plantmix))
suppressPackageStartupMessages(library(lme4))
2 Overview
In this vignette, varietal mixtures will be assembled from a panel of 200 genotypes. Only 1000 mixtures will be observed, among all 1.99^{4} possible ones, i.e., approximately 5%. This incomplete design will be optimized so that it is balanced: all genotype will be observed in the same number of mixtures. Phenotypes will then be simulated organized into a field trial in a randomized complete block design, according to all six GMA-SMA models that can be conceived: 1, 2, 2’, 2’‘, 3 and 3’. See the article by Forst et al (2019), especially for the first three models, as well as the introductory vignette.
3 Simulate data
set.seed(12345)3.1 Simulate the genotypes
3.1.1 Set the dimensions
nbGenos <- 200
levGenos <- sprintf(
fmt = paste0("g%0", floor(log10(nbGenos)) + 1, "i"),
1:nbGenos
)
nbSnps <- 1000
levSnps <- sprintf(
fmt = paste0("s%0", floor(log10(nbSnps)) + 1, "i"),
1:nbSnps
)3.1.2 SNP genotypes
nb_pops <- 10
weak_div_pops <- diag(nb_pops)
weak_div_pops[upper.tri(weak_div_pops)] <- 0.9
weak_div_pops[lower.tri(weak_div_pops)] <- weak_div_pops[upper.tri(weak_div_pops)]
tmp <- rep(nbGenos / nb_pops, nb_pops - 1)
tmp <- c(tmp, nbGenos - sum(tmp))
snpGenos <- simulGenosDoseStruct(
nb_genos = tmp,
nb_snps = nbSnps,
div_pops = weak_div_pops,
geno_IDs = levGenos,
snp_IDs = levSnps
)
dim(snpGenos)
## [1] 200 1000
snpGenos[1:3, 1:4]
## s0001 s0002 s0003 s0004
## g001 1 1 1 1
## g002 1 0 1 1
## g003 0 1 2 03.1.3 Additive genetic relationships
For simplicity, the first estimator of VanRaden (2008) is used, that assumes linkage equilibrium and Hardy-Weinberg equilibrium.
GRM <- estimGRM(snpGenos)
GRM <- as.matrix(Matrix::nearPD(GRM)$mat)image(Matrix(GRM), main = "GRM")
hist(diag(GRM), main = "GRM")
hist(GRM[upper.tri(GRM)], main = "GRM")
3.2 Simulate the phenotypes
3.2.1 Incomplete, yet balanced design
nbMixes <- 1000
design <- getDesignBinaryVarMix(levGenos, nbMixes = nbMixes)
tmp <- getMixturesPerGeno(getMixtureList(design$combs))
table(sapply(tmp, length)) # each genotype is in the same nb of mixtures -> balanced design
##
## 9 10 11
## 54 92 54
3.2.2 Field trial
monoStands <- paste(levGenos, levGenos, sep = "_")
pairs <- design$combs
mixStands <- paste(pairs[, 1], pairs[, 2], sep = "_")
IDs <- c(monoStands, mixStands)
nbBlocks <- 3
blocks <- LETTERS[1:nbBlocks]
dat <- do.call(rbind, lapply(blocks, function(block) {
data.frame(
ID = IDs,
block = block,
stringsAsFactors = TRUE
)
}))3.2.3 Design matrices
listContr <- list(block = "contr.sum")
X <- model.matrix(~ 1 + block, data = dat, contrasts = listContr)
allZ <- list()
allZ$GMA <- mkZGMA(dat, col = "ID", sep = "_")
system.time(
allZ <- append(
allZ,
mkAllZSMA(dat, col = "ID", sep = "_", verbose = 1)
)
)
## [1] "model 2'"
## [1] "model 2"
## [1] "model 2''"
## [1] "model 3"
## [1] "model 3'"
## user system elapsed
## 29.428 1.464 30.895
sapply(allZ, dim)
## GMA SMA_mod2p SMA_mod2 SMA_mod2pp_ij SMA_mod2pp_ii SMA_mod3 SMA_mod3p_ij
## [1,] 3600 3600 3600 3600 3600 3600 3600
## [2,] 200 1000 1200 1000 200 1200 1000
## SMA_mod3p_ii
## [1,] 3600
## [2,] 2003.2.4 Parameters
Several packages will be tested below, most notably lme4
and TMB. However, lme4 does not natively take
a GRM as input. Hence the GRM will be replaced by the identity matrix to
ensure a fair comparison between packages.
truth <- list(
"intercept" = 100,
"var_GMA" = 10,
"var_SMA" = 2,
"var_SMA_ij" = 1.5,
"var_SMA_ii" = 0.8,
"var_error" = 1
)
set.seed(1234)
truth[["blockEffs"]] <- sample(x = c(-1, 1), size = nbBlocks - 1, replace = TRUE) *
rnorm(n = nbBlocks - 1, mean = 3, sd = 2)
GRM <- diag(nbGenos)
truth[["GMAs"]] <- MASS::mvrnorm(
n = 1, mu = rep(0, nbGenos),
Sigma = truth$var_GMA * GRM
)
truth[["SMAs"]] <- rnorm(n = length(IDs), mean = 0, sd = sqrt(truth$var_SMA))
truth[["SMAs_ij"]] <- rnorm(n = length(mixStands), mean = 0, sd = sqrt(truth$var_SMA_ij))
truth[["SMAs_ii"]] <- rnorm(n = length(monoStands), mean = 0, sd = sqrt(truth$var_SMA_ii))
truth[["errors"]] <- rnorm(n = nrow(dat), mean = 0, sd = sqrt(truth$var_error))3.2.5 Phenotypes
y1 <- X %*% c(truth$intercept, truth$blockEffs) +
allZ$GMA %*% truth$GMAs +
truth$errors
dat$pheno1 <- y1[, 1]
y2 <- X %*% c(truth$intercept, truth$blockEffs) +
allZ$GMA %*% truth$GMAs +
allZ$SMA_mod2 %*% truth$SMAs +
truth$errors
dat$pheno2 <- y2[, 1]
y3 <- X %*% c(truth$intercept, truth$blockEffs) +
allZ$GMA %*% truth$GMAs +
allZ$SMA_mod3 %*% truth$SMAs +
truth$errors
dat$pheno3 <- y3[, 1]
y2p <- X %*% c(truth$intercept, truth$blockEffs) +
allZ$GMA %*% truth$GMAs +
allZ$SMA_mod2p %*% truth$SMAs_ij +
truth$errors
dat$pheno2p <- y2p[, 1]
y2pp <- X %*% c(truth$intercept, truth$blockEffs) +
allZ$GMA %*% truth$GMAs +
allZ$SMA_mod2pp_ij %*% truth$SMAs_ij +
allZ$SMA_mod2pp_ii %*% truth$SMAs_ii +
truth$errors
dat$pheno2pp <- y2pp[, 1]
y3p <- X %*% c(truth$intercept, truth$blockEffs) +
allZ$GMA %*% truth$GMAs +
allZ$SMA_mod3p_ij %*% truth$SMAs_ij +
allZ$SMA_mod3p_ii %*% truth$SMAs_ii +
truth$errors
dat$pheno3p <- y3p[, 1]4 Infer the parameters
Function fitGMASMA allows to fit the models with various
packages. In this vignette, only lme4 and TMB
will be used and, as you can see below, both return the same estimates.
(To save time, MM4LMM is commented.)
pkgs <- c("lme4", "TMB") # , "MM4LMM")
if ("MM4LMM" %in% pkgs) {
suppressPackageStartupMessages(library(MM4LMM))
}
runAllPkgs <- function(pkgs, form, dat, listZ, listVCov, listContr, ...) {
## mcMap(function(pkg) {
Map(function(pkg) {
print(paste0("fit with ", pkg, "..."))
fitGMASMA(form, dat, listZ, pkg, listVCov, listContr)
st <- system.time(
try(fit <- fitGMASMA(form, dat, listZ, pkg, listVCov, listContr, REML = TRUE, ...))
)
print(st)
fit
}, pkgs) # , mc.cores = nbCores)
}4.1 Model 1
4.1.1 Fit
form <- pheno1 ~ 1 + block
listZ <- list("GMA" = allZ$GMA)
listVCov <- list("GMA" = GRM)
listContr <- list(block = "contr.sum")
fits1 <- runAllPkgs(pkgs, form, dat, listZ, listVCov, listContr)
## [1] "fit with lme4..."
## user system elapsed
## 0.501 0.000 0.501
## [1] "fit with TMB..."
## user system elapsed
## 3.824 0.000 3.8244.1.2 Check
print("lme4:")
## [1] "lme4:"
tmp <- as.data.frame(VarCorr(fits1$lme4))
tmp <- setNames(tmp$vcov, tmp$grp)
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_error")]),
estim = tmp[c("GMA", "Residual")]
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10 10.294633 2.946325
## var_error 1 0.957673 -4.232704
print("TMB:")
## [1] "TMB:"
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_error")]),
estim = do.call(c, fits1$TMB$report[c("var_GMA", "var_err")])
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10 10.294632 2.946325
## var_error 1 0.957673 -4.232704
if ("MM4LMM" %in% names(fits1)) {
print(fits1$MM4LMM$Sigma2)
}4.2 Model 2
4.2.1 Fit
form <- pheno2 ~ 1 + block
listZ <- list(
"GMA" = allZ$GMA,
"SMA" = allZ$SMA_mod2
)
listVCov <- list(
"GMA" = GRM,
"SMA" = diag(ncol(listZ$SMA))
)
listContr <- list(block = "contr.sum")
fits2 <- runAllPkgs(pkgs, form, dat, listZ, listVCov, listContr)
## [1] "fit with lme4..."
## user system elapsed
## 1.497 0.001 1.497
## [1] "fit with TMB..."
## user system elapsed
## 11.226 0.036 11.2634.2.2 Check
print("lme4:")
## [1] "lme4:"
tmp <- as.data.frame(VarCorr(fits2$lme4))
tmp <- setNames(tmp$vcov, tmp$grp)
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA", "var_error")]),
estim = tmp[c("GMA", "SMA", "Residual")]
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10 10.1324793 1.3247926
## var_SMA 2 2.0113789 0.5689449
## var_error 1 0.9453839 -5.4616064
print("TMB:")
## [1] "TMB:"
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA", "var_error")]),
estim = do.call(c, fits2$TMB$report[c("var_GMA", "var_SMA1", "var_err")])
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10 10.1324647 1.3246471
## var_SMA 2 2.0114019 0.5700966
## var_error 1 0.9453792 -5.4620796
if ("MM4LMM" %in% names(fits2)) {
print(fits2$MM4LMM$Sigma2)
}4.3 Model 3
4.3.1 Fit
form <- pheno3 ~ 1 + block
listZ <- list(
"GMA" = allZ$GMA,
"SMA" = allZ$SMA_mod3
)
listVCov <- list(
"GMA" = GRM,
"SMA" = diag(ncol(listZ$SMA))
)
listContr <- list(block = "contr.sum")
fits3 <- runAllPkgs(pkgs, form, dat, listZ, listVCov, listContr)
## [1] "fit with lme4..."
## user system elapsed
## 2.237 0.000 2.238
## [1] "fit with TMB..."
## user system elapsed
## 11.189 0.064 11.2534.3.2 Check
print("lme4:")
## [1] "lme4:"
tmp <- as.data.frame(VarCorr(fits3$lme4))
tmp <- setNames(tmp$vcov, tmp$grp)
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA", "var_error")]),
estim = tmp[c("GMA", "SMA", "Residual")]
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10 9.9758558 -0.2414422
## var_SMA 2 2.0811855 4.0592759
## var_error 1 0.9443071 -5.5692907
print("TMB:")
## [1] "TMB:"
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA", "var_error")]),
estim = do.call(c, fits3$TMB$report[c("var_GMA", "var_SMA1", "var_err")])
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10 9.9758564 -0.2414362
## var_SMA 2 2.0811856 4.0592811
## var_error 1 0.9443071 -5.5692918
if ("MM4LMM" %in% names(fits3)) {
print(fits3$MM4LMM$Sigma2)
}4.4 Model 2’
4.4.1 Fit
form <- pheno2p ~ 1 + block
listZ <- list(
"GMA" = allZ$GMA,
"SMA" = allZ$SMA_mod2p
)
listVCov <- list(
"GMA" = GRM,
"SMA" = diag(ncol(listZ$SMA))
)
listContr <- list(block = "contr.sum")
fits2p <- runAllPkgs(pkgs, form, dat, listZ, listVCov, listContr)
## [1] "fit with lme4..."
## Warning in split.default(x, g): data length is not a multiple of split variable
## Warning in split.default(seq_along(x), f, drop = drop, ...): data length is not
## a multiple of split variable
## Warning in split.default(x, g): data length is not a multiple of split variable
## Warning in split.default(seq_along(x), f, drop = drop, ...): data length is not
## a multiple of split variable
## user system elapsed
## 1.448 0.000 1.447
## [1] "fit with TMB..."
## user system elapsed
## 9.622 0.068 9.6914.4.2 Check
print("lme4:")
## [1] "lme4:"
tmp <- as.data.frame(VarCorr(fits2p$lme4))
tmp <- setNames(tmp$vcov, tmp$grp)
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA", "var_error")]),
estim = tmp[c("GMA", "SMA", "Residual")]
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10 10.428119 4.281193
## var_SMA 2 1.544762 -22.761905
## var_error 1 0.946256 -5.374396
print("TMB:")
## [1] "TMB:"
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA", "var_error")]),
estim = do.call(c, fits2p$TMB$report[c("var_GMA", "var_SMA1", "var_err")])
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10 10.428119 4.281189
## var_SMA 2 1.544761 -22.761929
## var_error 1 0.946256 -5.374396
if ("MM4LMM" %in% names(fits2p)) {
print(fits2p$MM4LMM$Sigma2)
}4.5 Model 2’’
4.5.1 Fit
form <- pheno2pp ~ 1 + block
listZ <- list(
"GMA" = allZ$GMA,
"SMA_ij" = allZ$SMA_mod2pp_ij,
"SMA_ii" = allZ$SMA_mod2pp_ii
)
listVCov <- list(
"GMA" = GRM,
"SMA_ij" = diag(ncol(listZ$SMA_ij)),
"SMA_ii" = diag(ncol(listZ$SMA_ii))
)
listContr <- list(block = "contr.sum")
fits2pp <- runAllPkgs(pkgs, form, dat, listZ, listVCov, listContr)
## [1] "fit with lme4..."
## Warning in split.default(x, g): data length is not a multiple of split variable
## Warning in split.default(seq_along(x), f, drop = drop, ...): data length is not
## a multiple of split variable
## Warning in split.default(x, g): data length is not a multiple of split variable
## Warning in split.default(seq_along(x), f, drop = drop, ...): data length is not
## a multiple of split variable
## user system elapsed
## 4.454 0.000 4.460
## [1] "fit with TMB..."
## user system elapsed
## 11.405 0.052 11.4604.5.2 Check
print("lme4:")
## [1] "lme4:"
tmp <- as.data.frame(VarCorr(fits2pp$lme4))
tmp <- setNames(tmp$vcov, tmp$grp)
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA_ij", "var_SMA_ii", "var_error")]),
estim = tmp[c("GMA", "SMA_ij", "SMA_ii", "Residual")]
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10.0 10.1734023 1.734023
## var_SMA_ij 1.5 1.5577029 3.846860
## var_SMA_ii 0.8 0.8390411 4.880133
## var_error 1.0 0.9453835 -5.461654
print("TMB:")
## [1] "TMB:"
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA_ij", "var_SMA_ii", "var_error")]),
estim = do.call(c, fits2pp$TMB$report[c("var_GMA", "var_SMA1", "var_SMA2", "var_err")])
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10.0 10.1733882 1.733882
## var_SMA_ij 1.5 1.5577032 3.846877
## var_SMA_ii 0.8 0.8390382 4.879774
## var_error 1.0 0.9453838 -5.461621
if ("MM4LMM" %in% names(fits2pp)) {
print(fits2pp$MM4LMM$Sigma2)
}4.6 Model 3’
4.6.1 Fit
form <- pheno3p ~ 1 + block
listZ <- list(
"GMA" = allZ$GMA,
"SMA_ij" = allZ$SMA_mod3p_ij,
"SMA_ii" = allZ$SMA_mod3p_ii
)
listVCov <- list(
"GMA" = GRM,
"SMA_ij" = diag(ncol(listZ$SMA_ij)),
"SMA_ii" = diag(ncol(listZ$SMA_ii))
)
listContr <- list(block = "contr.sum")
fits3p <- runAllPkgs(pkgs, form, dat, listZ, listVCov, listContr)
## [1] "fit with lme4..."
## Warning in split.default(x, g): data length is not a multiple of split variable
## Warning in split.default(seq_along(x), f, drop = drop, ...): data length is not
## a multiple of split variable
## Warning in split.default(x, g): data length is not a multiple of split variable
## Warning in split.default(seq_along(x), f, drop = drop, ...): data length is not
## a multiple of split variable
## user system elapsed
## 3.792 0.000 3.792
## [1] "fit with TMB..."
## user system elapsed
## 13.325 0.164 13.4904.6.2 Check
print("lme4:")
## [1] "lme4:"
tmp <- as.data.frame(VarCorr(fits3p$lme4))
tmp <- setNames(tmp$vcov, tmp$grp)
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA_ij", "var_SMA_ii", "var_error")]),
estim = tmp[c("GMA", "SMA_ij", "SMA_ii", "Residual")]
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10.0 9.9892439 -0.1075612
## var_SMA_ij 1.5 1.5390438 2.6029207
## var_SMA_ii 0.8 0.9997089 24.9636064
## var_error 1.0 0.9440616 -5.5938365
print("TMB:")
## [1] "TMB:"
checks <- data.frame(
truth = do.call(c, truth[c("var_GMA", "var_SMA_ij", "var_SMA_ii", "var_error")]),
estim = do.call(c, fits3p$TMB$report[c("var_GMA", "var_SMA1", "var_SMA2", "var_err")])
)
checks$nBE <- normBiasError(checks$estim, checks$truth)
print(checks)
## truth estim nBE
## var_GMA 10.0 9.9892539 -0.1074613
## var_SMA_ij 1.5 1.5390442 2.6029496
## var_SMA_ii 0.8 0.9997069 24.9633674
## var_error 1.0 0.9440616 -5.5938432
if ("MM4LMM" %in% names(fits3p)) {
print(fits3p$MM4LMM$Sigma2)
}5 Appendix
t1 <- proc.time()
t1 - t0
## user system elapsed
## 191.285 2.615 193.915
print(sessionInfo(), locale = FALSE)
## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Debian GNU/Linux 12 (bookworm)
##
## Matrix products: default
## BLAS: /usr/local/R/4.5.1/lib/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0 LAPACK version 3.11.0
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plantmix_1.0.4 lme4_1.1-37 Matrix_1.7-3 ggplot2_3.5.2 knitr_1.50
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 jsonlite_2.0.0 dplyr_1.1.4 compiler_4.5.1
## [5] Rcpp_1.1.2 tidyselect_1.2.1 dichromat_2.0-0.1 jquerylib_0.1.4
## [9] splines_4.5.1 scales_1.4.0 boot_1.3-31 yaml_2.3.10
## [13] fastmap_1.2.0 lattice_0.22-7 R6_2.6.1 generics_0.1.4
## [17] igraph_2.1.4 rbibutils_2.3 MASS_7.3-65 tibble_3.3.0
## [21] nloptr_2.2.1 minqa_1.2.8 TMB_1.9.25 bslib_0.9.0
## [25] pillar_1.11.0 RColorBrewer_1.1-3 rlang_1.1.7 cachem_1.1.0
## [29] xfun_0.52 sass_0.4.10 cli_3.6.6 withr_3.0.2
## [33] magrittr_2.0.5 Rdpack_2.6.4 digest_0.6.39 grid_4.5.1
## [37] lifecycle_1.0.5 nlme_3.1-168 reformulas_0.4.4 vctrs_0.6.5
## [41] evaluate_1.0.5 glue_1.8.1 farver_2.1.2 rmarkdown_2.29
## [45] tools_4.5.1 pkgconfig_2.0.3 htmltools_0.5.8.1