3) Example of DBV-SBV models for species mixtures in a tester-based design
Preamble
Dependencies:
Simulate data
Simulate the genotypes
Set the dimensions
levSpecies <- c("S1", "S2")
nbGenos <- c("S1" = 500, "S2" = 500)
levGenos <- list(
"S1" = sprintf(
fmt = paste0("gS1_%0", floor(log10(nbGenos["S1"])) + 1, "i"),
1:nbGenos["S1"]
),
"S2" = sprintf(
fmt = paste0("gS2_%0", floor(log10(nbGenos["S2"])) + 1, "i"),
1:nbGenos["S2"]
)
)
nbSnps <- c("S1" = 1000, "S2" = 1000)
levSnps <- list(
"S1" = sprintf(
fmt = paste0("sS1_%0", floor(log10(nbSnps["S1"])) + 1, "i"),
1:nbSnps["S1"]
),
"S2" = sprintf(
fmt = paste0("sS2_%0", floor(log10(nbSnps["S2"])) + 1, "i"),
1:nbSnps["S2"]
)
)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)]
snpGenos <- lapply(levSpecies, function(species) {
tmp <- rep(nbGenos[species] / nb_pops, nb_pops - 1)
tmp <- c(tmp, nbGenos[species] - sum(tmp))
simulGenosDoseStruct(
nb_genos = tmp,
nb_snps = nbSnps[species],
div_pops = weak_div_pops,
geno_IDs = levGenos[[species]],
snp_IDs = levSnps[[species]]
)
})
names(snpGenos) <- levSpecies
sapply(snpGenos, dim)
## S1 S2
## [1,] 500 500
## [2,] 1000 1000
snpGenos$S1[1:3, 1:4]
## sS1_0001 sS1_0002 sS1_0003 sS1_0004
## gS1_001 0 1 1 1
## gS1_002 1 1 2 1
## gS1_003 1 0 1 1
table(snpGenos$S1)
##
## 0 1 2
## 154493 200572 144935Additive genetic relationships
For simplicity, the first estimator of VanRaden (2008) is used, that assumes linkage equilibrium and Hard-Weinberg equilibrium.
snpAFs <- lapply(snpGenos, function(M) {
colSums(M) / (2 * nrow(M))
})
GRMs_vr <- lapply(levSpecies, function(species) {
GRM <- estimGRM(snpGenos[[species]], snpAFs[[species]])
as.matrix(Matrix::nearPD(GRM)$mat)
})
names(GRMs_vr) <- levSpecies
Simulate the phenotypes
As in Salomon et al (2026), the design will be incomplete (sparse) but balanced, and tester-based, involving:
many genotypes of the first species, the focal one, whose breeding values will be statistically modeled as random, with kinship matrix \(K\);
and a small number of genotypes of the second species, the tester one, whose breeding values will be statistically modeled as fixed.
The yield data are generated according to the following equations:
intercrops: \(Y_{IC} = X_{IC} B_{IC} + Z_{DS_f} BV_f + Z_{D{\times}S} \, DBV{\times}SBV + E_{IC}\)
sole crops:
focal species: \(y_{SC_f} = X_{SC_f} \beta_{SC_f} + Z_{D_f} DBV_f + Z_{D_f} SIGV_f + e_{SC_f}\) where \(\beta_f\) only includes the contrasts for the “block” explanatory factor
tester species: \(y_{SC_t} = X_{SC_t} \beta_{SC_t} + e_{SC_t}\) where \(\beta_t\) includes the contrasts for the “block” and “DBV” explanatory factors (the “SIGV” explanatory factor for the tester species is ignored)
The parameter values correspond to cereal-legume mixtures such as winter wheat and pea, inspired from the papers of Moutier et al (2022) and Haug et al (2023).
nbGenosTrial <- c("S1" = 300, "S2" = 2)
levGenosTrial <- lapply(levSpecies, function(species) {
sample(levGenos[[species]], nbGenosTrial[species])
})
names(levGenosTrial) <- levSpecies
GRMs_vr_trial <- list()
GRMs_vr_trial$S1 <- GRMs_vr$S1[levGenosTrial$S1, levGenosTrial$S1]
GRMs_vr_trial$S2 <- diag(nbGenosTrial["S2"]) # diag because modeled as fixed
dimnames(GRMs_vr_trial$S2) <- list(levGenosTrial$S2, levGenosTrial$S2)
set.seed(12345)
out <- simulDBVSBVinter(GRMs_vr_trial) # see ?simulDBVSBVinter for more details
names(out)
## [1] "truth" "datW" "datL" "obsMC"
## [5] "sowingDensities" "props"
tmp <- list2env(out, envir = environment())
Indices of plant mixtures
Several indices can be used to compare sole crops and intercrops. Below some of them are computed using the true breeding values, i.e., with neither block effects nor experimental errors, to give an idea of what the simulated data correspond to.
RYT (LER) and RYP
## Reformat and compute
is_mix <- which(datW$type == "IC")
true_RYTs <- list()
true_RYPs <- list()
for (idx in is_mix) {
true_y1_IC <- as.vector(truth$mu[["S1"]]["IC"]) + datW$true_gen_yield_S1[idx]
true_y2_IC <- as.vector(truth$mu[["S2"]]["IC"]) + datW$true_gen_yield_S2[idx]
g1 <- as.character(datW$geno_S1[idx])
g2 <- as.character(datW$geno_S2[idx])
true_y1_SC <- as.vector(truth$mu[["S1"]]["SC"]) +
datW$true_gen_yield_S1[datW$geno_S1 == g1 &
is.na(datW$geno_S2)]
true_y2_SC <- as.vector(truth$mu[["S2"]]["SC"]) +
datW$true_gen_yield_S2[datW$geno_S2 == g2 &
is.na(datW$geno_S1)]
true_y2_SC <- unique(true_y2_SC)
yields <- data.frame(
SCcrop = c(true_y1_SC, true_y2_SC),
intercrop = c(true_y1_IC, true_y2_IC),
row.names = c(g1, g2)
)
tmp <- LER(yields)
mixId <- paste0(g1, "-", g2)
true_RYTs[[mixId]] <- c(
"RY_S1" = as.numeric(tmp$pLER[1]),
"RY_S2" = as.numeric(tmp$pLER[2]),
"RYT" = tmp$LER
)
true_RYPs[[mixId]] <- c(
"RYP_S1" = true_y1_IC /
(props[["S1"]] * true_y1_SC),
"RYP_S2" = true_y2_IC /
(props[["S2"]] * true_y2_SC)
)
}
true_RYTs <- data.frame(
mix = names(true_RYTs),
geno_S1 = sapply(strsplit(names(true_RYTs), "-"), `[`, 1),
geno_S2 = sapply(strsplit(names(true_RYTs), "-"), `[`, 2),
as.data.frame(do.call(rbind, true_RYTs)),
stringsAsFactors = TRUE
)
str(true_RYTs)
## 'data.frame': 300 obs. of 6 variables:
## $ mix : Factor w/ 300 levels "gS1_001-gS2_301",..: 24 285 118 267 158 234 140 288 96 279 ...
## $ geno_S1: Factor w/ 300 levels "gS1_001","gS1_002",..: 24 285 118 267 158 234 140 288 96 279 ...
## $ geno_S2: Factor w/ 2 levels "gS2_191","gS2_301": 1 1 1 1 1 1 1 1 1 1 ...
## $ RY_S1 : num 0.508 0.511 0.528 0.537 0.382 ...
## $ RY_S2 : num 0.853 0.93 0.808 0.787 1.056 ...
## $ RYT : num 1.36 1.44 1.34 1.32 1.44 ...
summary(true_RYTs[, grep("RY_", colnames(true_RYTs))])
## RY_S1 RY_S2
## Min. :0.2333 Min. :0.7094
## 1st Qu.:0.4570 1st Qu.:0.8527
## Median :0.4875 Median :0.9086
## Mean :0.4874 Mean :0.9061
## 3rd Qu.:0.5234 3rd Qu.:0.9536
## Max. :0.6410 Max. :1.1837
summary(true_RYTs[, grep("RYT", colnames(true_RYTs))])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.236 1.353 1.395 1.394 1.432 1.542
true_RYPs <- data.frame(
mix = names(true_RYPs),
geno_S1 = sapply(strsplit(names(true_RYPs), "-"), `[`, 1),
geno_S2 = sapply(strsplit(names(true_RYPs), "-"), `[`, 2),
as.data.frame(do.call(rbind, true_RYPs)),
stringsAsFactors = TRUE
)
str(true_RYPs)
## 'data.frame': 300 obs. of 5 variables:
## $ mix : Factor w/ 300 levels "gS1_001-gS2_301",..: 24 285 118 267 158 234 140 288 96 279 ...
## $ geno_S1: Factor w/ 300 levels "gS1_001","gS1_002",..: 24 285 118 267 158 234 140 288 96 279 ...
## $ geno_S2: Factor w/ 2 levels "gS2_191","gS2_301": 1 1 1 1 1 1 1 1 1 1 ...
## $ RYP_S1 : num 0.644 0.648 0.669 0.68 0.484 ...
## $ RYP_S2 : num 4.05 4.42 3.84 3.74 5.02 ...
summary(true_RYPs[, grep("RYP", colnames(true_RYPs))])
## RYP_S1 RYP_S2
## Min. :0.2955 Min. :3.370
## 1st Qu.:0.5789 1st Qu.:4.050
## Median :0.6175 Median :4.316
## Mean :0.6174 Mean :4.304
## 3rd Qu.:0.6630 3rd Qu.:4.530
## Max. :0.8119 Max. :5.622
if (FALSE) {
## using the RYT() function
keys <- unique(paste0(datL$focal, " in ", datL$standID))
tmp <- do.call(rbind, strsplit(keys, " in "))
datLavg <- data.frame(
key = keys,
focal = tmp[, 1],
standID = tmp[, 2],
stringsAsFactors = TRUE
)
datLavg$type <- "IC"
datLavg$type[as.character(datLavg$focal) == as.character(datLavg$standID)] <- "SC"
datLavg$focal_species <- "S1"
datLavg$focal_species[grep("^gS2_", datLavg$focal)] <- "S2"
datLavg$prop <- 1
datLavg$prop[datLavg$type == "IC" & datLavg$focal_species == "S1"] <- props["S1"]
datLavg$prop[datLavg$type == "IC" & datLavg$focal_species == "S2"] <- props["S2"]
for (i in 1:nrow(datLavg)) {
idx <- which(datL$focal == datLavg$focal[i] &
datL$standID == datLavg$standID[i])
datLavg$focal_yield[i] <- mean(datL$focal_yield[idx])
}
true_RYTs2 <- RYT(datLavg, "standID", "focal", "prop", "focal_yield")
true_RYTs2 <- droplevels(true_RYTs2[!is.na(true_RYTs2$RYT), ])
true_RYTs2 <- droplevels(true_RYTs2[!duplicated(true_RYTs2$standID), ])
}
## Plot
ggplot(true_RYTs) +
aes(x = RY_S1) +
geom_histogram(color = "white", bins = 30) +
geom_vline(
xintercept = sowingDensities$S1["IC"] /
sowingDensities$S1["SC"],
col = "red", linewidth = 2
) +
labs(
title = "True relative yields (partial land-equivalent ratios) of species 1 for all mixtures",
x = "RY (partial LER) of species 1"
) +
theme_bw()
ggplot(true_RYTs) +
aes(x = RY_S2) +
geom_histogram(color = "white", bins = 30) +
geom_vline(
xintercept = sowingDensities$S2["IC"] /
sowingDensities$S2["SC"],
col = "red", linewidth = 2
) +
labs(
title = "True partial land-equivalent ratio of species 2 for all mixtures",
x = "RY (partial LER) of species 2"
) +
theme_bw()
ggplot(true_RYTs) +
aes(x = geno_S2, y = RYT) +
geom_violin(aes(fill = geno_S2), trim = FALSE, show.legend = FALSE) +
geom_boxplot(width = 0.2) +
labs(
title = "True land-equivalent ratio for all mixtures"
) +
theme_bw()
p <- ggplot(true_RYTs) +
aes(x = RY_S1, y = RY_S2, color = geno_S2)
for (i in seq(0.75, 2, by = 0.25)) {
if (i == 1) {
p <- p + geom_abline(intercept = i, slope = -1, linetype = "solid", color = "black")
} else {
p <- p + geom_abline(intercept = i, slope = -1, linetype = "dotted", color = "black")
}
}
p + geom_abline(intercept = 0, slope = 1, linetype = "dotted", color = "black") +
geom_point(size = 2) +
labs(
title = "True relative yields (RYs, a.k.a. partial LERs)",
x = "relative yield (partial LER) of species 1",
y = "relative yiedl (partial LER) of species 2",
color = "Tester"
) +
## guides(color="none") +
scale_x_continuous(breaks = seq(0, 1.4, by = 0.1)) +
scale_y_continuous(breaks = seq(0, 1.4, by = 0.1)) +
coord_cartesian(xlim = c(0, 1.4), ylim = c(0, 1.4)) +
theme(aspect.ratio = 1) +
theme_bw()
ggplot(true_RYPs) +
aes(x = RYP_S1, y = RYP_S2, color = geno_S2) +
geom_abline(intercept = 0, slope = 1, linetype = "solid", color = "black") +
geom_hline(yintercept = 1) +
geom_vline(xintercept = 1) +
geom_point(size = 2) +
labs(
title = "True relative yields per plant (RYPs)",
x = "RYP of species 1",
y = "RYP of species 2",
color = "Tester"
) +
## guides(color="none") +
scale_x_continuous(breaks = seq(0, 6.5, by = 1)) +
scale_y_continuous(breaks = seq(0, 6.5, by = 1)) +
coord_cartesian(xlim = c(0, 6.5), ylim = c(0, 6.5)) +
theme(aspect.ratio = 1) +
theme_bw()
RYM
## Reformat and compute
tmp <- datW[, c("geno_S1", "geno_S2", "true_yield_S1", "true_yield_S2")]
tmp$ID <- NA
tmp$props <- NA
tmp$true_yield <- NA
## sole crop of species 1:
idx <- which(!is.na(tmp$geno_S1) & is.na(tmp$geno_S2))
tmp$ID[idx] <- as.character(tmp$geno_S1[idx])
tmp$props[idx] <- "1"
tmp$true_yield[idx] <- tmp$true_yield_S1[idx]
## sole crop of species 2:
idx <- which(is.na(tmp$geno_S1) & !is.na(tmp$geno_S2))
tmp$ID[idx] <- as.character(tmp$geno_S2[idx])
tmp$props[idx] <- "1"
tmp$true_yield[idx] <- tmp$true_yield_S2[idx]
## intercrops of species 1 and 2:
idx <- which(!is.na(tmp$geno_S1) & !is.na(tmp$geno_S2))
sep <- "|"
tmp$ID[idx] <- as.character(paste0(
tmp$geno_S1[idx], sep,
tmp$geno_S2[idx]
))
prop1 <- props["S1"]
prop2 <- props["S2"]
stopifnot(prop1 + prop2 == 1)
prop1 <- round(prop1, 2)
prop2 <- 1 - prop1
tmp$props[idx] <- paste0(prop1, sep, prop2)
tmp$true_yield[idx] <- tmp$true_yield_S1[idx] + tmp$true_yield_S2[idx]
stopifnot(all(!is.na(tmp$ID)))
tmp$ID <- factor(tmp$ID)
tmp$props <- factor(tmp$props)
## keep only one yield (the true one) per modality
dupIDs <- table(as.character(tmp$ID))
(dupIDs <- names(dupIDs)[dupIDs > 1])
## [1] "gS2_191" "gS2_301"
for (dupID in dupIDs) {
idx <- which(as.character(tmp$ID) == dupID)
tmp <- droplevels(tmp[-idx[2:length(idx)], ])
}
rm(dupIDs)
tmp <- RYM(tmp,
colIDstand = "ID", colIDcomps = "ID", colProps = "props",
colY = "true_yield", sep = "|"
)
summary(tmp$RYM)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
## 0.7757 0.9207 1.0352 1.0465 1.1540 1.5495 302
## Plot
ggplot(tmp) +
aes(x = RYM) +
geom_histogram(na.rm = TRUE, bins = 30, color = "white") +
geom_vline(xintercept = 1, linewidth = 2) +
geom_vline(xintercept = mean(tmp$RYM, na.rm = TRUE), linewidth = 2, color = "red") +
labs(title = "True relative yields of mixtures (RYMs)") +
theme_bw()
Explore the data
In this section, an exploratory data analysis is done on the data including block effects and experimental errors, so that it can be easily applied on real data, too.
str(datW)
## 'data.frame': 604 obs. of 20 variables:
## $ standID : Factor w/ 602 levels "gS1_001","gS1_001+gS2_301",..: 47 419 513 459 569 235 447 533 315 409 ...
## $ geno_S1 : Factor w/ 300 levels "gS1_001","gS1_002",..: 24 210 257 230 285 118 224 267 158 205 ...
## $ geno_S2 : Factor w/ 2 levels "gS2_191","gS2_301": NA NA NA NA NA NA NA NA NA NA ...
## $ type : Factor w/ 2 levels "SC","IC": 1 1 1 1 1 1 1 1 1 1 ...
## $ type2 : Factor w/ 3 levels "sole_S1","sole_S2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ block : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
## $ x : int 9 16 20 26 4 17 15 5 15 2 ...
## $ y : int 2 2 2 6 1 10 4 9 9 5 ...
## $ plot : Factor w/ 604 levels "10A1","10A10",..: 597 62 110 173 421 70 54 535 59 210 ...
## $ true_gen_yield_S1: num 8.5817 -0.3536 -12.2856 -0.689 0.0692 ...
## $ true_gen_yield_S2: num NA NA NA NA NA NA NA NA NA NA ...
## $ true_yield_S1 : num 70.4 61.4 49.5 61.1 61.9 ...
## $ true_yield_S2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ true_fix_yield_S1: num 61.8 61.8 61.8 61.8 61.8 ...
## $ true_fix_yield_S2: num NA NA NA NA NA NA NA NA NA NA ...
## $ true_rnd_yield_S1: num 8.5817 -0.3536 -12.2856 -0.689 0.0692 ...
## $ true_rnd_yield_S2: num NA NA NA NA NA NA NA NA NA NA ...
## $ yield_S1 : num 73.1 71.2 47.2 61.5 62.7 ...
## $ yield_S2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ tot_yield : num 73.1 71.2 47.2 61.5 62.7 ...
tapply(datW$type, datW$block, table)
## $A
##
## SC IC
## 152 150
##
## $B
##
## SC IC
## 152 150ggplot(datL) +
aes(x = block, y = focal_yield) +
geom_violin(aes(fill = block)) +
geom_boxplot(fill = "white", width = 0.2) +
theme_bw() +
facet_grid(focal_species ~ type)
is_mix <- datW$type == "IC"
subDatW <- droplevels(datW[is_mix, ])
ggplot(subDatW) +
aes(x = yield_S1, y = yield_S2, color = geno_S1, shape = geno_S2) +
geom_abline(intercept = seq(0, 200, by = 10), slope = -1, linetype = "dotted", color = "black") +
geom_point(size = 2) +
labs(
title = "Partial yields in intercrop",
x = "species 1 (in qt.ha-1)", y = "species 2 (in qt.ha-1)",
shape = "Tester (species S2)"
) +
guides(color = "none") +
theme_bw()
Yield map
## Add the empty micro-plots:
coords <- data.frame(
x = rep(sort(unique(datW$x)), each = length(unique(datW$y))),
y = sort(unique(datW$y)),
block = "A",
plot = NA
)
coords$block[coords$x >= min(datW$x[datW$block != "A"])] <- "B"
coords$plot <- paste0(coords$x, coords$block, coords$y)
length(idx <- which(!coords$plot %in% as.character(datW$plot)))
## [1] 16
tmp <- as.data.frame(matrix(nrow = length(idx), ncol = ncol(datW)))
colnames(tmp) <- colnames(datW)
tmp[, c("x", "y", "block", "plot")] <- coords[idx, ]
datWSupp <- rbind(
datW,
as.data.frame(tmp)
)
## Plot
xranges <- do.call(rbind, tapply(datW$x, datW$block, range))
p <- ggplot(datWSupp) +
aes(x = x, y = y) +
theme_bw() +
theme(
panel.grid.minor = element_blank(),
panel.grid.major = element_blank()
) +
scale_x_continuous(breaks = sort(unique(datW$x))) +
scale_y_continuous(
breaks = sort(unique(datW$y)),
sec.axis = dup_axis()
) +
guides(x = guide_axis(angle = 90)) +
geom_tile(na.rm = TRUE) +
geom_rect(aes(xmin = x - 0.5, xmax = x + 0.5, ymin = y - 0.5, ymax = y + 0.5),
color = "white"
) +
geom_text(
x = sum(xranges[1, ]) / 2, y = 10.7, label = "Block A",
hjust = 0, color = "black"
) +
geom_text(
x = sum(xranges[2, ]) / 2, y = 10.7, label = "Block B",
hjust = 0, color = "black"
) +
geom_vline(
xintercept = max(datW$x[datW$block == "A"]),
color = "black", linetype = "dashed", linewidth = 1
)
p + aes(fill = type) +
labs(title = "Types of microplots") +
scale_fill_discrete()
scaleCols <- c("#CB2027", "#ffec1b", "#b3e93e", "#60BD68", "#059748")
scaleLim <- range(datW$tot_yield)
p + aes(fill = tot_yield) +
labs(title = "Total yield for each microplot") +
scale_fill_continuous(type = "viridis")
Infer the parameters
Using intercrops only
Prepare the inputs
idxIC <- which(!is.na(datW$geno_S1) & !is.na(datW$geno_S2))
datW_IC <- droplevels(datW[idxIC, ])
listY <- list(Y_IC = datW_IC[, c("yield_S1", "yield_S2")])
listX <- list(X_IC = model.matrix(~ 1 + block + geno_S2, datW_IC,
contrasts.arg = list(
"block" = "contr.sum",
"geno_S2" = "contr.sum"
)
))
listZ <- list(Z_DS_f = model.matrix(~ 0 + geno_S1, datW_IC))
colnames(listZ$Z_DS_f) <- gsub("^geno_S1", "", colnames(listZ$Z_DS_f))
listVCov <- list(K = GRMs_vr_trial$S1[
levels(datW_IC$geno_S1),
levels(datW_IC$geno_S1)
])Fit the model
fitsTmb <- list()
i <- 1
for (REML in c(TRUE, FALSE)) {
print(paste0("fit model with ", ifelse(REML, "REML", "ML"), "..."))
st <- system.time(
fitTmb <- fitDBVSBVinter(listY, listX, listZ, listVCov,
lOptions = list(iter.max = 20, quantifUncertain = TRUE),
REML = REML, verbose = 0
)
)
print(st)
fitsTmb[[i]] <- fitTmb
i <- i + 1
break # skip ML to speed-up
}
## [1] "fit model with REML..."
## user system elapsed
## 21.608 0.112 8.463for (i in seq_along(fitsTmb)) {
fitTmb <- fitsTmb[[i]]
p <- ggplot(fitTmb$trace) +
aes(x = iter, y = objfn) +
geom_point() +
geom_line() +
labs(
title = "Optimization convergence",
subtitle = paste0("REML=", fitTmb$REML)
) +
theme_bw()
print(p)
}
Checks
for (i in seq_along(fitsTmb)) {
fitTmb <- fitsTmb[[i]]
print(paste0("REML=", fitTmb$REML))
print("Check fixed effects:")
checks <- data.frame(
param = rownames(fitTmb$B_IC),
truth = c(truth$B_IC),
fitTmb$B_IC,
row.names = NULL,
check.names = FALSE
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(checks)
print("Check (co)variances of random genetic effects:")
checks <- data.frame(
param = rownames(fitTmb$vcov_BV_f),
truth = c(
truth$var_DBV["S1"],
truth$var_SBV["S1"],
truth$cor_DBV_SBV["S1"]
),
fitTmb$vcov_BV_f,
row.names = NULL,
check.names = FALSE
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(checks)
print("Check (co)variances of residual errors:")
checks <- data.frame(
param = rownames(fitTmb$vcov_E_IC),
truth = c(truth$var_err_IC, truth$cor_err_IC),
fitTmb$vcov_E_IC,
row.names = NULL,
check.names = FALSE
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(checks)
print("Check random genetic effects of the focal species:")
checks <- data.frame(
param = c(
rep(c("DBV", "SBV"), each = nrow(truth$BV$S1)),
rep("BV_IC", length(truth$BV_IC$S1))
),
geno = levels(datW_IC$geno_S1),
truth = c(
truth$BV$S1[levels(datW_IC$geno_S1), ],
truth$BV_IC$S1[levels(datW_IC$geno_S1)]
),
rbind(
fitTmb$DBV_f[levels(datW_IC$geno_S1), ],
fitTmb$SBV_IC_f[levels(datW_IC$geno_S1), ],
fitTmb$BV_IC_f[levels(datW_IC$geno_S1), ]
),
row.names = NULL,
check.names = FALSE
)
checks$param <- factor(checks$param,
levels = c("BV_IC", "DBV", "SBV")
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(tapply(1:nrow(checks), checks$param, function(idx) {
cor(checks$truth[idx], checks$Estimate[idx])
}))
p <- ggplot(checks) +
aes(x = Estimate, y = truth) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_abline(slope = 1, intercept = 0, linetype = "dotted") +
geom_point() +
labs(
title = "Results with intercrop-only data",
subtitle = paste0("REML=", fitTmb$REML)
) +
theme_bw() +
facet_wrap(~param)
print(p)
}
## [1] "REML=TRUE"
## [1] "Check fixed effects:"
## param truth Estimate Std. Error z value
## 1 B_IC (Intercept) on species 1 32.0000000 31.4775615 0.22026611 142.906969
## 2 B_IC block1 on species 1 -3.2027775 -3.7124011 0.29227900 -12.701566
## 3 B_IC geno_S21 on species 1 -0.8802921 -0.8068663 0.28970081 -2.785171
## 4 B_IC (Intercept) on species 2 27.0000000 27.2483925 0.09649174 282.390921
## 5 B_IC block1 on species 2 -3.3969120 -3.1818882 0.13005989 -24.464793
## 6 B_IC geno_S21 on species 2 0.6347487 0.6603928 0.12944348 5.101785
## Pr(>|z^2|) nBE
## 1 0.000000e+00 -1.6326203
## 2 5.795807e-37 15.9119270
## 3 5.349953e-03 -8.3410816
## 4 0.000000e+00 0.9199722
## 5 3.502595e-132 -6.3299771
## 6 3.364652e-07 4.0400427
## [1] "Check (co)variances of random genetic effects:"
## param truth Estimate Std. Error z value Pr(>|z^2|) nBE
## 1 var_DBV_f 27.040 20.6557261 4.58011207 4.509874 6.486629e-06 -23.61048
## 2 var_SBV_IC_f 5.408 4.2616651 0.88329605 4.824730 1.401930e-06 -21.19702
## 3 cor_BV_f -0.900 -0.9255841 0.03601392 -25.700733 1.147042e-145 2.84268
## [1] "Check (co)variances of residual errors:"
## param truth Estimate Std. Error z value Pr(>|z^2|)
## 1 vars_E_IC of species 1 4.577666 8.417882 3.4385594 2.448084 0.01436181
## 2 vars_E_IC of species 2 0.975124 1.532500 0.6435353 2.381377 0.01724807
## 3 cor_E_IC -0.200000 -0.518689 0.1860633 -2.787701 0.00530835
## nBE
## 1 83.89028
## 2 57.15948
## 3 159.34448
## [1] "Check random genetic effects of the focal species:"
## BV_IC DBV SBV
## 0.8729402 0.9268083 0.9214229
Using both sole crops and intercrops
Prepare the inputs
idxIC <- which(!is.na(datW$geno_S1) & !is.na(datW$geno_S2))
datW_IC <- droplevels(datW[idxIC, ])
idxSCf <- which(datL$type == "SC" & datL$focal_species == "S1")
datL_SC_f <- droplevels(datL[idxSCf, ])
idxSCt <- which(datL$type == "SC" & datL$focal_species == "S2")
datL_SC_t <- droplevels(datL[idxSCt, ])
listY <- list()
listY$Y_IC <- datW_IC[, c("yield_S1", "yield_S2")]
listY$y_SC_f <- datL_SC_f$focal_yield
listY$y_SC_t <- datL_SC_t$focal_yield
sapply(listY[-1], length)
## y_SC_f y_SC_t
## 300 4
listX <- list()
listX$X_IC <- model.matrix(~ 1 + block + geno_S2,
data = datW_IC,
contrasts.arg = list(
"block" = "contr.sum",
"geno_S2" = "contr.sum"
)
)
listX$X_SC_f <- model.matrix(~ 1 + block, datL_SC_f,
contrasts.arg = list(block = "contr.sum")
)
listX$X_SC_t <- model.matrix(~ 1 + block + focal, datL_SC_t,
contrasts.arg = list(
block = "contr.sum",
focal = "contr.sum"
)
)
listZ <- list()
listZ$Z_DS_f <- model.matrix(~ 0 + geno_S1, datW_IC)
colnames(listZ$Z_DS_f) <- gsub("^geno_S1", "", colnames(listZ$Z_DS_f))
listZ$Z_D_f <- model.matrix(~ 0 + focal, datL_SC_f)
colnames(listZ$Z_D_f) <- gsub("^focal", "", colnames(listZ$Z_D_f))
listVCov <- list(K = GRMs_vr_trial$S1[
levels(datW_IC$geno_S1),
levels(datW_IC$geno_S1)
])Fit the model
fitsTmb <- list()
i <- 1
for (REML in c(TRUE, FALSE)) {
print(paste0("fit model with ", ifelse(REML, "REML", "ML"), "..."))
st <- system.time(
fitTmb <- fitDBVSBVinter(listY, listX, listZ, listVCov,
lOptions = list(iter.max = 20),
REML = REML, verbose = 0
)
)
print(st)
fitsTmb[[i]] <- fitTmb
i <- i + 1
break # skip ML to speed-up
}
## [1] "fit model with REML..."
## user system elapsed
## 38.070 0.145 16.877for (i in seq_along(fitsTmb)) {
fitTmb <- fitsTmb[[i]]
p <- ggplot(fitTmb$trace) +
aes(x = iter, y = objfn) +
geom_point() +
geom_line() +
labs(
title = "Optimization convergence",
subtitle = paste0("REML=", fitTmb$REML)
) +
theme_bw()
print(p)
}
Checks
for (i in seq_along(fitsTmb)) {
fitTmb <- fitsTmb[[i]]
print(paste0("REML=", fitTmb$REML))
print("Check fixed effects:")
checks <- data.frame(
param = rownames(fitTmb$B_IC),
truth = c(truth$B_IC),
fitTmb$B_IC,
row.names = NULL,
check.names = FALSE
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(checks)
checks <- data.frame(
param = rownames(fitTmb$beta_SC_f),
truth = obsMC$blObsContrs$S1[, "SC"],
fitTmb$beta_SC_f,
row.names = NULL,
check.names = FALSE
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(checks)
checks <- data.frame(
param = rownames(fitTmb$beta_SC_t),
truth = c(
obsMC$blObsContrs$S2[, "SC"],
obsMC$BVObsContrs$S2[-1, "SC", "DBV"]
),
fitTmb$beta_SC_t,
row.names = NULL,
check.names = FALSE
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(checks)
print("Check (co)variances of random genetic effects:")
checks <- data.frame(
param = c(
rownames(fitTmb$vcov_BV_f),
rownames(fitTmb$var_SIGV_f)
),
truth = c(
truth$var_DBV["S1"],
truth$var_SBV["S1"],
truth$cor_DBV_SBV["S1"],
truth$var_SIGV["S1"]
),
rbind(
fitTmb$vcov_BV_f,
fitTmb$var_SIGV_f
),
row.names = NULL,
check.names = FALSE
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(checks)
print("Check (co)variances of residual errors:")
checks <- data.frame(
param = c(
rownames(fitTmb$vcov_E_IC),
rownames(fitTmb$var_err_SC_f),
rownames(fitTmb$var_err_SC_t)
),
truth = c(
truth$var_err_IC, truth$cor_err_IC,
truth$var_err_SC
),
rbind(
fitTmb$vcov_E_IC,
fitTmb$var_err_SC_f,
fitTmb$var_err_SC_t
),
row.names = NULL,
check.names = FALSE
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(checks)
print("Check random genetic effects of the focal species:")
checks <- data.frame(
param = rep(c("DBV", "SBV", "BV_IC", "SIGV", "BV_SC"),
each = nrow(truth$BV$S1)
),
geno = levels(datW_IC$geno_S1),
truth = c(
truth$BV$S1[levels(datW_IC$geno_S1), ],
truth$BV_IC$S1[levels(datW_IC$geno_S1)],
truth$SIGVs$S1[levels(datW_IC$geno_S1)],
truth$BV_SC$S1[levels(datW_IC$geno_S1)]
),
rbind(
fitTmb$DBV_f[levels(datW_IC$geno_S1), ],
fitTmb$SBV_IC_f[levels(datW_IC$geno_S1), ],
fitTmb$BV_IC_f[levels(datW_IC$geno_S1), ],
fitTmb$SIGV_f[levels(datW_IC$geno_S1), ],
fitTmb$BV_SC_f[levels(datW_IC$geno_S1), ]
)
)
checks$param <- factor(checks$param,
levels = c("BV_SC", "BV_IC", "DBV", "SBV", "SIGV")
)
checks$nBE <- normBiasError(checks$Estimate, checks$truth)
print(tapply(1:nrow(checks), checks$param, function(idx) {
cor(checks$truth[idx], checks$Estimate[idx])
}))
p <- ggplot(checks) +
aes(x = Estimate, y = truth) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_abline(slope = 1, intercept = 0, linetype = "dotted") +
geom_point() +
labs(
title = "Results with both sole-crop and intercrop data",
subtitle = paste0("REML=", fitTmb$REML)
) +
theme_bw() +
facet_wrap(~param)
print(p)
}
## [1] "REML=TRUE"
## [1] "Check fixed effects:"
## param truth Estimate Std. Error z value
## 1 B_IC (Intercept) on species 1 32.0000000 31.4814126 0.19441322 161.930412
## 2 B_IC block1 on species 1 -3.2027775 -3.5509832 0.23450497 -15.142465
## 3 B_IC geno_S21 on species 1 -0.8802921 -0.8507660 0.23155259 -3.674180
## 4 B_IC (Intercept) on species 2 27.0000000 27.2486929 0.08572208 317.872529
## 5 B_IC block1 on species 2 -3.3969120 -3.2460221 0.10833010 -29.964176
## 6 B_IC geno_S21 on species 2 0.6347487 0.6770458 0.10706495 6.323693
## Pr(>|z^2|) nBE
## 1 0.000000e+00 -1.6205855
## 2 8.496846e-52 10.8719929
## 3 2.386141e-04 -3.3541289
## 4 0.000000e+00 0.9210849
## 5 2.876075e-197 -4.4419708
## 6 2.553851e-10 6.6636111
## param truth Estimate Std. Error z value Pr(>|z^2|)
## 1 beta_SC_f (Intercept) 64.316771 64.366390 0.3311259 194.38646 0.000000e+00
## 2 beta_SC_f block1 -3.765284 -3.637394 0.3338821 -10.89425 1.227779e-27
## nBE
## 1 0.0771481
## 2 -3.3965674
## param truth Estimate Std. Error z value Pr(>|z^2|)
## 1 beta_SC_t (Intercept) 29.7662043 29.7662043 0.3186491 93.413742 0.000000e+00
## 2 beta_SC_t block1 -1.9724155 -1.9724155 0.3186491 -6.189930 6.019102e-10
## 3 beta_SC_t focal1 -0.4266048 -0.4266048 0.3186491 -1.338792 1.806385e-01
## nBE
## 1 0.000000e+00
## 2 9.005997e-14
## 3 0.000000e+00
## [1] "Check (co)variances of random genetic effects:"
## param truth Estimate Std. Error z value Pr(>|z^2|) nBE
## 1 var_DBV_f 27.040 25.8376233 3.12811563 8.259804 1.459114e-16 -4.446660
## 2 var_SBV_IC_f 5.408 5.2285955 0.71172222 7.346399 2.036178e-13 -3.317392
## 3 cor_BV_f -0.900 -0.9246659 0.02298624 -40.226941 0.000000e+00 2.740660
## 4 var_SIGV_f 13.520 6.5281165 5.25134177 1.243133 2.138188e-01 -51.715115
## [1] "Check (co)variances of residual errors:"
## param truth Estimate Std. Error z value
## 1 vars_E_IC of species 1 4.577666 4.31535792 1.7003155 2.5379748
## 2 vars_E_IC of species 2 0.975124 0.78136956 0.3841625 2.0339559
## 3 cor_E_IC -0.200000 -0.07902282 0.3640070 -0.2170915
## 4 var_err_SC_f 17.382857 23.55855699 5.3581985 4.3967309
## 5 var_err_SC_t 2.468571 0.40614898 1.9304652 0.2103892
## Pr(>|z^2|) nBE
## 1 1.114960e-02 -5.730164
## 2 4.195605e-02 -19.869727
## 3 8.281371e-01 -60.488592
## 4 1.098935e-05 35.527530
## 5 8.333639e-01 -83.547206
## [1] "Check random genetic effects of the focal species:"
## BV_SC BV_IC DBV SBV SIGV
## 0.8677352 0.8898539 0.9428374 0.9339428 0.6573680
Reference
See the article for more details:
- Salomon J, Enjalbert J, Flutre T. 2026. Joint modeling of social genetic effects in mono- and pluri-specific groups: case study in intercrops. Genetics iyag177. https://doi.org/10.1093/genetics/iyag177
Appendix
t1 <- proc.time()
t1 - t0
## user system elapsed
## 71.667 0.517 35.513
print(sessionInfo(), locale = FALSE)
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plantmix_1.0.4 lme4_2.0-6 Matrix_1.7-6 ggplot2_4.0.3 knitr_1.52
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 generics_0.1.4 lattice_0.23-1 digest_0.6.39
## [5] magrittr_2.0.5 evaluate_1.0.5 grid_4.6.1 RColorBrewer_1.1-3
## [9] fastmap_1.2.0 jsonlite_2.0.0 viridisLite_0.4.3 scales_1.4.0
## [13] jquerylib_0.1.4 reformulas_0.4.4 Rdpack_2.6.6 cli_3.6.6
## [17] rlang_1.3.0 rbibutils_2.4.1 splines_4.6.1 withr_3.0.3
## [21] cachem_1.1.0 yaml_2.3.12 otel_0.2.0 tools_4.6.1
## [25] nloptr_2.2.1 minqa_1.2.8 dplyr_1.2.1 boot_1.3-32
## [29] buildtools_1.0.0 vctrs_0.7.3 R6_2.6.1 lifecycle_1.0.5
## [33] MASS_7.3-66 pkgconfig_2.0.3 pillar_1.11.1 bslib_0.12.0
## [37] gtable_0.3.6 glue_1.8.1 Rcpp_1.1.2 xfun_0.61
## [41] tibble_3.3.1 tidyselect_1.2.1 sys_3.4.3 farver_2.1.2
## [45] htmltools_0.5.9 nlme_3.1-171 igraph_2.3.4 labeling_0.4.3
## [49] rmarkdown_2.32 maketools_1.3.2 TMB_1.9.25 compiler_4.6.1
## [53] S7_0.2.2