1 Preamble

Dependencies:

suppressPackageStartupMessages(library(plantmix))
suppressPackageStartupMessages(library(ggplot2))

2 Simulate data

set.seed(12345)

2.1 Simulate the genotypes

2.1.1 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"]
  )
)

2.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)]
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        0        0        1        0
table(snpGenos$S1)
## 
##      0      1      2 
## 154021 201021 144958

2.1.3 Additive 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

species <- "S1"
GRM <- GRMs_vr[[species]]
image(Matrix(GRM), main = paste0("GRM for ", species))

hist(diag(GRM), main = paste0("GRM for ", species))

hist(GRM[upper.tri(GRM)], main = paste0("GRM for ", species))

2.2 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())

2.2.1 Design sparsity

plantmix:::plotAllocSchemeInterMixDesign(datW)

2.2.2 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.

2.2.2.1 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.526 0.434 0.521 0.579 0.324 ...
##  $ RY_S2  : num  0.894 0.93 0.838 0.875 1.084 ...
##  $ RYT    : num  1.42 1.36 1.36 1.45 1.41 ...
summary(true_RYTs[, grep("RY_", colnames(true_RYTs))])
##      RY_S1            RY_S2       
##  Min.   :0.2860   Min.   :0.6620  
##  1st Qu.:0.4578   1st Qu.:0.8522  
##  Median :0.4915   Median :0.9081  
##  Mean   :0.4884   Mean   :0.9068  
##  3rd Qu.:0.5200   3rd Qu.:0.9588  
##  Max.   :0.6546   Max.   :1.1258
summary(true_RYTs[, grep("RYT", colnames(true_RYTs))])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.231   1.357   1.392   1.395   1.430   1.531
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.666 0.55 0.66 0.734 0.41 ...
##  $ RYP_S2 : num  4.25 4.42 3.98 4.16 5.15 ...
summary(true_RYPs[, grep("RYP", colnames(true_RYPs))])
##      RYP_S1           RYP_S2     
##  Min.   :0.3622   Min.   :3.145  
##  1st Qu.:0.5799   1st Qu.:4.048  
##  Median :0.6226   Median :4.314  
##  Mean   :0.6186   Mean   :4.307  
##  3rd Qu.:0.6587   3rd Qu.:4.554  
##  Max.   :0.8292   Max.   :5.348

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()

2.2.2.2 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.    NA's 
##  0.8076  0.9229  1.0324  1.0481  1.1635  1.4996     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()

3 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  -2.538 0.893 -9.087 -5.424 3.982 ...
##  $ true_gen_yield_S2: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ true_yield_S1    : num  59.3 62.7 52.7 56.4 65.8 ...
##  $ 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  -2.538 0.893 -9.087 -5.424 3.982 ...
##  $ true_rnd_yield_S2: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ yield_S1         : num  62 72.5 50.4 56.8 66.6 ...
##  $ yield_S2         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ tot_yield        : num  62 72.5 50.4 56.8 66.6 ...
tapply(datW$type, datW$block, table)
## $A
## 
##  SC  IC 
## 152 150 
## 
## $B
## 
##  SC  IC 
## 152 150

ggplot(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()

3.1 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")

4 Infer the parameters

4.1 Using intercrops only

4.1.1 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)
])

4.1.2 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 
##  48.077   0.112  48.193

for (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)
}

4.1.3 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.4684790 0.21016257 149.733987
## 2      B_IC block1 on species 1 -3.2027775 -3.4803658 0.28771959 -12.096381
## 3    B_IC geno_S21 on species 1 -0.8802921 -0.9374902 0.28487500  -3.290882
## 4 B_IC (Intercept) on species 2 27.0000000 27.2755345 0.09771888 279.122472
## 5      B_IC block1 on species 2 -3.3969120 -3.3198444 0.13138060 -25.268909
## 6    B_IC geno_S21 on species 2  0.6347487  0.7107339 0.13015002   5.460882
##      Pr(>|z^2|)       nBE
## 1  0.000000e+00 -1.661003
## 2  1.103704e-33  8.667113
## 3  9.987363e-04  6.497618
## 4  0.000000e+00  1.020498
## 5 7.019581e-141 -2.268753
## 6  4.737745e-08 11.970920
## [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 23.3225601  4.6739948   4.989856  6.042439e-07 -13.747929
## 2 var_SBV_IC_f  5.408  4.6965672  1.0208149   4.600802  4.208673e-06 -13.155192
## 3     cor_BV_f -0.900 -0.9489233  0.0304501 -31.163228 3.357469e-213   5.435919
## [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  6.4159156  3.2937222  1.947923 0.05142423
## 2 vars_E_IC of species 2  0.975124  1.4827962  0.7434509  1.994478 0.04609987
## 3               cor_E_IC -0.200000 -0.3520019  0.2905397 -1.211545 0.22568667
##        nBE
## 1 40.15693
## 2 52.06232
## 3 76.00095
## [1] "Check random genetic effects of the focal species:"
##     BV_IC       DBV       SBV 
## 0.8696013 0.9348076 0.9323475

4.1.4 Parametric bootstrap

if (FALSE) { # slow
  system.time(
    pmTmb <- paramBoot4TMB(fitTmb, nb_boot = 5)
  )
  summary(do.call(rbind, pmTmb))
  fitTmb$sry_sdr[names(pmTmb[[1]]), ]
}

4.2 Using both sole crops and intercrops

4.2.1 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)
])

4.2.2 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 
## 113.088   0.191 113.285

for (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)
}

4.2.3 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.4731476 0.18679895 168.486748
## 2      B_IC block1 on species 1 -3.2027775 -3.3069796 0.22612544 -14.624536
## 3    B_IC geno_S21 on species 1 -0.8802921 -0.8457922 0.22486015  -3.761414
## 4 B_IC (Intercept) on species 2 27.0000000 27.2761557 0.08858871 307.896527
## 5      B_IC block1 on species 2 -3.3969120 -3.3850105 0.11078512 -30.554740
## 6    B_IC geno_S21 on species 2  0.6347487  0.6624723 0.10999834   6.022567
##      Pr(>|z^2|)        nBE
## 1  0.000000e+00 -1.6464137
## 2  1.959026e-48  3.2534930
## 3  1.689553e-04 -3.9191456
## 4  0.000000e+00  1.0227988
## 5 4.891513e-205 -0.3503627
## 6  1.716727e-09  4.3676609
##                   param    truth  Estimate Std. Error    z value   Pr(>|z^2|)
## 1 beta_SC_f (Intercept) 64.21732 64.236791  0.3280420 195.818783 0.000000e+00
## 2      beta_SC_f block1 -3.45987 -3.267909  0.3362193  -9.719575 2.488177e-22
##           nBE
## 1  0.03031617
## 2 -5.54823238
##                   param      truth   Estimate Std. Error   z value   Pr(>|z^2|)
## 1 beta_SC_t (Intercept) 29.7662043 29.7662043  0.3658129 81.370017 0.000000e+00
## 2      beta_SC_t block1 -1.9724155 -1.9724155  0.3658129 -5.391869 6.972846e-08
## 3      beta_SC_t focal1 -0.4266048 -0.4266048  0.3658129 -1.166183 2.435405e-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|)
## 1    var_DBV_f 27.040 27.7665367 3.01063147   9.222828  2.893654e-20
## 2 var_SBV_IC_f  5.408  5.5800529 0.69534080   8.024918  1.015938e-15
## 3     cor_BV_f -0.900 -0.9040776 0.02526856 -35.778749 2.364232e-280
## 4   var_SIGV_f 13.520 10.4906621 4.71838838   2.223357  2.619173e-02
##           nBE
## 1   2.6868959
## 2   3.1814516
## 3   0.4530655
## 4 -22.4063456
## [1] "Check (co)variances of residual errors:"
##                    param     truth    Estimate Std. Error   z value
## 1 vars_E_IC of species 1  4.577666  2.82382527  1.4191541 1.9897947
## 2 vars_E_IC of species 2  0.975124  0.80086910  0.3573848 2.2409154
## 3               cor_E_IC -0.200000  0.01887095  0.3931604 0.0479981
## 4           var_err_SC_f 17.382857 21.08533007  4.6360203 4.5481531
## 5           var_err_SC_t  2.468571  0.53527640  2.9207950 0.1832639
##     Pr(>|z^2|)        nBE
## 1 4.661356e-02  -38.31299
## 2 2.503155e-02  -17.87003
## 3 9.617178e-01 -109.43548
## 4 5.411876e-06   21.29956
## 5 8.545909e-01  -78.31635
## [1] "Check random genetic effects of the focal species:"
##     BV_SC     BV_IC       DBV       SBV      SIGV 
## 0.8902488 0.8960373 0.9460653 0.9390076 0.6416876

5 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

6 Appendix

t1 <- proc.time()
t1 - t0
##    user  system elapsed 
## 182.864   0.428 183.291
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           labeling_0.4.3    
## [17] generics_0.1.4     igraph_2.1.4       rbibutils_2.3      MASS_7.3-65       
## [21] tibble_3.3.0       nloptr_2.2.1       minqa_1.2.8        TMB_1.9.25        
## [25] bslib_0.9.0        pillar_1.11.0      RColorBrewer_1.1-3 rlang_1.1.7       
## [29] cachem_1.1.0       xfun_0.52          sass_0.4.10        viridisLite_0.4.2 
## [33] cli_3.6.6          withr_3.0.2        magrittr_2.0.5     Rdpack_2.6.4      
## [37] digest_0.6.39      grid_4.5.1         lifecycle_1.0.5    nlme_3.1-168      
## [41] reformulas_0.4.4   vctrs_0.6.5        evaluate_1.0.5     glue_1.8.1        
## [45] farver_2.1.2       rmarkdown_2.29     tools_4.5.1        pkgconfig_2.0.3   
## [49] htmltools_0.5.8.1

3) Example of DBV-SBV models for species mixtures in a tester-based design

J. Salomon, J. Enjalbert, T. Flutre

14/10/2026