| Title: | Genetic Study of Plant Mixtures |
| Version: | 1.0.2 |
| Maintainer: | Timothee Flutre <timothee.flutre@inrae.fr> |
| 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. |
| Depends: | R (≥ 4.0.0), ggplot2, lme4 |
| Imports: | graphics, igraph, MASS, Matrix, methods, Rcpp (≥ 1.0.13), stats, TMB (≥ 1.9.17), utils |
| License: | AGPL-3 |
| Copyright: | INRAE |
| Encoding: | UTF-8 |
| RoxygenNote: | 7.3.3 |
| Suggests: | breedR, INLA, emmeans, knitr, MM4LMM, rmarkdown, testthat |
| Additional_repositories: | https://famuvie.github.io/breedR, https://inla.r-inla-download.org/R/stable |
| VignetteBuilder: | knitr |
| LinkingTo: | TMB, Rcpp, RcppEigen |
| NeedsCompilation: | yes |
| Packaged: | 2026-06-04 19:04:01 UTC; tflutre |
| Author: | Timothee Flutre |
| Repository: | CRAN |
| Date/Publication: | 2026-06-09 15:50:07 UTC |
The plantmix package
Description
plantmix is a package to study plant mixtures.
Introduction
The main goal of the plantmix package is to provide a set of functions for pain-free analysis of both varietal (intraspecific) and crop (interspecific) mixtures, such as optimizing experimental designs, simulating data sets, fitting GMA-SMA and DBV-SBV models, computing diversity indices, etc. Work is ongoing on DBV-SBV models for varietal mixtures, spatial heterogeneity, etc.
Optimizing experimental designs
Fitting GMA-SMA models
Fitting DBV-SBV models
Simulating mixture data
Interaction indices
Miscellaneous
Author(s)
Maintainer: Timothee Flutre timothee.flutre@inrae.fr (ORCID) [contributor]
Other contributors:
INRAE [copyright holder, funder]
Jerome Enjalbert jerome.enjalbert@inrae.fr (ORCID) [contributor]
Emma Forst emma.forst@inra.fr (ORCID) [contributor]
Maxence Remerand maxence.remerand@ens-paris-saclay.fr (ORCID) [contributor]
Jemay Salomon jemaysalomon@gmail.com (ORCID) [contributor]
Change in contribution
Description
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).
Usage
CC(
dat,
colIDstand = "ID",
colIDfocal = "focal",
colProp = "prop",
colY = "yield"
)
Arguments
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 |
Value
input data.frame with an extra column named "CC"
Author(s)
Timothee Flutre
Examples
(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)
Land equivalent ratio (LER)
Description
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.
Usage
LER(x)
Arguments
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 |
Value
list with partial and total LERs
Author(s)
Timothee Flutre
See Also
Examples
(dat <- data.frame(solecrop=c(5, 15),
intercrop=c(4, 9),
row.names=c("grain", "fruit")))
LER(dat)
Relative interaction index (RII)
Description
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.
Usage
RII(
dat,
colIDstand = "ID",
colIDfocal = "focal",
colProp = "prop",
colY = "yield",
avgIfReps = FALSE
)
Arguments
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 |
Value
input data.frame with an extra column named "RII"
Author(s)
Timothee Flutre
See Also
Examples
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")
Net Relative interaction index (RIInet)
Description
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.
Usage
RIInet(
dat,
colIDstand = "ID",
colIDfocal = "focal",
colProp = "prop",
colRII = "RII"
)
Arguments
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 |
Value
data.frame with columns colIDstand and "RIInet"
Author(s)
Timothee Flutre
See Also
Examples
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)
Relative yield (RY)
Description
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.
Usage
RY(
dat,
colIDstand = "ID",
colIDfocal = "focal",
colProp = "prop",
colY = "yield",
avgIfReps = FALSE
)
Arguments
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 |
Value
input data.frame with an extra column named "RY"
Author(s)
Timothee Flutre
See Also
estimRYRep, RYP, RYT, RYM
Examples
(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)
Relative yield of mixture (RYM)
Description
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.
Usage
RYM(
dat,
colIDstand = NULL,
colIDcomps = "comps",
colProps = NULL,
colY = "yield",
sep = "-",
colOut = "RYM"
)
Arguments
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 |
Value
input data.frame with an extra column named from "colOut"
See Also
estimRYMRep, RYP, RY, RYT
Examples
(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)
Relative yield per plant (RYP)
Description
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.
Usage
RYP(
dat,
colIDstand = "ID",
colIDfocal = "focal",
colProp = "prop",
colY = "yield"
)
Arguments
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 |
Value
input data.frame with an extra column named "RYP"
Author(s)
Timothee Flutre
See Also
Examples
(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)
Relative yield total (RYT)
Description
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.
Usage
RYT(
dat,
colIDstand = "ID",
colIDfocal = "focal",
colProp = "prop",
colY = "yield"
)
Arguments
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 |
Value
input data.frame with an extra column named "RYT"
See Also
Examples
(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)
Estimate a GRM
Description
Estimates a genomic relationship matrix with the first estimator from VanRaden (2008).
Usage
estimGRM(M, AFs = NULL)
Arguments
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 |
Value
matrix
Author(s)
Timothee Flutre
See Also
Examples
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")
RYM with "rep" effect
Description
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.
Usage
estimRYMRep(
data,
colY,
colID,
colRep,
mix2genos,
conf_level = 0.95,
plotDiag = FALSE,
title = ""
)
Arguments
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 |
Value
data frame with the inference summaries of RYM per mixture
Author(s)
Timothee Flutre [aut], Arnaud Gauffreteau [ctb]
See Also
Examples
## 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)
RY with "rep" effect
Description
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.
Usage
estimRYRep(
data,
colY,
colIDstand,
colIDfocal,
colRep,
mix2genos,
conf_level = 0.95,
plotDiag = FALSE
)
Arguments
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 |
Value
data frame with the inference summaries of RY per mixture
Author(s)
Timothee Flutre [aut], Arnaud Gauffreteau [ctb]
See Also
Examples
(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")))
Fit DBV-SBV models for interspecific mixtures
Description
Fits DBV-SBV models for interspecific mixtures.
Usage
fitDBVSBVinter(
listY,
listX,
listZ,
listVCov,
REML = TRUE,
lOptions = NULL,
verbose = FALSE
)
Arguments
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 |
Value
list
Author(s)
Jemay Salomon, Timothee Flutre
See Also
Examples
## 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
Fit GMA-SMA models
Description
Fits GMA-SMA models.
Usage
fitGMASMA(
formFix,
data,
listZ,
pkg = "lme4",
listVCov,
contrasts = NULL,
REML = TRUE,
...
)
Arguments
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 |
Value
depends on the chosen package
Author(s)
Timothee Flutre
See Also
Examples
## see the first and second vignettes
Optimize design for binary crop mixtures
Description
Optimizes a design for binary crop mixtures
Usage
getDesignBinaryCropMix(
levGenos1,
levGenos2,
nbMixes,
seed = NULL,
showPlots = FALSE,
verbose = FALSE
)
Arguments
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 |
Value
list
Author(s)
Timothee Flutre
Examples
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)
Optimize design for binary varietal mixtures
Description
Optimizes a design for binary varietal mixtures
Usage
getDesignBinaryVarMix(
levGenos,
nbMixes,
seed = NULL,
showPlots = FALSE,
verbose = FALSE
)
Arguments
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 |
Value
list
Author(s)
Timothee Flutre
Examples
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()
Reformat composition of mixtures into a list
Description
Reformats composition of mixtures into a list.
Usage
getMixtureList(mixtures, sep = NULL)
Arguments
mixtures |
vector, matrix, data.frame or list containing of mixtures |
sep |
required only if |
Value
list of vectors, with one component per mixture, the vector containing the mixture components
Author(s)
Timothee Flutre
See Also
Examples
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)
Get mixture(s) per genotype
Description
Returns the list of mixtures per genotype.
Usage
getMixturesPerGeno(mix2genos)
Arguments
mix2genos |
list with one component per mixture, listing the genotypes it contains (output of |
Value
list with one vector per genotype
Author(s)
Timothee Flutre
See Also
Examples
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))
Image of a matrix
Description
Plots an image of a matrix.
Usage
imageMat(mat, title, col, breaks)
Arguments
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 |
Value
nothing
Author(s)
Timothee Flutre
Examples
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))
Information criterion
Description
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.
Usage
infoCriterion(k, lnLmax, n = NULL, type = "AIC")
Arguments
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 |
Value
numeric with attributes k and n (if not NULL)
Author(s)
Timothee Flutre
Inverse vec operator
Description
Applies the inverse vec operator to a vector.
Usage
invvec(x, n_col, n_row = NULL, sep = NULL)
Arguments
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 |
Value
matrix
Author(s)
Timothee Flutre
See Also
Examples
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 = "-")
Inverse vec for mixtures
Description
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.
Usage
invvecMixes(x, nb_comps)
Arguments
x |
vector; make sure that it is correctly sorted beforehand |
nb_comps |
number of components in each mixture |
Value
matrix with nb_comps columns
Author(s)
Timothee Flutre
Examples
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)
Sowing
Description
Computes the seed weight per genotype for each stand, assuming equal sowing proportions.
Usage
mixSowingWeight(pureSowingWeights, stands)
Arguments
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") |
Value
list which components are stands and values are named vectors with sowing weight per genotype for each stand
Author(s)
Timothee Flutre
See Also
Examples
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)
Design matrices of SMAs for GMA-SMA models
Description
Makes the design matrices with the contribution of specific mixing abilities (SMAs).
More details in the help of mkZSMA and in the vignette.
Usage
mkAllZSMA(
df,
col,
sep,
models = c("2p", "2", "2pp", "3", "3p"),
verbose = FALSE
)
Arguments
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 |
Value
list of SMA design matrices
Author(s)
Timothee Flutre
Examples
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="-"))
Design matrix of GMAs for GMA-SMA models
Description
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.
Usage
mkZGMA(df, col, sep = ",")
Arguments
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 |
Value
matrix
Author(s)
Timothee Flutre
See Also
Examples
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="-"))
Design matrix of SMAs for GMA-SMA models
Description
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.
Usage
mkZSMA(df, col, sep = ",", inc_SMA_ii = "no", skipUnusedCols = TRUE)
Arguments
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 |
Value
matrix
Author(s)
Timothee Flutre
See Also
Examples
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 3
Design matrices for DBV-SBV models in an inter-specific trial
Description
Makes the design matrices for direct and social breeding values per species when the trial includes binary intercrops and, possibly, sole crops.
Usage
mkZinterspe(df, genosPerSp, colIDfocal = "focal", colIDneighbors = "neighbors")
Arguments
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 |
Value
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
Author(s)
Timothee Flutre
Examples
(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"))))
Sowing
Description
Computes the number of seeds of a given genotype to be sowned as pure stand based on a given sowing area and density.
Usage
nbSeedsToSownInPure(sowingArea = 8.4, sowingDensity = 160)
Arguments
sowingArea |
area to be sowned in square meters |
sowingDensity |
number of seeds per square meter |
Value
number of seeds (in thousands)
Author(s)
Timothee Flutre
See Also
Examples
sowingArea <- 9.52
sowingDensity <- 160
nbSeedsToSownInPure(sowingArea, sowingDensity)
sowingArea <- 8.66
sowingDensity <- 200
nbSeedsToSownInPure(sowingArea, sowingDensity)
Normalized bias error
Description
Returns the normalized bias error (NBE).
Usage
normBiasError(estim, truth, perc = TRUE)
Arguments
estim |
numeric vector of estimates |
truth |
numeric vector of true values |
perc |
logical; if TRUE, the return value will be in percentage |
Value
numeric vector
Author(s)
Timothee Flutre
Parametric bootstrap with TMB
Description
Performs parametric bootstrap with TMB and nlminb on a fitted model to allow uncertainty quantification.
Usage
paramBoot4TMB(fit, nb_boot = 5, mc.cores = 1)
Arguments
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) |
Value
list with as many components as nb_boot
Author(s)
Jemay Salomon, Timothee Flutre
Examples
## see the example of `fitDBVSBVinter`
## then run (slow if nb_boot is high!): paramBoot4TMB(fitTmb)
Pivot mixture data into the long format
Description
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).
Usage
pivotMixData2Long(
df,
genos,
colC,
sep = "-",
prefixY = NULL,
sepY = "_",
colIDfocal = "focal",
colIDneighbors = "neighbor"
)
Arguments
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 |
Value
data.frame with more rows
See Also
Examples
## 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")))
Pivot mixture data into the wide format
Description
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.
Usage
pivotMixData2Wide(
df,
colIDstand = "ID",
colIDfocal = "focal",
colIDneighbors = "neighbor",
colPlot = c("x", "y"),
colY = "yield",
sepY = "_",
sepFocalNeighbors = NULL
)
Arguments
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 |
Value
data.frame
Author(s)
Timothee Flutre
See Also
Examples
## 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)
Plot design for crop mixtures
Description
Plots design for crop mixtures, as a graph and as a diallel.
Usage
plotDesignCropMix(graph, levGenos1, levGenos2, main = NULL)
Arguments
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 |
Value
nothing
Author(s)
Timothee Flutre
Examples
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)
Plot design for varietal mixtures
Description
Plots design for varietal mixtures, as a graph and as a diallel.
Usage
plotDesignVarMix(
graph,
levGenos,
main = NULL,
subplots = c("graph", "diallel")
)
Arguments
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 |
Value
nothing
Author(s)
Timothee Flutre
Examples
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)
Plot diallel
Description
Plots a diallel matrix.
Usage
plotDiallel(diallel, main = NULL)
Arguments
diallel |
matrix |
main |
title |
Value
nothing
Author(s)
Timothee Flutre
Examples
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)
Matrix-variate Normal distribution
Description
Random generation for the matrix-variate Normal distribution. See https://en.wikipedia.org/wiki/Matrix_normal_distribution.
Usage
rmatnorm(n = 1, M, U, V, pivot = c(U = "auto", V = "auto"))
Arguments
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 |
Value
array
Author(s)
Timothee Flutre
Examples
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 Sigma
Simulate w.r.t. a DBV-SBV model for interspecific mixtures
Description
Simulates an incomplete (sparse) yet balanced, tester-based design with respect to a DBV-SBV model for interspecific mixtures.
Usage
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
)
Arguments
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 |
Value
list
Author(s)
Timothee Flutre
See Also
Examples
## 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 SNP genotypes
Description
Simulates SNP genotypes as allele dose additively encoded, i.e. 0,1,2, using correlated allele frequencies to mimick genetic structure.
Usage
simulGenosDoseStruct(
nb_genos = c(100, 120, 80),
nb_snps = 1000,
div_pops = diag(3) * 0.5 + 0.5,
geno_IDs = NULL,
snp_IDs = NULL
)
Arguments
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) |
Value
matrix with genotypes in rows and SNPs in columns
Author(s)
Timothee Flutre thanks to code from Andres Legarra
See Also
Examples
## 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")
Model summary
Description
Summarizes the GMA-SMA models.
Usage
summarizeGMASMAs(fits)
Arguments
fits |
named list of "merMod" objects returned by |
Value
matrix
Author(s)
Timothee Flutre
See Also
Examples
## see the first vignette
Vec operator
Description
Applies the vec operator to a matrix, i.e., concatenate its columns, and save the names, too, if any.
Usage
vec(mat, sep = "-")
Arguments
mat |
matrix |
sep |
separator of column and row names; used only if |
Value
vector, possibly with names as <rowname><sep><colname>
Author(s)
Timothee Flutre
See Also
Examples
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)