Modeling field trials using statgenSTA

Example

The use of the package is demonstrated using maize data from the European Union project DROPS (https://cordis.europa.eu/project/id/244374). The data is available from https://doi.org/10.15454/IASSTN (E. J. Millet et al. 2019) and the relevant data set is included as data.frame in the statgenSTA package.

The first step is loading the data into R.

data(dropsRaw)

The object named dropsRaw contains data for 256 maize hybrids, grown with two water regimes (irrigated or rainfed), in seven fields in 2012 and 2013. The hybrids are divided in four families. A selection of ten experiments and eight traits is made from the full data set. These ten experiments form a good representation of the full set of experiments covering the five scenarios described in E. Millet et al. (2016). Throughout this vignette in all examples the trait grain.yield will be analyzed. For a more detailed description of the contents of the data see help(dropsRaw).

For the example, first a TD object is created for experiments conducted in 2012. The data for 2013 is then added later on. In practice all this could be done in one go.

## Create a TD object containing data for 2012.
dropsTD <- createTD(data = dropsRaw[dropsRaw$year == 2012, ],
                    genotype = "Variety_ID", 
                    trial = "Experiment",
                    loc = "Site",
                    repId = "Replicate", 
                    subBlock = "block",
                    rowCoord = "Row", 
                    colCoord = "Column", 
                    trLat = "Lat", 
                    trLong = "Long")

The TD object just created is a list with five items, one for each trial (combination of location, year and water regime) in the original data.frame. The column “Variety_ID” in the original data is renamed to “genotype” and converted to a factor. The columns “Replicate”, “block”, “loc” are renamed and converted likewise. The columns “Row” and “Column” are renamed to “rowCoord” and “colCoord” respectively. Simultaneously two columns “rowId” and “colId” are created containing the same information converted to a factor. This seemingly duplicate information is needed for spatial analysis. It is possible to define different columns as “rowId” and “colId” than the ones used as “rowCoord” and “colCoord”. Finally, the latitude and longitude of the trials is read from the column “Lat” and “Long” in the input and stored in a TD object named dropsTD. The information about which columns have been renamed, when creating a TD object, is stored as an attribute of each individual data.frame in the object.

Metadata

The metadata will be a data.frame with five rows, one for each trial in the object named dropsTD. The data.frame has the following columns:

## Extract metadata from the TD object. 
(dropsMeta <- getMeta(TD = dropsTD))
#>        trLocation trDate trDesign  trLat trLong trPlWidth trPlLength
#> Cam12R Campagnola   <NA>       NA 45.171  9.532        NA         NA
#> Cra12R    Craiova   <NA>       NA 44.116 24.348        NA         NA
#> Gai12W    Gaillac   <NA>       NA 43.910  1.890        NA         NA
#> Kar12W  Karlsruhe   <NA>       NA 49.100  8.320        NA         NA
#> Ner12R      Nerac   <NA>       NA 44.170  0.307        NA         NA

After extracting the metadata, it can be modified and then added back to the original TD object.

## Set trial data as 1-1-2012.
dropsMeta$trDate <- as.Date(rep("010112", times = 5), "%d%m%y")
dropsTD <- setMeta(TD = dropsTD, meta = dropsMeta)

Add extra data to a TD object

To add the data for the 2013 trials to the TD object the function addTD() can be used. This function is very similar to createTD(). The only exception is, that a TD object has to be specified to which the new data is added.

## Add the data for the 2013 trials to the TD object.
dropsTD <- addTD(TD = dropsTD, 
                 data = dropsRaw[dropsRaw$year == 2013, ],
                 genotype = "Variety_ID", 
                 trial = "Experiment",
                 loc = "Site", 
                 repId = "Replicate", 
                 subBlock = "block",
                 rowCoord = "Row", 
                 colCoord = "Column", 
                 trLat = "Lat", 
                 trLong = "Long")

## Inspect the metadata after the extra trial was added.
getMeta(TD = dropsTD)
#>         trLocation     trDate trDesign  trLat trLong trPlWidth trPlLength
#> Cam12R  Campagnola 2012-01-01       NA 45.171  9.532        NA         NA
#> Cra12R     Craiova 2012-01-01       NA 44.116 24.348        NA         NA
#> Gai12W     Gaillac 2012-01-01       NA 43.910  1.890        NA         NA
#> Kar12W   Karlsruhe 2012-01-01       NA 49.100  8.320        NA         NA
#> Ner12R       Nerac 2012-01-01       NA 44.170  0.307        NA         NA
#> Kar13R   Karlsruhe       <NA>       NA 49.100  8.320        NA         NA
#> Kar13W   Karlsruhe       <NA>       NA 49.100  8.320        NA         NA
#> Mar13R Martonvasar       <NA>       NA 47.311 18.778        NA         NA
#> Mur13R      Murony       <NA>       NA 46.780 21.050        NA         NA
#> Mur13W      Murony       <NA>       NA 46.780 21.050        NA         NA

Note that the metadata now contains information about all ten trials. Only for the 2012 trials the trial date is filled. Since this information is not used in any of the analyses in the package we leave this.

Summarizing a TD object

The summary() function can be used to get an idea of the content of the data in the TD object. Multiple traits can be summarized at once, but for clarity here a summary is only made for grain.yield of a single trial.

## Create a summary for grain yield in Gai12W.
summary(dropsTD, 
        trial = "Gai12W", 
        traits = "grain.yield")
#> 
#> Summary statistics for grain.yield in Gai12W  
#> 
#>                          grain.yield
#> Number of observations         483.0
#> Number of missing values        19.0
#> Mean                           11.23
#> Median                         11.22
#> Min                             6.71
#> Max                            14.68
#> First quantile                 10.38
#> Third quantile                 12.09
#> Variance                       1.845

Using the default argument settings nine summary statistics are printed, but many more are available. These can be accessed using the argument what in the summary() function. For a full list of available statistics, use help(summary.TD). It is also possible to output all statistics using what = "all".

It is possible to summarize the data in a TD object for different groups. This can be done by using the argument groupBy. It will display three main summary statistics per group. Again, more statistics can be displayed using the argument what.

## Create a summary per family in Gai12W
summary(dropsTD, 
        trial = "Gai12W", 
        traits = "grain.yield", 
        groupBy = "geneticGroup")
#> 
#> Summary statistics for grain.yield in Gai12W grouped by geneticGroup 
#> 
#>                        Stiff_Stalk Iodent Other Lancaster
#> Number of observations       107.0   79.0 210.0      87.0
#> Mean                         11.85  11.42 10.61     11.80
#> Standard deviation           1.320  0.964 1.337     1.110

Plotting a TD object

Several plots can be made to further investigate the contents of a TD object.

• a layout plot showing the layout of the trial
• a map plot showing the location of the trials on a geographical map
• a box plot showing the variability of the traits
• a correlation plot showing the correlation of values for a trait between trials
• a scatter plot showing the scatter plot for a trait between pairs of trials

All plots are described in detail in the sections below.

In all plots the default colors for both genotype groups and trial groups are chosen from a predefined color palette. For genotype groups the color palette is “Dark 2”, for trial groups it is “Alphabet”. See here for an overview of these colors.

It is possible to specify different colors for genotype groups and trial groups per plot using the arguments colGeno and colTrial respectively. Also, more conveniently, the default colors can be set using the options statgen.genoColors and statgen.trialColors.

## Set default colors for genotypes and trials.
options("statgen.genoColors" = c("blue", "green", "yellow"))
options("statgen.trialColors" = c("red", "brown", "purple"))

If a plot has more genotype groups than the number of colors specified as default colors, the default colors will be ignored and the function topo.colors() will be used to generate a color palette instead. For trial groups this is done similarly.

plot(dropsTD, 
     trials = "Gai12W")

## Plot the layout for Gai12W.
## Highlight genotypes A3 and 11430 in red and blue.
plot(dropsTD, 
     trials = "Gai12W", 
     highlight = c("A3", "11430"),
     colHighlight = c("red", "blue"))

## Plot the layout for Gai12W.
## Color sub blocks using polychrome colors for high contrast.
## Colors are specified here since this color palette is only available as such
## from R > 4.0.
cols <- c("#5A5156", "#E4E1E3", "#F6222E", "#FE00FA", "#16FF32", "#3283FE", 
          "#FEAF16", "#B00068", "#1CFFCE", "#90AD1C", "#2ED9FF", "#DEA0FD", 
          "#AA0DFE", "#F8A19F", "#325A9B", "#C4451C", "#1C8356", "#85660D", 
          "#B10DA1", "#FBE426", "#1CBE4F", "#FA0087", "#FC1CBF", "#F7E1A0", 
          "#C075A6", "#782AB6", "#AAF400", "#BDCDFF", "#822E1C", "#B5EFB5", 
          "#7ED7D1", "#1C7F93", "#D85FF7", "#683B79", "#66B0FF", "#3B00FB")
plot(dropsTD, 
     trials = "Gai12W", 
     colorSubBlock = TRUE,
     colSubBlock = cols)

## Plot the layout for Gai12W, label the genotypes.
plot(dropsTD, 
     trials = "Gai12W", 
     showGeno = TRUE)

## Plot the layout for Gai12W, show raw data for grain yield.
plot(dropsTD, 
     trials = "Gai12W", 
     traits = "grain.yield")

## Plot the locations of the trials on a map.
plot(dropsTD, 
     plotType = "map")

## Plot the locations of the trials on a map.
## Color the trials by water scenario.
plot(dropsTD, 
     plotType = "map",
     colorTrialBy = "scenarioWater",
     colTrial = c("red", "darkgreen"))

## Create a boxplot for grain yield.
plot(dropsTD, 
     plotType = "box", 
     traits = "grain.yield")

## Create a boxplot for grain yield with boxes grouped by year 
## Color the boxes by scenario within years.
plot(dropsTD, 
     plotType = "box", 
     traits = "grain.yield", 
     groupBy = "year", 
     colorTrialBy = "scenarioFull")

## Create a correlation plot for grain yield.
plot(dropsTD, 
     plotType = "cor", 
     traits = "grain.yield")

## Create a scatter plot matrix for grain yield.
## Color trials by scenario and genotypes by family.
plot(dropsTD, 
     plotType = "scatter", 
     traits = "grain.yield", 
     colorTrialBy = "scenarioFull", 
     colorGenoBy = "geneticGroup")

Model types

Models can be fitted for five different trial designs. These are listed in the following table with their respective model specifications.

In the models above, fixed effects are indicated in italics, whereas random effects are indicated in bold. The term genotype can be fitted as fixed or as random effect depending on the value of the argument what. Extra fixed effects may be fitted using the argument covariates.
If "SpATS" is used as modeling engine, an extra spatial term is always included in the model (see Spatial models). A spatial term is also included when the modeling engine is "asreml" and the function argument spatial is set to TRUE.

Using the TD object named dropsTD from the previous section, a model for the trial Gai12W and trait grain.yield can now be fitted on the data. The trial was set up using a resolvable row column design. This is specified in fitTD() using the argument design. Since engine is not specified as an argument, "SpATS" is used as engine for fitting the model.

## Fit a single trial model using a model based on a resolvable row column design.
modDropsSp <- fitTD(TD = dropsTD, 
                    trials = "Gai12W", 
                    traits = "grain.yield",
                    design = "res.rowcol")

Note that by not supplying the what argument in the fitTD() function, two models are fitted. One for genotype as a fixed effect and one for genotype as a random effect. The results of both these models are stored in the STA object named modDropsSp. This is very useful for extracting different results from the model later on. A trade-off is that fitting two models takes more time than fitting only one. Therefore, when fitting models to large data sets it is sensible to explicitly define the what argument if only a specific subset of the results is needed as output.

## Fit a single trial model with genotype as random effect.
modDropsSp2 <- fitTD(TD = dropsTD, 
                     trials = "Gai12W", 
                     traits = "grain.yield", 
                     what = "random", 
                     design = "res.rowcol")

Spatial models

When using "SpATS" as a modeling engine for fitting a model, an extra spatial component is always included in the model. This spatial component is composed using the PSANOVA() function in the SpATS package, which uses 2-dimensional smoothing with P-splines as described in Lee, Durbán, and Eilers (2013) and in Rodríguez-Álvarez et al. (2018). See help(PSANOVA, SpATS) for a detailed description. The arguments nseg and nest.div of the PSANOVA() function can be modified using the control argument in the fitTD() function.
The default number for the number of segments is (number of columns / 2, number of rows / 2). Fitting the model in the previous section specifying the number of segments for columns and rows as 28 and 18 respectively, works as follows:

## Fit a spatial single trial model using SpATS. 
## Manually specify the number of segments for rows and columns.
modDropsSp3 <- fitTD(TD = dropsTD, 
                     trials = "Gai12W", 
                     traits = "grain.yield", 
                     design = "res.rowcol", 
                     control = list(nSeg = c(28, 18)))

Alternatively, spatial models can be fitted using "asreml" as modeling engine and setting the argument spatial = TRUE. In this case seven models are fitted and the best model, based on a goodness-of-fit criterion, either AIC or BIC. The default is AIC, this can be modified with the control argument in the fitTD() function.

The seven models fitted largely depend on the trial design. On top of the model described in the previous section extra random terms are added. These extra random terms depend on the structure of the data. If the trial has a regular structure, i.e. all replicates appear the same amount of times in the trial, the following combinations of random and spatial terms are fitted:

If the design is not regular the following combinations of random and spatial terms are fitted:

Fitting a model similar to the one above using asreml with BIC as goodness-of-fit criterion works as follows:

if (requireNamespace("asreml", quietly = TRUE)) {
  ## Fit a spatial single trial model using asreml.
  modDropsAs <- fitTD(TD = dropsTD, 
                      trials = "Gai12W", 
                      traits = "grain.yield", 
                      design = "res.rowcol", 
                      spatial = TRUE, 
                      engine = "asreml",
                      control = list(criterion = "BIC"))
}

The fitted models and the best model are stored in the output together with a summary table with details on the fitted models.

if (exists("modDropsAs")) {
  ## Overview of fitted models.
  print(modDropsAs$Gai12W$sumTab$grain.yield, digits = 2, row.names = FALSE)
}  
#>    spatial random AIC BIC   H2    row   col   error correlated error converge
#>       none   <NA> 610 627 0.83     NA    NA 4.1e-01               NA     TRUE
#>   AR1(x)id   <NA> 611 633 0.83 -0.048    NA 4.0e-01               NA     TRUE
#>   id(x)AR1   <NA> 612 633 0.83     NA 0.039 4.1e-01               NA     TRUE
#>  AR1(x)AR1   <NA> 613 639 0.83 -0.042 0.026 4.1e-01               NA     TRUE
#>   AR1(x)id  units 611 633 0.83 -0.048    NA 6.5e-07            0.404     TRUE
#>   id(x)AR1  units 611 632 0.78     NA 0.402 2.8e-01            0.038     TRUE
#>  AR1(x)AR1  units 582 607 0.84  0.980 0.850 4.0e-01            0.177     TRUE

As the overview shows, the best model, the model with the lowest BIC, is AR1(x)AR1 with units in the random part of the model.

Model summary

Since genotype has been fitted both as fixed and as random factor in the object named modDropsSp, it is possible to calculate both the Best Linear Unbiased Estimators (BLUEs) and the Best Linear Unbiased Predictors (BLUPs). Therefore, both are printed in the summary of the model together with their respective standard errors.

## Set nBest to 5 to limit the output to the best 5 genotypes.
summary(modDropsSp, 
        nBest = 5)
#> Summary statistics 
#> ==================
#> 
#> Summary statistics for grain.yield in Gai12W  
#> 
#>                          grain.yield
#> Number of observations         483.0
#> Number of missing values        19.0
#> Mean                           11.23
#> Median                         11.22
#> Min                             6.71
#> Max                            14.68
#> First quantile                 10.38
#> Third quantile                 12.09
#> Variance                       1.845
#> 
#> 
#> Estimated heritability 
#> ======================
#> 
#> Heritability: 0.84 
#> 
#> Predicted means (BLUEs & BLUPs) 
#> ===============================
#> Best 5 genotypes
#>            BLUEs       SE    BLUPs     SE
#> Lo1261  13.98285  0.47144 13.65848 0.4329
#> Lo1253  13.88760  0.47244 13.44540 0.4350
#> Lo1223  13.71634  0.47173 13.28001 0.4337
#> DKMBST  13.63513  0.47242 13.23844 0.4345
#> HMV5422 13.56305  0.68282 13.02194 0.5785

Model plots

Two types of plots can be made for fitted model objects of class STA.

## Base plots for the model with genotype fitted as random effect.
plot(modDropsSp, 
     plotType = "base", 
     what = "random")

## Spatial plot for the model with genotype fitted as fixed effect.
plot(modDropsSp, 
     plotType = "spatial")

## Spatial plot for the model with genotype fitted as fixed effect.
## Display the spatial trend as a percentage.
plot(modDropsSp, 
     plotType = "spatial", 
     spaTrend = "percentage")

Outlier detection

After fitting a model, it is possible to perform an outlier detection on the results. This outlier detection checks the residuals in the fitted model and compares them to a limit. Observations with a residual larger than this limit are marked as outliers. The default limit is calculated based on the number of observations in the data. Setting a custom limit is also possible.

## Outlier detection for the model with genotype fitted as random.
outliers <- outlierSTA(modDropsSp, 
                       traits = "grain.yield",
                       what = "random")
#> No large standardized residuals.

The output of the outlier detection function is a data.frame with the outliers and of vector of indices of the outliers in the original data. This vector can be used for removing outliers from the data.

It is possible to perform outlier detection on a model with genotype fitted as fixed, however, in doing so it is only possible to detect outliers in genotypes that are replicated in the data. For unreplicated genotypes the residual will always be 0 and, therefore, these genotypes will never be marked as outliers.

To get an idea of the magnitude of the outliers, the argument commonFactors can be used to see the values of observations with similar characteristics. In the example below, for all outliers also the other observations for the same genotype are shown.

## Outlier detection for the model with genotype fitted as random.
## A custom limit is used and commonFactors set to genotype.
outliers <- outlierSTA(modDropsSp, 
                       traits = "grain.yield", 
                       what = "random",
                       rLimit = 2.7, 
                       commonFactors = "genotype")
#> Large standardized residuals.
#> 
#>   trial genotype       trait     value       res outlier
#>  Gai12W     A347 grain.yield  8.767218 -2.766018    TRUE
#>  Gai12W     A347 grain.yield 10.448529  1.853837   FALSE
#>  Gai12W     B105 grain.yield  7.513174 -2.875952    TRUE
#>  Gai12W     B105 grain.yield  9.331618  1.438732   FALSE

Model reports

For objects of class STA there is a report() function available in the statgenSTA package. This function creates a pdf report summarizing the main results of the fitted model. Also, the tex file and figures used for generating the pdf report are saved. By editing the tex file it is possible to modify the report to one’s needs, providing high flexibility.

When no outfile is specified, reports will be created with a default names, e.g. “modelReport_trial_trait_fixed_timestamp.pdf”, in the current working directory. The argument outfile can be used to change the name and location of the report. The value of this argument should be a valid location and name for a pdf file, i.e. including the postfix “.pdf”. Non-existing directories are created by the report() function. When an outfile is specified, trial, trait and modeltype (fixed or random) are concatenated to it in the name of the report.

The reports contain general information on the model fitted, a summary of the results, the plots described in the previous section, a list of best (highest BLUEs or BLUPs) genotypes and a scatter plot of all genotypes and their BLUEs or BLUPs. For some traits a low value might mean a genotype is performing well. To correctly show the best genotypes in the report in this case, set the argument descending to TRUE in the report() function.

## Create a report in the current working directory
report(modDropsSp)
## Create a report for the model with genotype fitted as random.
report(modDropsSp, 
       outfile = "./myReports/dropsReport.pdf", 
       what = "random")

Reporting for an object of class STA can be done for multiple trials, traits and model types simultaneously. For each combination of trial, trait and model type, a separate pdf report will be created. The arguments traits, trials and what can be used for specifying the models for which the reports should be generated.

Adding extra variables

The data.frame named BLUEsDrops consists of only three columns: genotype, trial, and grain.yield. If the model would have been fitted for multiple traits all these traits would become columns in the data.frame. It might be useful to add extra columns from the data used to fit the model to the output for use in further analysis of the data. This can be achieved using the argument keep in the extractSTA() function. To include the water scenario in the output, useful for using the BLUEs for a GxE analysis, use the following command:

## Extract BLUEs from the fitted model.
BLUEsDrops2 <- extractSTA(STA = modDropsSp, 
                          what = "BLUEs", 
                          keep = "scenarioWater")
head(BLUEsDrops2)
#>   genotype  trial scenarioWater grain.yield
#> 1    11430 Gai12W            WW   11.707337
#> 2       A3 Gai12W            WW   11.564098
#> 3     A310 Gai12W            WW   11.056274
#> 4     A347 Gai12W            WW    9.499485
#> 5     A374 Gai12W            WW   10.954845
#> 6     A375 Gai12W            WW   10.808332

Multiple columns can be included in the output in this way. Note that not every column from the original TD object can be included in the extracted data. Only columns that, for each genotype, have only one value per trial can be included. For example, the column repId containing replicates, that has several different values for a single genotype within each trial, cannot be included. When trying to do so it will be ignored with a warning.
It is, however, possible to include repId when extracting fitted values, since for each observation in the original data a fitted value is computed.

## Extract fitted values from the model.
## Add repId and family to the output.
fitVals <- extractSTA(STA = modDropsSp, 
                      what = "fitted", 
                      keep = c("repId", "geneticGroup"))
head(fitVals)
#>   genotype  trial repId geneticGroup grain.yield
#> 1    DE811 Gai12W     1  Stiff_Stalk    9.767671
#> 2  UH_P089 Gai12W     1       Iodent   11.133769
#> 3    W182E Gai12W     1        Other    8.929808
#> 4      Lp5 Gai12W     1  Stiff_Stalk   11.967323
#> 5    F7082 Gai12W     1        Other    8.921387
#> 6   Lo1087 Gai12W     1  Stiff_Stalk   13.213437

Prepare data for GxE analysis

Performing a GxE analysis on can be done using the statgenGxE package. This package uses the same input format for its analyses as the statgenSTA package.

To use the BLUEs or BLUPs from the fitted model in a GxE analysis they can easily be converted into a TD object using the function STAtoTD(). This function creates a TD object from a fitted model including one or more of the following: BLUEs, standard errors of BLUEs, BLUPs and standard errors of BLUPs. Optionally, a column wt with weights (calculated as \((1 / seBLUEs)^2\)) can be added as well. In the same way as described in the previous section extra columns can be added to the output using the argument keep.

## Fit a model for all trials with genotype as fixed factor.
modDropsSpTot <- fitTD(TD = dropsTD, 
                       traits = "grain.yield", 
                       what = "fixed", 
                       design = "res.rowcol")
## Create a TD object containing BLUEs and standard errors of BLUEs.
TDGxE <- STAtoTD(STA = modDropsSpTot, 
                 what = c("BLUEs", "seBLUEs"))
## Add weights and water scenario to the output.
TDGxE2 <- STAtoTD(STA = modDropsSpTot, 
                  what = c("BLUEs", "seBLUEs"), 
                  keep = "scenarioWater",
                  addWt = TRUE)