Pipe operator

Description

See magrittr::%>% for details.

Usage

lhs %>% rhs

Arguments

Value

The result of calling 'rhs(lhs)'.

Check connectivity between trials

Description

Check connectivity between trials

Usage

check_connectivity(
  data = NULL,
  genotype = "line",
  trial = "Experiment",
  response = NULL,
  all = FALSE,
  return_matrix = FALSE
)

Arguments

Value

A data.frame with the genotype connectivity. If return_matrix is TRUE, it will return a n_trial x n_trial matrix with the amount of genotypes shared between each pair of trial.

Examples

library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
head(
  check_connectivity(
    data = dat,
    genotype = "gen",
    trial = "county",
    response = "yield",
    all = TRUE,
    return_matrix = FALSE
  )
)

Check Experimental Design

Description

This function helps to identify the experimental design of each trial, filters the data and then provide a summary for the traits and the experimental design. This works as a quality check before we fit any model. Returns an object of class checkAgri.

Usage

check_design_met(
  data = NULL,
  genotype = NULL,
  trial = NULL,
  traits = NULL,
  rep = NULL,
  block = NULL,
  row = NULL,
  col = NULL
)

Arguments

Value

An object of class checkAgri, with a list of:

Examples

library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
print(results)
plot(results, type = "connectivity")
plot(results, type = "missing")

Correlation Covariance Heatmap

Description

Correlation Covariance Heatmap

Usage

covcor_heat(
  matrix,
  corr = TRUE,
  size = 4,
  digits = 3,
  legend = c(0.6, 0.7),
  upper_tri = FALSE,
  reorder = FALSE
)

Arguments

Value

A ggplot object showing the upper triangular elements of the matrix.

Examples

library(agriutilities)
data(iris)
M <- cor(iris[, -5])
covcor_heat(matrix = M, corr = TRUE)

Extract Residual Variance-Covariance from ASReml-R

Description

This function is specially useful for extracting residual variance covariance matrices from ASReml-R when running repeated measurements analysis.

Usage

extract_rcov(model = NULL, time = NULL, plot = NULL, vc_error = NULL)

Arguments

Details

The expected residual variance covariance structure must be of the form: '~id(Plot):corv(Time)', where 'Plot' is a unique identifier of each experimental unit, and 'Time', represents the variable that contains the time when the experimental units were measured. This form also requires that the levels of the factor Time are nested in the levels of the factor Plot. If it is not in that form you can sort the dataset by using the following command 'arrange(grassUV, Plant, Time)'.

Value

An object with a list of:

Examples

## Not run: 
library(ggpubr)
library(agriutilities)
library(tidyverse)
library(asreml)

head(grassUV)
str(grassUV)

# Exploration -------------------------------------------------------------

grassUV %>%
  ggplot(
    aes(x = Time, y = y, group = Plant, color = Plant)
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(~Tmt) +
  theme_minimal(base_size = 15)

tmp <- grassUV %>%
  group_by(Time, Plant) %>%
  summarise(mean = mean(y, na.rm = TRUE)) %>%
  spread(Time, mean) %>%
  column_to_rownames("Plant")

gg_cor(tmp, label_size = 5)

tmp %>%
  cor(use = "pairwise.complete.obs") %>%
  as.data.frame() %>%
  rownames_to_column(var = "Time") %>%
  gather("DAP2", "corr", -1) %>%
  type.convert(as.is = FALSE) %>%
  mutate(corr = ifelse(Time < DAP2, NA, corr)) %>%
  mutate(DAP2 = as.factor(DAP2)) %>%
  ggplot(
    aes(x = Time, y = corr, group = DAP2, color = DAP2)
  ) +
  geom_point() +
  geom_line() +
  theme_minimal(base_size = 15) +
  color_palette(palette = "jco") +
  labs(color = "Time", y = "Pearson Correlation")

# Modeling ----------------------------------------------------------------

# Identity variance model.
model_0 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):idv(Time),
  data = grassUV
)

# Simple correlation model; homogeneous variance form.
model_1 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):corv(Time),
  data = grassUV
)

# Exponential (or power) model; homogeneous variance form.
model_2 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):expv(Time),
  data = grassUV
)

# Exponential (or power) model; heterogeneous variance form.
model_3 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):exph(Time),
  data = grassUV
)

# Antedependence variance model of order 1
model_4 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):ante(Time),
  data = grassUV
)

# Autoregressive model of order 1; homogeneous variance form.
model_5 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):ar1v(Time),
  data = grassUV
)

# Autoregressive model of order 1; heterogeneous variance form.
model_6 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):ar1h(Time),
  data = grassUV
)

# Unstructured variance model.
model_7 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):us(Time),
  data = grassUV
)

# Model Comparison --------------------------------------------------------

models <- list(
  "id" = model_0,
  "cor" = model_1,
  "exp" = model_2,
  "exph" = model_3,
  "ante" = model_4,
  "ar1" = model_5,
  "ar1h" = model_6,
  "us" = model_7
)

summary_models <- data.frame(
  model = names(models),
  aic = unlist(lapply(models, \(x) summary(x)$aic)),
  bic = unlist(lapply(models, \(x) summary(x)$bic)),
  loglik = unlist(lapply(models, \(x) summary(x)$loglik)),
  nedf = unlist(lapply(models, \(x) summary(x)$nedf)),
  param = unlist(lapply(models, \(x) attr(summary(x)$aic, "param"))),
  row.names = NULL
)

summary_models %>%
  ggplot(
    aes(x = reorder(model, -bic), y = bic, group = 1)
  ) +
  geom_point(size = 2) +
  geom_text(aes(x = model, y = bic + 5, label = param)) +
  geom_line() +
  theme_minimal(base_size = 15) +
  labs(x = NULL)

# Extracting Variance Covariance Matrix -----------------------------------

extract_rcov(model_4)

covcor_heat(
  matrix = extract_rcov(model_1)$corr,
  legend = "none",
  size = 5
) + ggtitle(label = "Uniform Correlation (corv)")
covcor_heat(
  matrix = extract_rcov(model_2)$corr,
  legend = "none",
  size = 5
) + ggtitle(label = "Exponetial (expv)")

## End(Not run)

Extract Variance-Covariance from ASReml-R

Description

Extract Variance-Covariance from ASReml-R

Usage

extract_vcov(model = NULL, gen = "genotype", env = "trial", vc_model = "corv")

Arguments

Value

An object with a list of:

Examples

## Not run: 
library(agridat)
library(agriutilities)

data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
met_results <- met_analysis(out, progress = FALSE)
extract_vcov(
  model = met_results$met_models$yield,
  vc_model = "us"
)

## End(Not run)

Factor Analytic Summary

Description

Factor Analytic Summary

Usage

fa_summary(
  model = NULL,
  trial = "trial",
  genotype = "genotype",
  BLUEs_trial = NULL,
  mult_fa1 = -1,
  mult_fa2 = 1,
  filter_score = 1.5,
  k_biplot = 1,
  size_label_var = 2,
  alpha_label_var = 0.2,
  size_label_ind = 2,
  alpha_label_ind = 0.8,
  size_arrow = 0.2,
  alpha_arrow = 0.2,
  base_size = 12
)

Arguments

Value

An object with a list of:

Examples

## Not run: 
library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
met_results <- met_analysis(out, vcov = "fa2", progress = FALSE)

pp <- met_results$trial_effects
model <- met_results$met_models$yield
fa_summary(
  model = model,
  trial = "trial",
  genotype = "genotype",
  BLUEs_trial = pp,
  mult_fa1 = -1,
  mult_fa2 = -1,
  filter_score = 1,
  k_biplot = 10,
  size_label_var = 3,
  alpha_label_var = 0.5,
  size_label_ind = 3,
  alpha_label_ind = 0.8,
  size_arrow = 0.2,
  alpha_arrow = 0.1
)

## End(Not run)

Triangular Correlation Plot

Description

Return a ggplot object to plot a triangular correlation figure between 2 or more variables.

Usage

gg_cor(
  data,
  colours = c("#db4437", "white", "#4285f4"),
  blackLabs = c(-0.7, 0.7),
  show_signif = TRUE,
  p_breaks = c(0, 0.001, 0.01, 0.05, Inf),
  p_labels = c("***", "**", "*", "ns"),
  show_diagonal = FALSE,
  diag = NULL,
  return_table = FALSE,
  return_n = FALSE,
  adjusted = TRUE,
  label_size = 3,
  method = "pearson"
)

Arguments

Value

A ggplot object containing a triangular correlation figure with all numeric variables in data. If return_table is TRUE, the table used to produce the figure is returned instead.

Author(s)

Daniel Ariza, Johan Aparicio.

Examples

library(agriutilities)
data(iris)
gg_cor(
  data = iris,
  colours = c("#db4437", "white", "#4285f4"),
  label_size = 6
)

Cullis heritability for lme4 models

Description

Cullis heritability for lme4 models

Usage

h_cullis(model, genotype, re_MME = FALSE)

Arguments

Value

A numerical value of the Cullis heritability estimate. If re_MME is TRUE, a list with matrices of the mixed models equations is returned.

Author(s)

Paul Schmidt, Johan Aparicio.

Examples

library(lme4)
library(agridat)
library(agriutilities)
dat <- john.alpha
g.ran <- lmer(
  formula = yield ~ rep + (1 | gen) + (1 | rep:block),
  data = dat
)
h_cullis(model = g.ran, genotype = "gen")

Calculate Cullis heritabilities from SpATS objects

Description

Calculate Cullis heritabilities from SpATS objects

Usage

h_cullis_spt(model)

Arguments

Value

A data frame. The data frame has the following components

• trait : Character string with the trait being analyzed

• H2Cullis : Generalized heritability proposed by Cullis (2006)

• H2Oakey : Generalized heritability proposed by Oakey (2006)

• reBLUP_avg : Average BLUP reliability

• vdBLUP_avg : Average pairwise prediction error variance of genotype effects

• PEV_avg : Average predictive error variance (PEV) of genotype effects

• var_G : Genotypic Variance

Author(s)

Johan Aparicio

References

Cullis, B. R., Smith, A. B., & Coombes, N. E. (2006). On the design of early generation variety trials with correlated data. Journal of agricultural, biological, and environmental statistics, 11, 381-393.

Oakey, H., A. Verbyla, W. Pitchford, B. Cullis, and H. Kuchel (2006). Joint modeling of additive and non-additive genetic line effects in single field trials. Theoretical and Applied Genetics, 113, 809 - 819.

Examples

library(SpATS)
library(agriutilities)
data(wheatdata)
wheatdata$R <- as.factor(wheatdata$row)
wheatdata$C <- as.factor(wheatdata$col)

m1 <- SpATS(
  response = "yield",
  spatial = ~ PSANOVA(col, row, nseg = c(10, 20), nest.div = 2),
  genotype = "geno",
  genotype.as.random = TRUE,
  fixed = ~ colcode + rowcode,
  random = ~ R + C,
  data = wheatdata,
  control = list(tolerance = 1e-03, monitoring = 0)
)

h_cullis_spt(m1)

Heritability for Factor Analytic Models in ASReml-R

Description

Heritability for Factor Analytic Models in ASReml-R

Usage

heritability_fa(
  model_fa = NULL,
  genotype = "line",
  env = "loc",
  vc_model = c("fa2"),
  diag = FALSE
)

Arguments

Value

An object with a list of:

Examples

## Not run: 
library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
met_results <- met_analysis(out, progress = FALSE)
model <- met_results$met_models$yield
heritability_fa(
   model_fa = model,
   genotype = "genotype",
   env = "trial",
   vc_model = "us"
 )

## End(Not run)

Find the AIC and BIC for a set of models fitted using asreml

Description

The function is a wrapper for 'icREML' function described in Verbyla (2019).

Usage

ic_reml_asr(fm, scale = 1)

Arguments

Value

A data frame. The data frame has the following components

• model : the names of the models

• loglik : the full log-likelihood for each model

• p : the number of fixed effects parameters for each model

• q : the number of (non-zero) variance parameters for each model.

• b : the number of variance parameters that are fixed or on the boundary. These parameters are not counted in the AIC or BIC.

• AIC : the AIC for each model

• BIC : the BIC for each model

• logdet : the log-determinant used in adjusting the residual log-likelihood for each model

Author(s)

Ari Verbyla (averbyla at avdataanalytics.com.au)

References

Verbyla, A. P. (2019). A note on model selection using information criteria for general linear models estimated using REML. Australian & New Zealand Journal of Statistics, 61(1), 39-50.

Find the AIC and BIC for a model fitted using SpATS

Description

This function calculates the AIC and BIC for a model fitted in SpATS following the methodology proposed by Verbyla (2019).

Usage

ic_reml_spt(model, scale = 1, k = 2, label = "spats")

Arguments

Value

A data frame. The data frame has the following components

• model : the name of the models

• loglik : the full log-likelihood for each model

• p : the number of fixed effects parameters for each model

• q : the number of (non-zero) variance parameters for each model.

• b : the number of variance parameters that are fixed or on the boundary. These parameters are not counted in the AIC or BIC.

• AIC : the AIC for each model

• BIC : the BIC for each model

• logdet : the log-determinant used in adjusting the residual log-likelihood for each model

Author(s)

Johan Aparicio

References

Verbyla, A. P. (2019). A note on model selection using information criteria for general linear models estimated using REML. Australian & New Zealand Journal of Statistics, 61(1), 39-50.

Examples

library(SpATS)
library(agriutilities)
data(wheatdata)
wheatdata$R <- as.factor(wheatdata$row)
wheatdata$C <- as.factor(wheatdata$col)

m1 <- SpATS(
  response = "yield",
  spatial = ~ PSANOVA(col, row, nseg = c(10, 20), nest.div = 2),
  genotype = "geno",
  genotype.as.random = TRUE,
  fixed = ~ colcode + rowcode,
  random = ~ R + C,
  data = wheatdata,
  control = list(tolerance = 1e-03, monitoring = 0)
)

m2 <- SpATS(
  response = "yield",
  spatial = ~ PSANOVA(col, row, nseg = c(10, 20), nest.div = 2),
  genotype = "geno",
  genotype.as.random = TRUE,
  fixed = ~colcode,
  random = ~ R + C,
  data = wheatdata,
  control = list(tolerance = 1e-03, monitoring = 0)
)

rbind.data.frame(
  ic_reml_spt(m1, label = "colcode_rowcode"),
  ic_reml_spt(m2, label = "colcode_no_rowcode")
)

rbind.data.frame(
  h_cullis_spt(m1),
  h_cullis_spt(m2)
)

Multi-Environmental Trial Analysis

Description

The results of the single_trial_analysis() function are used in met_analysis() to fit multi-environmental trial models. Returns an object of class metAgri, with a list of trial effects, BLUPs, heritability, variance components, stability and the models fitted.

Usage

met_analysis(
  sma_output = NULL,
  h2_filter = 0.2,
  workspace = "1gb",
  vcov = NULL,
  filter_traits = NULL,
  remove_trials = NULL,
  progress = TRUE
)

Arguments

Value

An object of class metAgri, with a list of:

Examples

## Not run: 
library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
met_results <- met_analysis(out, progress = FALSE)
print(met_results)
covcor_heat(matrix = met_results$VCOV$yield$CORR)

## End(Not run)

Genetic Gain Parameters

Description

Genetic Gain Parameters

Usage

parameters_gg(model, trait = "trait")

Arguments

Value

A data.frame with some parameters from the linear regression (Slope, se_Slope, Intercept, r2, Pr(>F)) and the percentage of Genetic Gain.

Examples

library(ggplot2)
library(agridat)
library(magrittr)
library(agriutilities)

data(baker.barley.uniformity)
dat <- baker.barley.uniformity
head(dat)

model <- lm(yield ~ year, dat)

dat %>%
  na.omit() %>%
  ggplot(
    aes(x = year, y = yield)
  ) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_bw()

parameters_gg(model = model, trait = "yield")

Plot an object of class checkAgri

Description

Create several plots for an object of class checkAgri

Usage

## S3 method for class 'checkAgri'
plot(
  x,
  type = c("connectivity", "missing", "boxplot"),
  axis_size = 15,
  text_size = 5,
  ...
)

Arguments

Value

A ggplot object.

Author(s)

Johan Aparicio [aut]

Examples

library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
plot(results, type = "missing")
plot(results, type = "boxplot")

Plot an object of class metAgri

Description

Create several plots for an object of class metAgri

Usage

## S3 method for class 'metAgri'
plot(
  x,
  type = c("correlation", "covariance", "multi_traits"),
  filter_traits = NULL,
  text_size = 4,
  ...
)

Arguments

Value

A ggplot object.

Author(s)

Johan Aparicio [aut]

Examples

## Not run: 
library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
met_results <- met_analysis(out, progress = FALSE)
print(met_results)
plot(met_results, type = "correlation")
plot(met_results, type = "covariance")

## End(Not run)

Plot an object of class smaAgri

Description

Create several plots for an object of class smaAgri

Usage

## S3 method for class 'smaAgri'
plot(
  x,
  type = c("summary", "correlation", "spatial"),
  filter_traits = NULL,
  nudge_y_cv = 3,
  nudge_y_h2 = 0.07,
  horizontal = FALSE,
  theme_size = 15,
  axis_size = 8,
  text_size = 4,
  ...
)

Arguments

Value

A ggplot object.

Author(s)

Johan Aparicio [aut]

Examples

library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
print(out)
plot(out, type = "summary", horizontal = TRUE)
plot(out, type = "correlation")
plot(out, type = "spatial")

Print an object of class checkAgri

Description

Prints information about check_design_met() function.

Usage

## S3 method for class 'checkAgri'
print(x, ...)

Arguments

Value

an object inheriting from class checkAgri.

Author(s)

Johan Aparicio [aut]

Examples

library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
print(results)

Print an object of class metAgri

Description

Prints information about met_analysis() function.

Usage

## S3 method for class 'metAgri'
print(x, ...)

Arguments

Value

an object inheriting from class metAgri.

Author(s)

Johan Aparicio [aut]

Examples

## Not run: 
library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
met_results <- met_analysis(out, progress = FALSE)
print(met_results)

## End(Not run)

Print an object of class smaAgri

Description

Prints information about single_trial_analysis function.

Usage

## S3 method for class 'smaAgri'
print(x, ...)

Arguments

Value

an object inheriting from class smaAgri.

Author(s)

Johan Aparicio [aut]

Examples

library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
print(out)

Single Trial Analysis

Description

The results of the check_design_met() function are used in single_trial_analysis() to fit single trial models. This function can fit, Completely Randomized Designs (CRD), Randomized Complete Block Designs (RCBD), Resolvable Incomplete Block Designs (res-IBD), Non-Resolvable Row-Column Designs (Row-Col) and Resolvable Row-Column Designs (res-Row-Col).

Returns an object of class smaAgri, with a list of trial summary, BLUEs, BLUPs, heritability, variance components, potential extreme observations, residuals, the models fitted and the data used. This function will generate the required output to be used in the two-stage analysis.

Usage

single_trial_analysis(
  results = NULL,
  progress = TRUE,
  engine = "asreml",
  remove_outliers = FALSE
)

Arguments

Value

An object of class smaAgri, with a list of:

Examples

library(agridat)
library(agriutilities)
data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
print(out)

Stability Coefficients

Description

Stability Coefficients

Usage

stability(
  predictions = NULL,
  genotype = NULL,
  trial = NULL,
  response = NULL,
  best = "max"
)

Arguments

Value

A data.frame with several stability measures. "superiority" (cultivar-superiority measure), "static" (Shukla's stability variance) and "wricke" (Wricke's ecovalence).

Examples

## Not run: 
library(agridat)
library(agriutilities)

data(besag.met)
dat <- besag.met
results <- check_design_met(
  data = dat,
  genotype = "gen",
  trial = "county",
  traits = c("yield"),
  rep = "rep",
  block = "block",
  col = "col",
  row = "row"
)
out <- single_trial_analysis(results, progress = FALSE)
met_results <- met_analysis(out, progress = FALSE)

head(
  stability(
    predictions = met_results$BLUPs_GxE,
    genotype = "genotype",
    trial = "trial",
    response = "predicted.value"
  )
)

## End(Not run)