R: Penalized factorial regression using glmnet

GnE {factReg}

R Documentation

Penalized factorial regression using glmnet

Description

Based on multi-environment field trials, fits the factorial regression model $Y_{ij} = \mu + e_j + g_i + \sum_{k=1}^s \beta_{ik} x_{ij} + \epsilon_{ij},$ with environmental main effects $e_j$, genotypic main effects $g_{i}$ and genotype-specific environmental sensitivities $\beta_{ik}$. See e.g. Millet et al 2019 and Bustos-Korts et al 2019. There are $s$ environmental indices with values $x_{ij}$. Optionally, predictions can be made for a set of test environments, for which environmental indices are available. The new environments must contain the same set of genotypes, or a subset.

Penalization: the model above is fitted using glmnet, simultaneously penalizing $e_j$, $g_i$ and $\beta_{ik}$. If penG = 0 and penE = 0, the main effects $g_i$ and $e_j$ are not penalized. If these parameters are 1, the the main effects are penalized to the same degree as the sensitivities. Any non negative values are allowed. Cross validation is performed with each fold containing a number of environments (details below).

After fitting the model, it is possible to replace the estimated genotypic main effects and sensitivities by their predicted genetic values. Specifically, if a kinship matrix K is assigned the function performs genomic prediction with g-BLUP for the genotypic main effect and each of the sensitivities in turn.

Predictions for the test environments are first constructed using the estimated genotypic main effects and sensitivities; next, predicted environmental main effects are added. The latter are obtained by regressing the estimated environmental main effects for the training environments on the average values of the indices in these environments, as in Millet et al. 2019.

Usage

GnE(
  dat,
  Y,
  G,
  E,
  K = NULL,
  indices = NULL,
  indicesData = NULL,
  testEnv = NULL,
  weight = NULL,
  outputFile = NULL,
  corType = c("pearson", "spearman"),
  partition = data.frame(),
  nfolds = 10,
  alpha = 1,
  lambda = NULL,
  penG = 0,
  penE = 0,
  scaling = c("train", "all", "no"),
  quadratic = FALSE,
  verbose = FALSE
)

Arguments

`dat`	A `data.frame` with data from multi-environment trials. Each row corresponds to a particular genotype in a particular environment. The data do not need to be balanced, i.e. an environment does not need to contain all genotypes. `dat` should contain the training as well as the test environments (see testEnv)
`Y`	The trait to be analyzed: either of type character, in which case it should be one of the column names in `dat`, or numeric, in which case the Yth column of `dat` will be analyzed.
`G`	The column in `dat` containing the factor genotype (either character or numeric).
`E`	The column in `dat` containing the factor environment (either character or numeric).
`K`	A kinship matrix. Used for replacing the estimated genotypic main effect and each of the sensitivities by genomic prediction from a g-BLUP model for each of them. If `NULL`, the estimated effects from the model are returned and used for constructing predictions.
`indices`	The columns in `dat` containing the environmental indices (vector of type character). Alternatively, if the indices are always constant within environments (i.e. not genotype dependent), the environmental data can also be provided using the argument `indicesData` (see below).
`indicesData`	An optional `data.frame` containing environmental indices (covariates); one value for each environment and index. It should have the environment names as row names (corresponding to the names contained in `dat$E`); the column names are the indices. If `indices` (see before) is also provided, the latter will be ignored.
`testEnv`	vector (character). Data from these environments are not used for fitting the model. Accuracy is evaluated for training and test environments separately. The default is `NULL`, i.e. no test environments, in which case the whole data set is training. It is also possible that there are test environments, but without any data; in this case, no accuracy is reported for test environments (CHECK correctness).
`weight`	Numeric vector of length `nrow(dat)`, specifying the weight (inverse variance) of each observation, used in glmnet. Default `NULL`, giving constant weights.
`outputFile`	The file name of the output files, without .csv extension which is added by the function. If not `NULL` the prediction accuracies for training and test environments are written to separate files. If `NULL` the output is not written to file.
`corType`	type of correlation: Pearson (default) or Spearman rank sum.
`partition`	`data.frame` with columns E and partition. The column E should contain the training environments (type character); partition should be of type integer. Environments in the same fold should have the same integer value. Default is `data.frame()`, in which case the function uses a leave-one-environment out cross-validation. If `NULL`, the (inner) training sets used for cross-validation will be drawn randomly from all observations, ignoring the environment structure. In the latter case, the number of folds (nfolds) can be specified.
`nfolds`	Default `NULL`. If `partition == NULL`, this can be used to specify the number of folds to be used in glmnet.
`alpha`	Type of penalty, as in glmnet (1 = LASSO, 0 = ridge; in between = elastic net). Default is 1.
`lambda`	Numeric vector; defines the grid over which the penalty lambda is optimized in cross validation. Default: NULL (defined by glmnet). Important special case: lambda = 0 (no penalty).
`penG`	numeric; default 0. If 1, genotypic main effects are penalized. If 0, they are not. Any non negative real number is allowed.
`penE`	numeric; default 0. If 1, environmental main effects are penalized. If 0, they are not. Any non negative real number is allowed.
`scaling`	determines how the environmental variables are scaled. "train" : all data (test and training environments) are scaled using the mean and and standard deviation in the training environments. "all" : using the mean and standard deviation of all environments. "no" : No scaling.
`quadratic`	boolean; default `FALSE`. If `TRUE`, quadratic terms (i.e., squared indices) are added to the model.
`verbose`	boolean; default `FALSE`. If `TRUE`, the accuracies per environment are printed on screen.

Value

A list with the following elements:

predTrain: A data.frame with predictions for the training set
predTest: A data.frame with predictions for the test set
resTrain: A data.frame with residuals for the training set
resTest: A data.frame with residuals for the test set
mu: the estimated overall mean
envInfoTrain: The estimated environmental main effects, and the predicted effects, obtained when the former are regressed on the averaged indices, using penalized regression
envInfoTest: The predicted environmental main effects for the test environments, obtained from penalized regression using the estimated main effects for the training environments and the averaged indices
parGeno: data.frame containing the estimated genotypic main effects (first column) and sensitivities (subsequent columns)
trainAccuracyEnv: a data.frame with the accuracy (r) for each training environment, as well as the root mean square error (RMSE), mean absolute deviation (MAD) and rank (the latter is a proportion: how many of the best 5 genotypes are in the top 10). To be removed or further developed. All these quantities are also evaluated for a model with only genotypic and environmental main effects (columns rMain, RMSEMain and rankMain)
testAccuracyEnv: A data.frame with the accuracy for each test environment, with the same columns as trainAccuracyEnv
trainAccuracyGeno: a data.frame with the accuracy (r) for each genotype, averaged over the training environments
testAccuracyGeno: a data.frame with the accuracy (r) for each genotype, averaged over the test environments
lambda: The value of lambda selected using cross validation
lambdaSequence: The values of lambda used in the fits of glmnet. If lambda was provided as input, the value of lambda is returned
RMSEtrain: The root mean squared error on the training environments
RMSEtest: The root mean squared error on the test environments
Y: The name of the trait that was predicted, i.e. the column name in dat that was used
G: The genotype label that was used, i.e. the argument G that was used
E: The environment label that was used, i.e. the argument E that was used
indices: The indices that were used, i.e. the argument indices that was used
quadratic: The quadratic option that was used

References

Millet, E.J., Kruijer, W., Coupel-Ledru, A. et al. Genomic prediction of maize yield across European environmental conditions. Nat Genet 51, 952–956 (2019). doi:10.1038/s41588-019-0414-y

Examples

## load the data, which are contained in the package
data(drops_GE)
data(drops_GnE)

## We remove identifiers that we don't need.
drops_GE_GnE <- rbind(drops_GE[, -c(2, 4)], drops_GnE[, -c(2, 4)])

## Define indeces.
ind <- colnames(drops_GE)[13:23]

## Define test environments.
testenv <- levels(drops_GnE$Experiment)

## Additive model, only main effects (set the penalty parameter to a large value).
Additive_model <- GnE(drops_GE_GnE, Y = "grain.yield", lambda = 100000,
                      G = "Variety_ID", E = "Experiment", testEnv = testenv,
                      indices = ind, penG = FALSE, penE = FALSE,
                      alpha = 0.5, scaling = "train")

## Full model, no penalization (set the penalty parameter to zero).
Full_model <- GnE(drops_GE_GnE, Y = "grain.yield", lambda = 0,
                  G = "Variety_ID", E = "Experiment", testEnv = testenv,
                  indices = ind, penG = FALSE, penE = FALSE,
                  alpha = 0.5, scaling = "train")

## Elastic Net model, set alpha parameter to 0.5.
Elnet_model <- GnE(drops_GE_GnE, Y = "grain.yield", lambda = NULL,
                   G = "Variety_ID", E = "Experiment", testEnv = testenv,
                   indices = ind, penG = FALSE, penE = FALSE,
                   alpha = 0.5, scaling = "train")

## Lasso model, set alpha parameter to 1.
Lasso_model <- GnE(drops_GE_GnE, Y = "grain.yield", lambda = NULL,
                   G = "Variety_ID", E = "Experiment", testEnv = testenv,
                   indices = ind, penG = FALSE, penE = FALSE,
                   alpha = 1, scaling = "train")

## Ridge model, set alpha parameter to 0.
Ridge_model <- GnE(drops_GE_GnE, Y = "grain.yield", lambda = NULL,
                   G = "Variety_ID", E = "Experiment", testEnv = testenv,
                   indices = ind, penG = FALSE, penE = FALSE,
                   alpha = 0, scaling = "train")

[Package factReg version 1.0.0 Index]