| Title: | Genetic Study of Plant Mixtures |
|---|---|
| Description: | Fit linear mixed models dedicated to the genetic study of plant mixtures, such as those based on general and specific mixing abilities (GMA-SMA) as well as direct and social breeding values (DBV-SBV), also known as direct and indirect genetic effects (DGE-IGE). More details in Forst et al (2019, <doi:10.1016/j.fcr.2019.107571>) for GMA-SMA models, and Salomon et al (2026, <doi:10.64898/2026.03.27.714849>) for DBV-SBV models. The package also provides functions to optimize experimental designs, simulate data sets and compute interaction indices. |
| Authors: | Timothee Flutre [aut, ctb, cre] (ORCID: <https://orcid.org/0000-0003-4489-4782>), INRAE [cph, fnd], Jerome Enjalbert [ctb] (ORCID: <https://orcid.org/0000-0003-0408-2253>), Emma Forst [ctb] (ORCID: <https://orcid.org/0000-0001-9557-3589>), Maxence Remerand [ctb] (ORCID: <https://orcid.org/0009-0007-3217-5613>), Jemay Salomon [ctb] (ORCID: <https://orcid.org/0000-0003-4454-7916>) |
| Maintainer: | Timothee Flutre <[email protected]> |
| License: | AGPL-3 |
| Version: | 1.0.2 |
| Built: | 2026-06-10 08:21:23 UTC |
| Source: | https://github.com/cran/plantmix |
Computes the change in contribution (CC) proposed by Williams and McCarthy (2001, doi:10.1046/j.1440-1703.2001.00368.x) (this function implements the equation for mixtures with any number of components as written in appendix V of the paper).
CC( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield" )CC( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield" )
dat |
data frame with one yield value per row and with at least four columns: the identifier of each stand, the identifier of the focal genotype or species in each stand, the sowing (plant) proportion of the focal in each stand, and the focal yield in each stand (see the example below); a stand with a single proportion equals to 1 will be interpreted as a monovarietal culture ("pure stand"), a mixture otherwise (and its proportions should sum to 1) |
colIDstand |
column name for the stand identifiers |
colIDfocal |
column name for the focal identifiers |
colProp |
column name for the proportions |
colY |
column name for the yield values |
input data.frame with an extra column named "CC"
Timothee Flutre
(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield=c(50, 40, 25, 22))) CC(dat)(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield=c(50, 40, 25, 22))) CC(dat)
Estimates a genomic relationship matrix with the first estimator from VanRaden (2008).
estimGRM(M, AFs = NULL)estimGRM(M, AFs = NULL)
M |
matrix of SNP genotypes, with individuals in rows and SNPs in columns, encoded in allele dose in \[0,2\] |
AFs |
vector of allele frequencies; if NULL, will be estimated from |
matrix
Timothee Flutre
strong_div_pops <- diag(3) strong_div_pops[upper.tri(strong_div_pops)] <- 0.5 strong_div_pops[lower.tri(strong_div_pops)] <- strong_div_pops[upper.tri(strong_div_pops)] strong_div_pops M <- simulGenosDoseStruct(div_pops=strong_div_pops) A <- estimGRM(M) imageMat(A, "Strong divergence")strong_div_pops <- diag(3) strong_div_pops[upper.tri(strong_div_pops)] <- 0.5 strong_div_pops[lower.tri(strong_div_pops)] <- strong_div_pops[upper.tri(strong_div_pops)] strong_div_pops M <- simulGenosDoseStruct(div_pops=strong_div_pops) A <- estimGRM(M) imageMat(A, "Strong divergence")
Computes the relative yield of mixtures (RYMs) per replicate (e.g., blocks), then fits a linear model correcting for this rep effect, and quantifies the uncertainty around it.
estimRYMRep( data, colY, colID, colRep, mix2genos, conf_level = 0.95, plotDiag = FALSE, title = "" )estimRYMRep( data, colY, colID, colRep, mix2genos, conf_level = 0.95, plotDiag = FALSE, title = "" )
data |
data frame |
colY |
column name specifying the response variable |
colID |
column name specifying the stand identifiers |
colRep |
column name specifying the replicate |
mix2genos |
named list with one component pe rmixture, each component being a vector of genotype names |
conf_level |
confidence level for the uncertainty intervals |
plotDiag |
if TRUE, diagnostics will be plotted |
title |
title for the diagnostics plots |
data frame with the inference summaries of RYM per mixture
Timothee Flutre [aut], Arnaud Gauffreteau [ctb]
## generate fake data set.seed(1234) nbGenos <- 25 genos <- sprintf("g%02i", 1:nbGenos) pairs <- t(combn(x=genos, m=2)) mixIDs <- sample(paste(pairs[,1], pairs[,2], sep="_"), size=75) nbBlocks <- 3 blocks <- LETTERS[1:nbBlocks] dat <- do.call(rbind, lapply(blocks, function(block){ data.frame(ID=c(genos, mixIDs), block=block, stringsAsFactors=TRUE) })) listContr <- list(block="contr.sum") X <- model.matrix(~ 1 + block, data=dat, contrasts=listContr) Z_GMA <- mkZGMA(dat, "ID", "_") Z_SMA <- mkZSMA(dat, "ID", "_", inc_SMA_ii="only_pur") truth <- list("intercept"=100, "var_GMA"=10, "var_SMA"=4, "var_error"=1) truth[["blockEffs"]] <- rnorm(n=nbBlocks - 1, mean=0, sd=2) truth[["GMAs"]] <- rnorm(n=nbGenos, mean=0, sd=sqrt(truth$var_GMA)) truth[["SMAs"]] <- rnorm(n=nlevels(dat$ID), mean=0, sd=sqrt(truth$var_SMA)) truth[["errors"]] <- rnorm(n=nrow(dat), mean=0, sd=sqrt(truth$var_error)) y <- X %*% c(truth$intercept, truth$blockEffs) + Z_GMA %*% truth$GMAs + Z_SMA %*% truth$SMAs + truth$errors dat$yield <- y[,1] hist(dat$yield, breaks=20, las=1, main="Simulated data") boxplot(yield ~ block, data=dat, las=1, main="Simulated data") ## compute the average yield per ID over blocks, and then compute the RYMs: avg <- tapply(dat$yield, dat$ID, mean) avg <- data.frame(ID=names(avg), yield=avg, stringsAsFactors=TRUE) avg <- RYM(avg, colIDcomps="ID", colY="yield", sep="_") ## compute the RYMs per block, and then correct for the "block" effect: mix2genos <- strsplit(levels(dat$ID), "_") names(mix2genos) <- levels(dat$ID) out <- estimRYMRep(dat, "yield", "ID", "block", mix2genos) head(out) ## compare both approaches: op <- par(mfrow=c(1,2)) hist(avg$RYM, las=1, main="Estimated RYM by averaging over blocks", xlab="RYM") abline(v=1, lwd=2); abline(v=mean(avg$RYM, na.rm=TRUE), col="red", lwd=2) hist(out$estim, las=1, main="Estimated RYM after correcting the 'block' effect", xlab="RYM") abline(v=1, lwd=2); abline(v=mean(out$estim), col="red", lwd=2) par(op) out$avgRYM <- avg[rownames(out), "RYM"] plot(out$avgRYM, out$estim, las=1, type="n", xlab="RYM averaged over blocks", ylab="RYM after correcting the 'block' effect", main="Estimated RYMs") idx <- which(out$pv > 0.05) points(out$avgRYM[idx], out$estim[idx], pch=19, col="black") idx <- which(out$pv <= 0.05) points(out$avgRYM[idx], out$estim[idx], pch=19, col="red") abline(a=0, b=1, h=1, v=1, lty=2) segments(x0=as.numeric(out$avgRYM), y0=out$cil, x1=as.numeric(out$avgRYM), y1=out$ciu, col="grey", lty=2, lwd=2)## generate fake data set.seed(1234) nbGenos <- 25 genos <- sprintf("g%02i", 1:nbGenos) pairs <- t(combn(x=genos, m=2)) mixIDs <- sample(paste(pairs[,1], pairs[,2], sep="_"), size=75) nbBlocks <- 3 blocks <- LETTERS[1:nbBlocks] dat <- do.call(rbind, lapply(blocks, function(block){ data.frame(ID=c(genos, mixIDs), block=block, stringsAsFactors=TRUE) })) listContr <- list(block="contr.sum") X <- model.matrix(~ 1 + block, data=dat, contrasts=listContr) Z_GMA <- mkZGMA(dat, "ID", "_") Z_SMA <- mkZSMA(dat, "ID", "_", inc_SMA_ii="only_pur") truth <- list("intercept"=100, "var_GMA"=10, "var_SMA"=4, "var_error"=1) truth[["blockEffs"]] <- rnorm(n=nbBlocks - 1, mean=0, sd=2) truth[["GMAs"]] <- rnorm(n=nbGenos, mean=0, sd=sqrt(truth$var_GMA)) truth[["SMAs"]] <- rnorm(n=nlevels(dat$ID), mean=0, sd=sqrt(truth$var_SMA)) truth[["errors"]] <- rnorm(n=nrow(dat), mean=0, sd=sqrt(truth$var_error)) y <- X %*% c(truth$intercept, truth$blockEffs) + Z_GMA %*% truth$GMAs + Z_SMA %*% truth$SMAs + truth$errors dat$yield <- y[,1] hist(dat$yield, breaks=20, las=1, main="Simulated data") boxplot(yield ~ block, data=dat, las=1, main="Simulated data") ## compute the average yield per ID over blocks, and then compute the RYMs: avg <- tapply(dat$yield, dat$ID, mean) avg <- data.frame(ID=names(avg), yield=avg, stringsAsFactors=TRUE) avg <- RYM(avg, colIDcomps="ID", colY="yield", sep="_") ## compute the RYMs per block, and then correct for the "block" effect: mix2genos <- strsplit(levels(dat$ID), "_") names(mix2genos) <- levels(dat$ID) out <- estimRYMRep(dat, "yield", "ID", "block", mix2genos) head(out) ## compare both approaches: op <- par(mfrow=c(1,2)) hist(avg$RYM, las=1, main="Estimated RYM by averaging over blocks", xlab="RYM") abline(v=1, lwd=2); abline(v=mean(avg$RYM, na.rm=TRUE), col="red", lwd=2) hist(out$estim, las=1, main="Estimated RYM after correcting the 'block' effect", xlab="RYM") abline(v=1, lwd=2); abline(v=mean(out$estim), col="red", lwd=2) par(op) out$avgRYM <- avg[rownames(out), "RYM"] plot(out$avgRYM, out$estim, las=1, type="n", xlab="RYM averaged over blocks", ylab="RYM after correcting the 'block' effect", main="Estimated RYMs") idx <- which(out$pv > 0.05) points(out$avgRYM[idx], out$estim[idx], pch=19, col="black") idx <- which(out$pv <= 0.05) points(out$avgRYM[idx], out$estim[idx], pch=19, col="red") abline(a=0, b=1, h=1, v=1, lty=2) segments(x0=as.numeric(out$avgRYM), y0=out$cil, x1=as.numeric(out$avgRYM), y1=out$ciu, col="grey", lty=2, lwd=2)
Computes the relative yields (RYs) per replicate (e.g., blocks), then fits a linear model correcting for this rep effect, and quantifies the uncertainty around it.
estimRYRep( data, colY, colIDstand, colIDfocal, colRep, mix2genos, conf_level = 0.95, plotDiag = FALSE )estimRYRep( data, colY, colIDstand, colIDfocal, colRep, mix2genos, conf_level = 0.95, plotDiag = FALSE )
data |
data frame |
colY |
column name specifying the response variable |
colIDstand |
column name specifying the stand identifiers |
colIDfocal |
column name for the focal identifiers |
colRep |
column name specifying the replicate |
mix2genos |
named list with one component per mixture, each component being a vector of genotype names; see |
conf_level |
confidence level for the uncertainty intervals |
plotDiag |
if TRUE, diagnostic plots |
data frame with the inference summaries of RY per mixture
Timothee Flutre [aut], Arnaud Gauffreteau [ctb]
(dat <- data.frame(ID=c("geno1", "geno1", "geno2", "geno2", "mixg1g2", "mixg1g2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno1", "geno2", "geno2", "geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5), block=c("A","B","A","B","A","A","B","B"), yield=c(51,45, 39,43, 25,22,26,25))) estimRYRep(dat, "yield", "ID", "focal", "block", list("mixg1g2"=c("geno1","geno2")))(dat <- data.frame(ID=c("geno1", "geno1", "geno2", "geno2", "mixg1g2", "mixg1g2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno1", "geno2", "geno2", "geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 1, 1, 0.5, 0.5, 0.5, 0.5), block=c("A","B","A","B","A","A","B","B"), yield=c(51,45, 39,43, 25,22,26,25))) estimRYRep(dat, "yield", "ID", "focal", "block", list("mixg1g2"=c("geno1","geno2")))
Fits DBV-SBV models for interspecific mixtures.
fitDBVSBVinter( listY, listX, listZ, listVCov, REML = TRUE, lOptions = NULL, verbose = FALSE )fitDBVSBVinter( listY, listX, listZ, listVCov, REML = TRUE, lOptions = NULL, verbose = FALSE )
listY |
named list of response variables (a matrix Y_IC and, optionally, a vector y_SC_f and a vector y_SC_t) |
listX |
named list of design matrices for the fixed-effect explanatory factors (a matrix X_IC and, optionally, a matrix X_SC_f and a matrix X_SC_t) |
listZ |
named list of design matrices of random-effect explanatory factors |
listVCov |
named list of square, symmetric matrices used, after re-scaling, as variance-covariance matrices for the random-effect factors; named "K" for the BVs (DBVs and SBVs) and "Kmixpair" for the DBVxSBVs; these matrices must have dimnames (both rows and columns) coherent with the column names of the design matrices in |
REML |
logical |
lOptions |
named list of options (for experts) |
verbose |
verbosity level |
list
Jemay Salomon, Timothee Flutre
## simulate a data set with both sole crops and intercrops: GRMs <- list("S1" = diag(100), "S2" = diag(2)) dimnames(GRMs$S1) <- list(paste0("gS1_", 1:100), paste0("gS1_", 1:100)) dimnames(GRMs$S2) <- list(paste0("gS2_", 1:2), paste0("gS2_", 1:2)) out <- simulDBVSBVinter(GRMs) names(out) datW <- out$datW str(datW) ## fit the model using intercrops only: 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$S1[levels(datW_IC$geno_S1), levels(datW_IC$geno_S1)]) fitTmb <- fitDBVSBVinter(listY, listX, listZ, listVCov, lOptions = list(iter.max = 20), REML = TRUE) names(fitTmb) fitTmb$sdr ## see the third vignette for more details## simulate a data set with both sole crops and intercrops: GRMs <- list("S1" = diag(100), "S2" = diag(2)) dimnames(GRMs$S1) <- list(paste0("gS1_", 1:100), paste0("gS1_", 1:100)) dimnames(GRMs$S2) <- list(paste0("gS2_", 1:2), paste0("gS2_", 1:2)) out <- simulDBVSBVinter(GRMs) names(out) datW <- out$datW str(datW) ## fit the model using intercrops only: 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$S1[levels(datW_IC$geno_S1), levels(datW_IC$geno_S1)]) fitTmb <- fitDBVSBVinter(listY, listX, listZ, listVCov, lOptions = list(iter.max = 20), REML = TRUE) names(fitTmb) fitTmb$sdr ## see the third vignette for more details
Fits GMA-SMA models.
fitGMASMA( formFix, data, listZ, pkg = "lme4", listVCov, contrasts = NULL, REML = TRUE, ... )fitGMASMA( formFix, data, listZ, pkg = "lme4", listVCov, contrasts = NULL, REML = TRUE, ... )
formFix |
formula whose right-hand side should only include the fixed effects, the random effects being deduced based on the names of |
data |
data frame; missing data will be removed |
listZ |
named list of design matrices for the random effects; its names should be |
pkg |
package used to fit the model among lme4, MM4LMM, TMB, INLA and breedR; lme4 and INLA do not use |
listVCov |
named list of variance-covariance matrices for the random effects; the names of this list should be the same as the names of |
contrasts |
|
REML |
logical |
... |
additional arguments specific to each package |
depends on the chosen package
Timothee Flutre
## see the first and second vignettes## see the first and second vignettes
Optimizes a design for binary crop mixtures
getDesignBinaryCropMix( levGenos1, levGenos2, nbMixes, seed = NULL, showPlots = FALSE, verbose = FALSE )getDesignBinaryCropMix( levGenos1, levGenos2, nbMixes, seed = NULL, showPlots = FALSE, verbose = FALSE )
levGenos1 |
character vector with the names of the genotypes of the first species |
levGenos2 |
character vector with the names of the genotypes of the second species |
nbMixes |
number of mixtures |
seed |
seed for the pseudo-random number generator |
showPlots |
logical indicating if plots should be showed |
verbose |
verbosity level |
list
Timothee Flutre
nbGenos1 <- 30 levGenos1 <- sprintf(fmt=paste0("S1_%0", floor(log10(nbGenos1))+1, "i"), 1:nbGenos1) nbGenos2 <- 8 levGenos2 <- sprintf(fmt=paste0("S2_%0", floor(log10(nbGenos2))+1, "i"), 1:nbGenos2) nbMixes <- 60 # only binary and balanced design <- getDesignBinaryCropMix(levGenos1, levGenos2, nbMixes, seed=12345) plotDesignCropMix(design$graph, levGenos1, levGenos2) ggplot(design$entropies) + aes(x=iteration, y=entropy, group=type) + geom_point() + geom_line() + geom_hline(yintercept = design$max_ent) + theme_bw() + facet_wrap(~ type)nbGenos1 <- 30 levGenos1 <- sprintf(fmt=paste0("S1_%0", floor(log10(nbGenos1))+1, "i"), 1:nbGenos1) nbGenos2 <- 8 levGenos2 <- sprintf(fmt=paste0("S2_%0", floor(log10(nbGenos2))+1, "i"), 1:nbGenos2) nbMixes <- 60 # only binary and balanced design <- getDesignBinaryCropMix(levGenos1, levGenos2, nbMixes, seed=12345) plotDesignCropMix(design$graph, levGenos1, levGenos2) ggplot(design$entropies) + aes(x=iteration, y=entropy, group=type) + geom_point() + geom_line() + geom_hline(yintercept = design$max_ent) + theme_bw() + facet_wrap(~ type)
Optimizes a design for binary varietal mixtures
getDesignBinaryVarMix( levGenos, nbMixes, seed = NULL, showPlots = FALSE, verbose = FALSE )getDesignBinaryVarMix( levGenos, nbMixes, seed = NULL, showPlots = FALSE, verbose = FALSE )
levGenos |
character vector with the names of the genotypes |
nbMixes |
number of mixtures |
seed |
seed for the pseudo-random number generator |
showPlots |
logical indicating if plots should be showed |
verbose |
verbosity level |
list
Timothee Flutre
nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) plotDesignVarMix(design$graph, levGenos) ggplot(design$entropies) + aes(x=iteration, y=entropy) + geom_point() + geom_line() + geom_hline(yintercept = design$max_ent) + theme_bw()nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) plotDesignVarMix(design$graph, levGenos) ggplot(design$entropies) + aes(x=iteration, y=entropy) + geom_point() + geom_line() + geom_hline(yintercept = design$max_ent) + theme_bw()
Reformats composition of mixtures into a list.
getMixtureList(mixtures, sep = NULL)getMixtureList(mixtures, sep = NULL)
mixtures |
vector, matrix, data.frame or list containing of mixtures |
sep |
required only if |
list of vectors, with one component per mixture, the vector containing the mixture components
Timothee Flutre
nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) mixtures <- getMixtureList(design$combs) str(mixtures, list.len=10)nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) mixtures <- getMixtureList(design$combs) str(mixtures, list.len=10)
Returns the list of mixtures per genotype.
getMixturesPerGeno(mix2genos)getMixturesPerGeno(mix2genos)
mix2genos |
list with one component per mixture, listing the genotypes it contains (output of |
list with one vector per genotype
Timothee Flutre
nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) mixtures <- getMixtureList(design$combs) str(mixtures, list.len=10) geno2mixes <- getMixturesPerGeno(mixtures) geno2mixes[c(1,2)] table(sapply(geno2mixes, length))nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) mixtures <- getMixtureList(design$combs) str(mixtures, list.len=10) geno2mixes <- getMixturesPerGeno(mixtures) geno2mixes[c(1,2)] table(sapply(geno2mixes, length))
Plots an image of a matrix.
imageMat(mat, title, col, breaks)imageMat(mat, title, col, breaks)
mat |
matrix |
title |
optional title of the plot; if missing, it will be the matrix name and dimensions |
col |
optional vector of colors; if missing, it will be a grey color gradient |
breaks |
optional vector of breaks to go with col |
nothing
Timothee Flutre
M <- diag(c(rep(2, 3), rep(15, 5))) imageMat(M) # see how the the "2" values are ignored! imageMat(M, col=c("lightgrey","red","blue"), breaks=c(-1,1,3,20))M <- diag(c(rep(2, 3), rep(15, 5))) imageMat(M) # see how the the "2" values are ignored! imageMat(M, col=c("lightgrey","red","blue"), breaks=c(-1,1,3,20))
Returns an information criterion, among AIC, AICc and BIC.
Caution: choosing k and n may not be straightforward for certain models, such as mixed models.
infoCriterion(k, lnLmax, n = NULL, type = "AIC")infoCriterion(k, lnLmax, n = NULL, type = "AIC")
k |
number of parameters, sometimes also called "effective degrees of freedom" |
lnLmax |
value of the log-likelihood when maximized |
n |
sample size |
type |
specific information criterion to be returned |
numeric with attributes k and n (if not NULL)
Timothee Flutre
Applies the inverse vec operator to a vector.
invvec(x, n_col, n_row = NULL, sep = NULL)invvec(x, n_col, n_row = NULL, sep = NULL)
x |
vector |
n_col |
number of columns of the output matrix |
n_row |
number of rows of the output matrix |
sep |
separator used to retrieve row and column names, only if |
matrix
Timothee Flutre
mat <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE, dimnames = list(as.character(1:2), letters[1:3])) mat (x <- vec(mat)) invvec(x, n_col = 3, sep = "-")mat <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE, dimnames = list(as.character(1:2), letters[1:3])) mat (x <- vec(mat)) invvec(x, n_col = 3, sep = "-")
Tranforms a vector of component observations into a matrix with one row per mixture and components in colums. Requires all mixtures to have the same number of components, and all the obervationss from the same mixture to be one after the others, in the same order for each mixture.
invvecMixes(x, nb_comps)invvecMixes(x, nb_comps)
x |
vector; make sure that it is correctly sorted beforehand |
nb_comps |
number of components in each mixture |
matrix with nb_comps columns
Timothee Flutre
datL <- data.frame(ID = c("g1+g2", "g1+g2", "g1+g2", "g1+g2"), focal = c("g1", "g2", "g1", "g2"), neighbor = c("g2", "g1", "g2", "g1"), block = c("A", "A", "B", "B"), yield = c(1, 2, 3, 4) ) invvecMixes(datL$yield, 2)datL <- data.frame(ID = c("g1+g2", "g1+g2", "g1+g2", "g1+g2"), focal = c("g1", "g2", "g1", "g2"), neighbor = c("g2", "g1", "g2", "g1"), block = c("A", "A", "B", "B"), yield = c(1, 2, 3, 4) ) invvecMixes(datL$yield, 2)
Computes the land equivalent ratio (LER) of Willey and Osiru (1972, doi:10.1017/S0021859600025909) interpreted as the relative land area required under sole cropping to produce the same yield as under intercropping. It corresponds to equation 28 in Weigelt and Jolliffe (2003, doi:10.1046/j.1365-2745.2003.00805.x), and is equivalent to the index called "relative yield total" (RYT) corresponding to equation 26 of the same paper. Missing values are forbidden.
LER(x)LER(x)
x |
data frame with species in rows (with species names as row names) and at least two columns, the first being the performance under sole-cropping and the second under inter-cropping |
list with partial and total LERs
Timothee Flutre
(dat <- data.frame(solecrop=c(5, 15), intercrop=c(4, 9), row.names=c("grain", "fruit"))) LER(dat)(dat <- data.frame(solecrop=c(5, 15), intercrop=c(4, 9), row.names=c("grain", "fruit"))) LER(dat)
Computes the seed weight per genotype for each stand, assuming equal sowing proportions.
mixSowingWeight(pureSowingWeights, stands)mixSowingWeight(pureSowingWeights, stands)
pureSowingWeights |
vector which names are genotypes and values are sowing weights in pure stands |
stands |
vector which names are stand identifiers and values are genotype(s), a single one if it is a pure stand, several ones if it is a mixed stand (separated by a dash "-", e.g., "var1-var8") |
list which components are stands and values are named vectors with sowing weight per genotype for each stand
Timothee Flutre
sowingArea <- 9.52 sowingDensity <- 160 (tmp <- nbSeedsToSownInPure(sowingArea, sowingDensity)) TKWs <- c("var1"=38.605, "var2"=40.051, "var6"=36.251, "var8"=33.368) pureSowingWeights <- TKWs * tmp stands <- c("var1"="var1", "mix8"="var1-var2", "mix34"="var1-var8", "mix50"="var1-var6-var8") mixSowingWeight(pureSowingWeights, stands)sowingArea <- 9.52 sowingDensity <- 160 (tmp <- nbSeedsToSownInPure(sowingArea, sowingDensity)) TKWs <- c("var1"=38.605, "var2"=40.051, "var6"=36.251, "var8"=33.368) pureSowingWeights <- TKWs * tmp stands <- c("var1"="var1", "mix8"="var1-var2", "mix34"="var1-var8", "mix50"="var1-var6-var8") mixSowingWeight(pureSowingWeights, stands)
Makes the design matrices with the contribution of specific mixing abilities (SMAs).
More details in the help of mkZSMA and in the vignette.
mkAllZSMA( df, col, sep, models = c("2p", "2", "2pp", "3", "3p"), verbose = FALSE )mkAllZSMA( df, col, sep, models = c("2p", "2", "2pp", "3", "3p"), verbose = FALSE )
df |
data frame |
col |
name of the column containing the names of pure and mixed stands |
sep |
separator used to separate the genotype names in |
models |
vector of model identifiers; see the vignette |
verbose |
verbosity level |
list of SMA design matrices
Timothee Flutre
dat <- data.frame(mix=paste0("mix", 1:3), varieties=c("var3","var1-var3","var1-var2"), pheno=c(10, 7, 9)) (listZSMA <- mkAllZSMA(df=dat, col="varieties", sep="-"))dat <- data.frame(mix=paste0("mix", 1:3), varieties=c("var3","var1-var3","var1-var2"), pheno=c(10, 7, 9)) (listZSMA <- mkAllZSMA(df=dat, col="varieties", sep="-"))
Makes a design matrix with the contribution of general mixing abilities (GMAs). The design is based on the names of pure and mixed stands, in which the genotype names are separated by a specific symbol.
mkZGMA(df, col, sep = ",")mkZGMA(df, col, sep = ",")
df |
data frame |
col |
name of the column containing the names of pure and mixed stands |
sep |
separator used to separate the genotype names in |
matrix
Timothee Flutre
dat <- data.frame(mix=paste0("mix", 1:3), varieties=c("var3","var1-var3","var1-var2"), pheno=c(10, 7, 9)) (Z_GMA <- mkZGMA(df=dat, col="varieties", sep="-"))dat <- data.frame(mix=paste0("mix", 1:3), varieties=c("var3","var1-var3","var1-var2"), pheno=c(10, 7, 9)) (Z_GMA <- mkZGMA(df=dat, col="varieties", sep="-"))
Makes the design matrices for direct and social breeding values per species when the trial includes binary intercrops and, possibly, sole crops.
mkZinterspe(df, genosPerSp, colIDfocal = "focal", colIDneighbors = "neighbors")mkZinterspe(df, genosPerSp, colIDfocal = "focal", colIDneighbors = "neighbors")
df |
data frame with one row per sole-crop plot or two rows per binary intercrop plot, with a factor column corresponding to the focal identifiers ("DBV" a.k.a. "DGE" or "Pr") and a column corresponding to the "neighbor(s)" identifiers ("SBV" a.k.a. "IGE" or "As") |
genosPerSp |
list of two components, one per species, each being a vector with the genotype identifiers of the given species (they will be sorted) |
colIDfocal |
name of the column containing the identifiers of the focal genotypes |
colIDneighbors |
name of the column containing the identifiers of the neighbour genotypes |
list of two components, one per species, each being a list of two matrices, the first for the DBVs and the second for the SBVs; the column names of these matrices correspond to the genotype identifiers after sorting
Timothee Flutre
(dat <- data.frame(focal=c("wheat01", "pea02", "wheat01", "pea01", "wheat01"), neighbors=c("pea02", "wheat01", "pea01", "wheat01", "wheat01"), name=c("wheat01-pea02", "wheat01-pea02", "wheat01-pea01", "wheat01-pea01", "wheat01-wheat01"), stand = c("inter", "inter", "inter", "inter", "sole"), species = c("wheat", "pea", "wheat", "pea", "wheat"))) (out <- mkZinterspe(dat, list("wheat"="wheat01", "pea"=c("pea01","pea02"))))(dat <- data.frame(focal=c("wheat01", "pea02", "wheat01", "pea01", "wheat01"), neighbors=c("pea02", "wheat01", "pea01", "wheat01", "wheat01"), name=c("wheat01-pea02", "wheat01-pea02", "wheat01-pea01", "wheat01-pea01", "wheat01-wheat01"), stand = c("inter", "inter", "inter", "inter", "sole"), species = c("wheat", "pea", "wheat", "pea", "wheat"))) (out <- mkZinterspe(dat, list("wheat"="wheat01", "pea"=c("pea01","pea02"))))
Makes a design matrix with the contribution of specific mixing abilities (SMAs). Before Forst et al (2019), there was only one kind of SMA, the inter-genotypic SMA (SMA_ij). Forst et al (2019) introduced a second kind of SMA, the intra-genotypic SMA (SMA_ii). As a result, various design matrices can be made depending on which kinds of SMAs are modeled. The design is based on the names of pure and mixed stands, in which the genotype names are separated by a specific symbol.
mkZSMA(df, col, sep = ",", inc_SMA_ii = "no", skipUnusedCols = TRUE)mkZSMA(df, col, sep = ",", inc_SMA_ii = "no", skipUnusedCols = TRUE)
df |
data frame |
col |
name of the column containing the names of pure and mixed stands |
sep |
separator used to separate the genotype names in |
inc_SMA_ii |
specify how the intra-genotypic SMA should be included, with "no" meaning no SMA_ii, "only_pur" meaning that the SMA_ii is only included in the pure stands (this is model 2 of Forst et al, 2019), "pur_mix" meaning that the SMA_ii is included in both the pure and the mixed stands (this is model 3 of Forst et al, 2019) |
skipUnusedCols |
skip unused columns, i.e., full of zeroes and not present in the data frame |
matrix
Timothee Flutre
dat <- data.frame(mix=paste0("mix", 1:3), varieties=c("var3","var1-var3","var1-var2"), pheno=c(10, 7, 9)) (Z_SMA <- mkZSMA(df=dat, col="varieties", sep="-", inc_SMA_ii="only_pur")) # model 2 (Z_SMA <- mkZSMA(df=dat, col="varieties", sep="-", inc_SMA_ii="pur_mix")) # model 3dat <- data.frame(mix=paste0("mix", 1:3), varieties=c("var3","var1-var3","var1-var2"), pheno=c(10, 7, 9)) (Z_SMA <- mkZSMA(df=dat, col="varieties", sep="-", inc_SMA_ii="only_pur")) # model 2 (Z_SMA <- mkZSMA(df=dat, col="varieties", sep="-", inc_SMA_ii="pur_mix")) # model 3
Computes the number of seeds of a given genotype to be sowned as pure stand based on a given sowing area and density.
nbSeedsToSownInPure(sowingArea = 8.4, sowingDensity = 160)nbSeedsToSownInPure(sowingArea = 8.4, sowingDensity = 160)
sowingArea |
area to be sowned in square meters |
sowingDensity |
number of seeds per square meter |
number of seeds (in thousands)
Timothee Flutre
sowingArea <- 9.52 sowingDensity <- 160 nbSeedsToSownInPure(sowingArea, sowingDensity) sowingArea <- 8.66 sowingDensity <- 200 nbSeedsToSownInPure(sowingArea, sowingDensity)sowingArea <- 9.52 sowingDensity <- 160 nbSeedsToSownInPure(sowingArea, sowingDensity) sowingArea <- 8.66 sowingDensity <- 200 nbSeedsToSownInPure(sowingArea, sowingDensity)
Returns the normalized bias error (NBE).
normBiasError(estim, truth, perc = TRUE)normBiasError(estim, truth, perc = TRUE)
estim |
numeric vector of estimates |
truth |
numeric vector of true values |
perc |
logical; if TRUE, the return value will be in percentage |
numeric vector
Timothee Flutre
Performs parametric bootstrap with TMB and nlminb on a fitted model to allow uncertainty quantification.
paramBoot4TMB(fit, nb_boot = 5, mc.cores = 1)paramBoot4TMB(fit, nb_boot = 5, mc.cores = 1)
fit |
list corresponding to a model fitted with TMB |
nb_boot |
number of parametric bootstraps to perform |
mc.cores |
the number of cores to use, i.e. at most how many child processes will be run simultaneously (see the "parallel" package) |
list with as many components as nb_boot
Jemay Salomon, Timothee Flutre
## see the example of `fitDBVSBVinter` ## then run (slow if nb_boot is high!): paramBoot4TMB(fitTmb)## see the example of `fitDBVSBVinter` ## then run (slow if nb_boot is high!): paramBoot4TMB(fitTmb)
Pivots mixture data into the long format, from one row per stand into one row per component of each stand (i.e., several rows if the stand is a mixture).
pivotMixData2Long( df, genos, colC, sep = "-", prefixY = NULL, sepY = "_", colIDfocal = "focal", colIDneighbors = "neighbor" )pivotMixData2Long( df, genos, colC, sep = "-", prefixY = NULL, sepY = "_", colIDfocal = "focal", colIDneighbors = "neighbor" )
df |
data frame |
genos |
named list of vectors, one per species, containing all the identifiers of the genotypes of the given species |
colC |
name of the column containing the component(s) of each stand |
sep |
separator used to separate the genotype names in |
prefixY |
prefix of the columns containing the yield data (one column per mixture component); it is assumed that the components in |
sepY |
separator used in the name of the new "yield" columns made by concatenating |
colIDfocal |
name of the column in the output that will contain the identifiers of the focal genotypes |
colIDneighbors |
name of the column in the output that will contain the identifiers of the neighbor genotypes; for monovarietal stands, the entries in this neighbor column will be identical to the entries in the focal column |
data.frame with more rows
## example of species mixtures: (dat <- data.frame(name=c("wheat01-pea02","wheat01-pea01", "wheat01"), stand=c("inter","inter","sole"), x=c(1, 1, 1), y=c(1, 2, 3), "yield_cereal" = c(10, 11, 20), "yield_legume" = c(9, 8, NA), check.names = FALSE, stringsAsFactors=TRUE)) pivotMixData2Long(df=dat, colC="name", sep="-", prefixY="yield", genos=list("cereal"=c("wheat01"), "legume"=c("pea01","pea02"))) ## example of varietal mixtures: (dat <- data.frame(name=c("g1+g2","g1+g2","g1","g2"), "yield_focal" = c(10, 12, 20, 18), "yield_neighbor" = c(11, 9, NA, NA), check.names = FALSE, stringsAsFactors=TRUE)) pivotMixData2Long(df=dat, colC="name", sep="+", prefixY="yield", genos=list(c("g1","g2")))## example of species mixtures: (dat <- data.frame(name=c("wheat01-pea02","wheat01-pea01", "wheat01"), stand=c("inter","inter","sole"), x=c(1, 1, 1), y=c(1, 2, 3), "yield_cereal" = c(10, 11, 20), "yield_legume" = c(9, 8, NA), check.names = FALSE, stringsAsFactors=TRUE)) pivotMixData2Long(df=dat, colC="name", sep="-", prefixY="yield", genos=list("cereal"=c("wheat01"), "legume"=c("pea01","pea02"))) ## example of varietal mixtures: (dat <- data.frame(name=c("g1+g2","g1+g2","g1","g2"), "yield_focal" = c(10, 12, 20, 18), "yield_neighbor" = c(11, 9, NA, NA), check.names = FALSE, stringsAsFactors=TRUE)) pivotMixData2Long(df=dat, colC="name", sep="+", prefixY="yield", genos=list(c("g1","g2")))
Pivots mixture data into the wide format, from one row per component of each stand (i.e., several rows if the stand is a mixture) into one row per stand.
pivotMixData2Wide( df, colIDstand = "ID", colIDfocal = "focal", colIDneighbors = "neighbor", colPlot = c("x", "y"), colY = "yield", sepY = "_", sepFocalNeighbors = NULL )pivotMixData2Wide( df, colIDstand = "ID", colIDfocal = "focal", colIDneighbors = "neighbor", colPlot = c("x", "y"), colY = "yield", sepY = "_", sepFocalNeighbors = NULL )
df |
data frame |
colIDstand |
name of the column containing the identifier of each stand |
colIDfocal |
name of the column containing the identifiers of the focal genotypes |
colIDneighbors |
optional name of the column containing the identifiers of the neighbor genotypes |
colPlot |
name of the column(s) allowing to uniquely identify a plot |
colY |
name of the column containing the yield data |
sepY |
separator used in the name of the new "yield" columns made by concatenating |
sepFocalNeighbors |
in case |
data.frame
Timothee Flutre
## only binary mixtures: dat0 <- data.frame( focal = c("wheat01", "pea02", "wheat01", "pea01"), neighbor = c("pea02", "wheat01", "pea01", "wheat01"), name = c( rep("wheat01-pea02", 2), rep("wheat01-pea01", 2) ), stand = rep("inter", 4), x = 1, y = c(1, 1, 2, 2), yield = c(10, 11, 12, 13), stringsAsFactors = TRUE ) dat0 (dat1 <- pivotMixData2Wide(dat0, colIDstand="name")) ## binary mixtures and a monovarietal: dat0 <- data.frame( focal = c("wheat01", "pea02", "wheat01", "pea01", "wheat01"), neighbor = c("pea02", "wheat01", "pea01", "wheat01", "wheat01"), name = c( rep("wheat01-pea02", 2), rep("wheat01-pea01", 2), "wheat01" ), stand = c(rep("inter", 4), "sole"), x = 1, y = c(1, 1, 2, 2, 3), yield = c(10, 11, 12, 13, 20.5), stringsAsFactors = TRUE ) dat0 (dat1 <- pivotMixData2Wide(dat0, colIDstand="name")) ## reverse conversion dat1v2 <- dat1 colnames(dat1v2)[colnames(dat1v2) == "yield_focal"] <- "yield-cereal" colnames(dat1v2)[colnames(dat1v2) == "yield_neighbor"] <- "yield-legume" dat1v2$yield <- apply(dat1v2[,c(7,8)], 1, sum, na.rm=TRUE) (dat1v2 <- dat1v2[,-c(1,2)]) (dat2 <- pivotMixData2Long(dat1v2, colC="name", sepY="-", genos=list("cereal"=c("wheat01","wheat02"), "legume"=c("pea01","pea02")))) all.equal(dat2, dat0)## only binary mixtures: dat0 <- data.frame( focal = c("wheat01", "pea02", "wheat01", "pea01"), neighbor = c("pea02", "wheat01", "pea01", "wheat01"), name = c( rep("wheat01-pea02", 2), rep("wheat01-pea01", 2) ), stand = rep("inter", 4), x = 1, y = c(1, 1, 2, 2), yield = c(10, 11, 12, 13), stringsAsFactors = TRUE ) dat0 (dat1 <- pivotMixData2Wide(dat0, colIDstand="name")) ## binary mixtures and a monovarietal: dat0 <- data.frame( focal = c("wheat01", "pea02", "wheat01", "pea01", "wheat01"), neighbor = c("pea02", "wheat01", "pea01", "wheat01", "wheat01"), name = c( rep("wheat01-pea02", 2), rep("wheat01-pea01", 2), "wheat01" ), stand = c(rep("inter", 4), "sole"), x = 1, y = c(1, 1, 2, 2, 3), yield = c(10, 11, 12, 13, 20.5), stringsAsFactors = TRUE ) dat0 (dat1 <- pivotMixData2Wide(dat0, colIDstand="name")) ## reverse conversion dat1v2 <- dat1 colnames(dat1v2)[colnames(dat1v2) == "yield_focal"] <- "yield-cereal" colnames(dat1v2)[colnames(dat1v2) == "yield_neighbor"] <- "yield-legume" dat1v2$yield <- apply(dat1v2[,c(7,8)], 1, sum, na.rm=TRUE) (dat1v2 <- dat1v2[,-c(1,2)]) (dat2 <- pivotMixData2Long(dat1v2, colC="name", sepY="-", genos=list("cereal"=c("wheat01","wheat02"), "legume"=c("pea01","pea02")))) all.equal(dat2, dat0)
Plots design for crop mixtures, as a graph and as a diallel.
plotDesignCropMix(graph, levGenos1, levGenos2, main = NULL)plotDesignCropMix(graph, levGenos1, levGenos2, main = NULL)
graph |
graph |
levGenos1 |
character vector with the names of the genotypes of the first species |
levGenos2 |
character vector with the names of the genotypes of the second species |
main |
main title for the graph plot |
nothing
Timothee Flutre
nbGenos1 <- 30 levGenos1 <- sprintf(fmt=paste0("S1_%0", floor(log10(nbGenos1))+1, "i"), 1:nbGenos1) nbGenos2 <- 8 levGenos2 <- sprintf(fmt=paste0("S2_%0", floor(log10(nbGenos2))+1, "i"), 1:nbGenos2) nbMixes <- 60 # only binary and balanced design <- getDesignBinaryCropMix(levGenos1, levGenos2, nbMixes, seed=12345) plotDesignCropMix(design$graph, levGenos1, levGenos2)nbGenos1 <- 30 levGenos1 <- sprintf(fmt=paste0("S1_%0", floor(log10(nbGenos1))+1, "i"), 1:nbGenos1) nbGenos2 <- 8 levGenos2 <- sprintf(fmt=paste0("S2_%0", floor(log10(nbGenos2))+1, "i"), 1:nbGenos2) nbMixes <- 60 # only binary and balanced design <- getDesignBinaryCropMix(levGenos1, levGenos2, nbMixes, seed=12345) plotDesignCropMix(design$graph, levGenos1, levGenos2)
Plots design for varietal mixtures, as a graph and as a diallel.
plotDesignVarMix( graph, levGenos, main = NULL, subplots = c("graph", "diallel") )plotDesignVarMix( graph, levGenos, main = NULL, subplots = c("graph", "diallel") )
graph |
graph |
levGenos |
character vector with the names of the genotypes |
main |
main title for the graph plot |
subplots |
character vector indicating which object(s) to plot |
nothing
Timothee Flutre
nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) plotDesignVarMix(design$graph, levGenos)nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) plotDesignVarMix(design$graph, levGenos)
Plots a diallel matrix.
plotDiallel(diallel, main = NULL)plotDiallel(diallel, main = NULL)
diallel |
matrix |
main |
title |
nothing
Timothee Flutre
nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) plotDiallel(design$diallel)nbGenos <- 25 levGenos <- sprintf(fmt=paste0("geno%0", floor(log10(nbGenos))+1, "i"), 1:nbGenos) nbMixes <- 75 # only binary and balanced design <- getDesignBinaryVarMix(levGenos, nbMixes, seed=12345) plotDiallel(design$diallel)
Computes the relative interaction index (RII) as defined in Armas et al (2004, doi:10.1890/03-0650): RII_i = (Y_i(j) - Y_i) / (Y_i(j) + Y_i) where Y_i(j) is the yield per plant of i when mixed with j, and Y_i is the yield per plant of i when growned in isolation. Missing values are forbidden.
RII( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield", avgIfReps = FALSE )RII( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield", avgIfReps = FALSE )
dat |
data.frame with one yield value per row and with at least four columns: the identifier of each stand, the identifier of the focal genotype or species in each stand, the sowing (plant) proportion of the focal in each stand, and the focal yield in each stand (see the example below); a stand with a single proportion equals to 1 will be interpreted as a monovarietal culture ("pure stand"), a mixture otherwise (and its proportions should sum to 1) |
colIDstand |
column name for the stand identifiers |
colIDfocal |
column name for the focal identifiers |
colProp |
column name for the proportions |
colY |
column name for the yield values |
avgIfReps |
if TRUE, replicates of monovarietals will be averaged; if FALSE, the presence of replicates will return an error |
input data.frame with an extra column named "RII"
Timothee Flutre
sow_density <- 200 (dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield_qt_ha=c(50, 40, 25, 22))) dat$yield_g_m2 <- (dat$yield_qt_ha / 10^4) * 10^6 dat$yield_g_plant <- dat$yield_g_m2 / (sow_density * dat$prop) RII(dat, colY = "yield_g_plant")sow_density <- 200 (dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield_qt_ha=c(50, 40, 25, 22))) dat$yield_g_m2 <- (dat$yield_qt_ha / 10^4) * 10^6 dat$yield_g_plant <- dat$yield_g_m2 / (sow_density * dat$prop) RII(dat, colY = "yield_g_plant")
Computes the net relative interaction index (RII) as defined in Stefan et al (2022): RIInet = sum_i prop_i RII_i. Missing values are forbidden.
RIInet( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colRII = "RII" )RIInet( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colRII = "RII" )
dat |
data.frame with one yield value per row and with at least four columns: the identifier of each stand, the identifier of the focal genotype or species in each stand, the sowing (plant) proportion of the focal in each stand, and the relative interaction index of the focal. |
colIDstand |
column name for the stand identifiers |
colIDfocal |
column name for the focal identifiers |
colProp |
column name for the proportions |
colRII |
column name for the RII values |
data.frame with columns colIDstand and "RIInet"
Timothee Flutre
sow_density <- 200 (dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield_qt_ha=c(50, 40, 25, 22))) dat$yield_g_m2 <- (dat$yield_qt_ha / 10^4) * 10^6 dat$yield_g_plant <- dat$yield_g_m2 / (sow_density * dat$prop) dat <- RII(dat, colY = "yield_g_plant") RIInet(dat)sow_density <- 200 (dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield_qt_ha=c(50, 40, 25, 22))) dat$yield_g_m2 <- (dat$yield_qt_ha / 10^4) * 10^6 dat$yield_g_plant <- dat$yield_g_m2 / (sow_density * dat$prop) dat <- RII(dat, colY = "yield_g_plant") RIInet(dat)
Random generation for the matrix-variate Normal distribution. See https://en.wikipedia.org/wiki/Matrix_normal_distribution.
rmatnorm(n = 1, M, U, V, pivot = c(U = "auto", V = "auto"))rmatnorm(n = 1, M, U, V, pivot = c(U = "auto", V = "auto"))
n |
number of observations |
M |
mean matrix; if missing, will be replaced by a matrix of zeroes; if |
U |
between-row covariance matrix |
V |
between-column covariance matrix |
pivot |
2-element vector with values TRUE/FALSE/"auto", where TRUE (FALSE) means using pivoting (or not) for Choleski decomposition of U and/or V (see |
array
Timothee Flutre
set.seed(1859) Sigma <- matrix(c(3,2,2,4), nrow=2, ncol=2) rho <- Sigma[2,1] / prod(sqrt(diag(Sigma))) samples <- rmatnorm(n=100, M=matrix(0, nrow=10^3, ncol=2), U=diag(10^3), V=Sigma) tmp <- t(apply(samples, 3, function(mat){ c(var(mat[,1]), var(mat[,2]), cor(mat[,1], mat[,2])) })) summary(tmp) # corresponds well to Sigmaset.seed(1859) Sigma <- matrix(c(3,2,2,4), nrow=2, ncol=2) rho <- Sigma[2,1] / prod(sqrt(diag(Sigma))) samples <- rmatnorm(n=100, M=matrix(0, nrow=10^3, ncol=2), U=diag(10^3), V=Sigma) tmp <- t(apply(samples, 3, function(mat){ c(var(mat[,1]), var(mat[,2]), cor(mat[,1], mat[,2])) })) summary(tmp) # corresponds well to Sigma
Computes the "relative yield" (RY) as defined in Fowler (1982, doi:10.2307/2259865) citing de Wit (1960): RY_ij = Y_ij / Y_i where Y_ij is the yield of i when mixed with j and Y_i is the yield of i in monovarietal culture. Note that equation 1 of table 2 in Williams and McCarthy (2001, doi:10.1046/j.1440-1703.2001.00368.x) is called "RY" but in fact it corresponds to "RYP" (relative yield per plant). Note also that the equation for "RY" in Reiss and Drinkwater (2018, doi:10.1002/eap.1629) in fact corresponds to "RYM" (relative yield of mixture). Missing values are forbidden.
RY( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield", avgIfReps = FALSE )RY( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield", avgIfReps = FALSE )
dat |
data.frame with one yield value per row and with at least four columns: the identifier of each stand, the identifier of the focal genotype or species in each stand, the sowing (plant) proportion of the focal in each stand, and the focal yield in each stand (see the example below); a stand with a single proportion equals to 1 will be interpreted as a monovarietal culture ("pure stand") and a mixture otherwise; note that contrary to |
colIDstand |
column name for the stand identifiers |
colIDfocal |
column name for the focal identifiers |
colProp |
column name for the proportions |
colY |
column name for the yield values |
avgIfReps |
if TRUE, replicates of monovarietals will be averaged; if FALSE, the presence of replicates will return an error |
input data.frame with an extra column named "RY"
Timothee Flutre
estimRYRep, RYP, RYT, RYM
(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield=c(50, 40, 25, 22))) RY(dat)(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield=c(50, 40, 25, 22))) RY(dat)
Computes the "relative yield of mixture" (RYM) as defined initially by Wilson (1988, doi:10.2307/2403626) for binary mixtures sowed in equal proportions, and extended by Williams and McCarthy (2001, doi:10.1046/j.1440-1703.2001.00368.x) to binary mixtures of unequal proportions. A RYM above 1 means that the mixture yielded better than the sum of the pure-stand yields, each weighted by their proportion in the mixture. It corresponds to equation 35 in Weigelt and Jolliffe (2003, doi:10.1046/j.1365-2745.2003.00805.x), as well as equation 25b (and not to 25a!) of the same paper where the "control" treatment is the expected mixture yield calculated based on the weighted component monovarietal yields. Note that the RYM was unfortunately called "relative yield" (RY) by Reiss and Drinkwater (2018, doi:10.1002/eap.1629). Missing values are forbidden.
RYM( dat, colIDstand = NULL, colIDcomps = "comps", colProps = NULL, colY = "yield", sep = "-", colOut = "RYM" )RYM( dat, colIDstand = NULL, colIDcomps = "comps", colProps = NULL, colY = "yield", sep = "-", colOut = "RYM" )
dat |
data.frame with one yield value per row and with at least four columns: the identifier of each stand, the identifiers of the components of each stand (separated by a specific symbol in the case of mixtures), the sowing (plant) proportion of the components in each stand (separated by the same symbol in the case of mixtures, and the yield of each stand (see the example below); a stand with a single component should have its proportion equal to 1 and will be interpreted as a monovarietal culture ("pure stand"), a mixture otherwise (and its proportions should sum to 1) |
colIDstand |
column name for the stand identifiers (they should be unique); if NULL, the code will attempt to use |
colIDcomps |
column name for the identifiers of the stand components (separated by |
colProps |
column name for the stand proportions (separated by |
colY |
column name for the yield values |
sep |
separator |
colOut |
column name for the output RYM |
input data.frame with an extra column named from "colOut"
estimRYMRep, RYP, RY, RYT
(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2"), comps=c("geno1", "geno2", "geno1-geno2"), props=c("1", "1", "0.6-0.4"), yield=c(50, 40, 47))) RYM(dat)(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2"), comps=c("geno1", "geno2", "geno1-geno2"), props=c("1", "1", "0.6-0.4"), yield=c(50, 40, 47))) RYM(dat)
Computes the relative yield per plant (RYP) as defined in Fowler (1982, doi:10.2307/2259865): RYP_ij = Y_ij / (p Y_i) where Y_ij is the yield of i when mixed with j at a proportion of p, and Y_i is the yield of i in monovarietal culture. Note that, unfortunately, the equation of RYP was called "relative yield" (RY) in table 2 of Williams and McCarthy (2001, doi:10.1046/j.1440-1703.2001.00368.x), although these authors do have the right equation, i.e., they only dropped the "per plant". Missing values are forbidden.
RYP( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield" )RYP( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield" )
dat |
data.frame with one yield value per row and with at least four columns: the identifier of each stand, the identifier of the focal genotype or species in each stand, the sowing (plant) proportion of the focal in each stand, and the focal yield in each stand (see the example below); a stand with a single proportion equals to 1 will be interpreted as a monovarietal culture ("pure stand"), a mixture otherwise (and its proportions should sum to 1); note that contrary to |
colIDstand |
column name for the stand identifiers |
colIDfocal |
column name for the focal identifiers |
colProp |
column name for the proportions |
colY |
column name for the yield values |
input data.frame with an extra column named "RYP"
Timothee Flutre
(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield=c(50, 40, 25, 22))) RYP(dat)(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield=c(50, 40, 25, 22))) RYP(dat)
Computes the "relative yield total" (RYT) of a mixture from de Wit and Den Bergh (1965) as the sum of all relative yields (RY) of the mixture. It corresponds to equation 26 (based on 25a and not 25b!) in Weigelt and Jolliffe (2003, doi:10.1046/j.1365-2745.2003.00805.x), and is equivalent to the index called "land equivalent ratio" (LER) corresponding to equation 28 of the same paper. Missing values are forbidden.
RYT( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield" )RYT( dat, colIDstand = "ID", colIDfocal = "focal", colProp = "prop", colY = "yield" )
dat |
data.frame with one yield value per row and with at least four columns: the identifier of each stand, the identifier of the focal genotype or species in each stand, the sowing (plant) proportion of the focal in each stand, and the focal yield in each stand (see the example below); a stand with a single proportion equals to 1 will be interpreted as a monovarietal culture ("pure stands") and a mixture otherwise; note that contrary to |
colIDstand |
column name for the stand identifiers |
colIDfocal |
column name for the focal identifiers |
colProp |
column name for the proportions |
colY |
column name for the yield values |
input data.frame with an extra column named "RYT"
(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield=c(50, 40, 25, 22))) RYT(dat)(dat <- data.frame(ID=c("geno1", "geno2", "mixg1g2", "mixg1g2"), focal=c("geno1", "geno2", "geno1", "geno2"), prop=c(1, 1, 0.5, 0.5), yield=c(50, 40, 25, 22))) RYT(dat)
Simulates an incomplete (sparse) yet balanced, tester-based design with respect to a DBV-SBV model for interspecific mixtures.
simulDBVSBVinter( GRMs, parsGenetics = list(mu = list(S1 = c(SC = 65, IC = 32), S2 = c(SC = 30, IC = 27)), sdBlocks = 4, CV_g = c(S1 = 0.08, S2 = 0.08), prop_var_SBV = c(S1 = 0.2, S2 = 0.2), prop_var_SIGV = c(S1 = 0.5, S2 = 0), cor_DBV_SBV = c(S1 = -0.9, S2 = -0.9), prop_var_DBVxSBV = c(S1 = 0, S2 = 0), use_GRM_for_DBVxSBV = TRUE, cor_DxS = 0, H2_SC = c(S1 = 0.7, S2 = 0.7), T2 = c(S1 = 0.7, S2 = 0.7), cor_err_IC = -0.2), parsDesign = list(nbBlocks = 2, nbPlotsPerBl = 500, nbYs = 10, sowingDensities = list(S1 = c(SC = 300, IC = 150), S2 = c(SC = 40, IC = 40)), testersS2 = TRUE, incomplete_balanced = TRUE), sep = "+", seed = NULL )simulDBVSBVinter( GRMs, parsGenetics = list(mu = list(S1 = c(SC = 65, IC = 32), S2 = c(SC = 30, IC = 27)), sdBlocks = 4, CV_g = c(S1 = 0.08, S2 = 0.08), prop_var_SBV = c(S1 = 0.2, S2 = 0.2), prop_var_SIGV = c(S1 = 0.5, S2 = 0), cor_DBV_SBV = c(S1 = -0.9, S2 = -0.9), prop_var_DBVxSBV = c(S1 = 0, S2 = 0), use_GRM_for_DBVxSBV = TRUE, cor_DxS = 0, H2_SC = c(S1 = 0.7, S2 = 0.7), T2 = c(S1 = 0.7, S2 = 0.7), cor_err_IC = -0.2), parsDesign = list(nbBlocks = 2, nbPlotsPerBl = 500, nbYs = 10, sowingDensities = list(S1 = c(SC = 300, IC = 150), S2 = c(SC = 40, IC = 40)), testersS2 = TRUE, incomplete_balanced = TRUE), sep = "+", seed = NULL )
GRMs |
named list with a genomic relationship matrix per species, named S1 and S2; should be a diagonal for S2 if it is a tester (see |
parsGenetics |
list of genetic parameters |
parsDesign |
list of design parameters |
sep |
character separating the names of a mixture components |
seed |
seed for the generation of pseudo-random numbers |
list
Timothee Flutre
## simulate a data set with both sole crops and intercrops: GRMs <- list("S1" = diag(100), "S2" = diag(2)) dimnames(GRMs$S1) <- list(paste0("gS1_", 1:100), paste0("gS1_", 1:100)) dimnames(GRMs$S2) <- list(paste0("gS2_", 1:2), paste0("gS2_", 1:2)) out <- simulDBVSBVinter(GRMs) names(out) str(out$datW) str(out$datL) ## see the third vignette for more details## simulate a data set with both sole crops and intercrops: GRMs <- list("S1" = diag(100), "S2" = diag(2)) dimnames(GRMs$S1) <- list(paste0("gS1_", 1:100), paste0("gS1_", 1:100)) dimnames(GRMs$S2) <- list(paste0("gS2_", 1:2), paste0("gS2_", 1:2)) out <- simulDBVSBVinter(GRMs) names(out) str(out$datW) str(out$datL) ## see the third vignette for more details
Simulates SNP genotypes as allele dose additively encoded, i.e. 0,1,2, using correlated allele frequencies to mimick genetic structure.
simulGenosDoseStruct( nb_genos = c(100, 120, 80), nb_snps = 1000, div_pops = diag(3) * 0.5 + 0.5, geno_IDs = NULL, snp_IDs = NULL )simulGenosDoseStruct( nb_genos = c(100, 120, 80), nb_snps = 1000, div_pops = diag(3) * 0.5 + 0.5, geno_IDs = NULL, snp_IDs = NULL )
nb_genos |
number of genotypes (i.e. individuals) per population |
nb_snps |
number of SNPs |
div_pops |
matrix of divergence among populations, with a diagonal of 1's; the closer off-diagonal values are from 1; the weaker the divergence; the further, the stronger |
geno_IDs |
vector of genotype identifiers (if NULL, will be "geno001", etc) |
snp_IDs |
vector of SNP identifiers (if NULL, will be "snp001", etc) |
matrix with genotypes in rows and SNPs in columns
Timothee Flutre thanks to code from Andres Legarra
## weak divergences among populations: weak_div_pops <- diag(3) 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)] weak_div_pops ## strong divergences among populations: strong_div_pops <- diag(3) strong_div_pops[upper.tri(strong_div_pops)] <- 0.5 strong_div_pops[lower.tri(strong_div_pops)] <- strong_div_pops[upper.tri(strong_div_pops)] strong_div_pops M <- simulGenosDoseStruct(div_pops=weak_div_pops) A <- estimGRM(M) imageMat(A, "Weak divergence") M <- simulGenosDoseStruct(div_pops=strong_div_pops) A <- estimGRM(M) imageMat(A, "Strong divergence")## weak divergences among populations: weak_div_pops <- diag(3) 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)] weak_div_pops ## strong divergences among populations: strong_div_pops <- diag(3) strong_div_pops[upper.tri(strong_div_pops)] <- 0.5 strong_div_pops[lower.tri(strong_div_pops)] <- strong_div_pops[upper.tri(strong_div_pops)] strong_div_pops M <- simulGenosDoseStruct(div_pops=weak_div_pops) A <- estimGRM(M) imageMat(A, "Weak divergence") M <- simulGenosDoseStruct(div_pops=strong_div_pops) A <- estimGRM(M) imageMat(A, "Strong divergence")
Summarizes the GMA-SMA models.
summarizeGMASMAs(fits)summarizeGMASMAs(fits)
fits |
named list of "merMod" objects returned by |
matrix
Timothee Flutre
## see the first vignette## see the first vignette
Applies the vec operator to a matrix, i.e., concatenate its columns, and save the names, too, if any.
vec(mat, sep = "-")vec(mat, sep = "-")
mat |
matrix |
sep |
separator of column and row names; used only if |
vector, possibly with names as <rowname><sep><colname>
Timothee Flutre
mat <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE, dimnames = list(as.character(1:2), letters[1:3])) mat vec(mat)mat <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE, dimnames = list(as.character(1:2), letters[1:3])) mat vec(mat)