Nothing
#' Eberhart and Russell's regression model
#' @description
#' `r badge('stable')`
#'
#' Regression-based stability analysis using the Eberhart and Russell (1966) model.
#'
#' @param .data The dataset containing the columns related to Environments, Genotypes,
#' replication/block and response variable(s)
#' @param env The name of the column that contains the levels of the
#' environments.
#' @param gen The name of the column that contains the levels of the genotypes.
#' @param rep The name of the column that contains the levels of the
#' replications/blocks
#' @param resp The response variable(s). To analyze multiple variables in a
#' single procedure use, for example, `resp = c(var1, var2, var3)`.
#' @param verbose Logical argument. If `verbose = FALSE` the code will run
#' silently.
#' @return An object of class `ge_reg` with the folloing items for each
#' variable:
#' * data: The data with means for genotype and environment combinations and the
#' environment index
#' * anova: The analysis of variance for the regression model.
#' * regression: A data frame with the following columns: `GEN`, the genotypes;
#' `b0` and `b1` the intercept and slope of the regression, respectively;
#' `t(b1=1)` the calculated t-value; `pval_t` the p-value for the t test; `s2di`
#' the deviations from the regression (stability parameter); `F(s2di=0)` the
#' F-test for the deviations; `pval_f` the p-value for the F test; `RMSE` the
#' root-mean-square error; `R2` the determination coefficient of the regression.
#' * b0_variance: The variance of b0.
#' * b1_variance: The variance of b1.
#' @md
#' @seealso [superiority()], [ecovalence()], [ge_stats()]
#' @author Tiago Olivoto, \email{tiagoolivoto@@gmail.com}
#' @export
#' @examples
#' \donttest{
#' library(metan)
#'reg <- ge_reg(data_ge2,
#' env = ENV,
#' gen = GEN,
#' rep = REP,
#' resp = PH)
#'plot(reg)
#'
#'}
#' @references Eberhart, S.A., and W.A. Russell. 1966. Stability parameters for comparing Varieties.
#' Crop Sci. 6:36-40. \doi{10.2135/cropsci1966.0011183X000600010011x}
ge_reg = function(.data,
env,
gen,
rep,
resp,
verbose = TRUE){
factors <-
.data %>%
select({{env}}, {{gen}}, {{rep}}) %>%
mutate(across(everything(), as.factor))
vars <-
.data %>%
select({{resp}}, -names(factors)) %>%
select_numeric_cols()
factors %<>% set_names("ENV", "GEN", "REP")
listres <- list()
nvar <- ncol(vars)
if (verbose == TRUE) {
pb <- progress(max = nvar, style = 4)
}
for (var in 1:nvar) {
data <-
factors %>%
mutate(Y = vars[[var]])
if(has_na(data)){
data <- remove_rows_na(data)
has_text_in_num(data)
}
data2 <-
data %>%
mean_by(ENV, GEN) %>%
as.data.frame()
model1 <- lm(Y ~ GEN + ENV + ENV/REP + ENV * GEN, data = data)
modav <- anova(model1)
mydf <-
data %>%
mean_by(GEN, ENV)
iamb <-
data %>%
mean_by(ENV) %>%
add_cols(IndAmb = Y - mean(Y)) %>%
select_cols(ENV, IndAmb)
iamb2 <-
data %>%
mean_by(ENV, GEN) %>%
left_join(iamb, by = "ENV")
meandf <- make_mat(mydf, GEN, ENV, Y) %>% rownames_to_column("GEN")
matx <- make_mat(mydf, GEN, ENV, Y) %>% as.matrix()
iij <- apply(matx, 2, mean, na.rm = TRUE) - mean(matx, na.rm = TRUE)
sumij2 <- sum((iij)^2)
YiIj <- matx %*% iij
if(has_na(matx)){
missing <- which(apply(is.na(matx), 1, function(x){any(x) == TRUE}) == TRUE)
YiIj_complete <- NULL
for(i in seq_along(missing)){
YiIj_complete[i] <- matx[missing[i],][!is.na(matx[missing[i],])] %*% iij[!is.na(matx[missing[i],])]
}
YiIj[which(is.na(YiIj))] <- YiIj_complete
warning("Genotypes ", paste(names(missing), collapse = ", "), " missing in some environments", call. = FALSE)
warning("Regression parameters computed after removing missing values", call. = FALSE)
}
bij <- YiIj/sumij2
svar <- (apply(matx^2, 1, sum, na.rm = TRUE)) - (((apply(matx, 1, sum, na.rm = TRUE))^2)/ncol(matx))
bYijIj <- bij * YiIj
dij <- svar - bYijIj
pred <- apply(matx, 1, mean, na.rm = TRUE) + bij %*% iij
gof <- function(x, y){
R2 = NULL
RMSE = NULL
for (i in 1:nrow(x)){
R2[i] = cor(x[i, ], y[i, ], use = "complete.obs")^2
RMSE[i] = sqrt(sum((x[i, ] - y[i, ])^2, na.rm = TRUE)/ncol(x))
}
return(list(R2 = R2, RMSE = RMSE))
}
gof <- gof(pred, matx)
S2e <- modav$"Mean Sq"[5]
nrep <- length(levels(data$REP))
en <- length(levels(data$ENV))
ge <- length(levels(mydf$GEN))
S2di <- (dij/(en - 2)) - (S2e/nrep)
amod2 <- anova(lm(Y ~ GEN + ENV, data = data2))
# amod2 <- anova(model2)
SSL <- amod2$"Sum Sq"[2]
SSGxL <- amod2$"Sum Sq"[3]
SS.L.GxL <- SSL + SSGxL
SSL.Linear <- (1/length(levels(data$GEN))) * (colSums(matx, na.rm = TRUE) %*% iij)^2/sum(iij^2)
SS.L.GxL.linear <- sum(bYijIj) - SSL.Linear
Df <- c(en * ge - 1,
ge - 1,
ge * (en - 1),
1,
ge - 1,
ge * (en - 2),
replicate(length(dij), en - 2),
en*(nrep - 1) * (ge - 1))
poolerr <- modav$"Sum Sq"[5]/nrep
sigma2 <- modav$"Mean Sq"[5]
dferr <- modav$"Df"[5]
vbo <- sigma2 / (en * nrep)
vb1 <- sigma2 / (nrep * sumij2)
tcal <- (bij - 1) / sqrt(vb1)
ptcal <- 2 * pt(abs(tcal), dferr, lower.tail = FALSE)
SSS <- c(sum(amod2$"Sum Sq"),
amod2$"Sum Sq"[1],
SSL + SSGxL,
SSL.Linear,
SS.L.GxL.linear,
sum(dij),
dij,
poolerr) * nrep
MSSS <- (SSS/Df)
FVAL <- c(NA,
MSSS[2]/MSSS[6],
NA,
NA,
MSSS[5]/MSSS[6],
NA,
MSSS[7:(length(MSSS) - 1)]/MSSS[length(MSSS)],
NA)
PLINES <- 1 - pf(FVAL[7:(length(MSSS) - 1)], Df[7], Df[length(Df)])
pval <- c(NA,
1 - pf(FVAL[2], Df[2], Df[6]),
NA,
NA,
1 - pf(FVAL[5], Df[5], Df[6]),
NA,
PLINES,
NA)
anovadf <- data.frame(Df,
`Sum Sq` = SSS,
`Mean Sq` = MSSS,
`F value` = FVAL,
`Pr(>F)` = pval,
check.names = FALSE)
rownames(anovadf) <- c("Total", "GEN", "ENV + (GEN x ENV)", "ENV (linear)",
" GEN x ENV (linear)", "Pooled deviation",
levels(data$GEN), "Pooled error")
temp <- structure(list(data = iamb2,
anova = as_tibble(rownames_to_column(anovadf, "SV")),
regression = tibble(GEN = levels(mydf$GEN),
b0 = apply(matx, 1, mean, na.rm = TRUE),
b1 = as.numeric(bij),
`t(b1=1)` = tcal,
pval_t = ptcal,
s2di = as.numeric(S2di),
`F(s2di=0)` = FVAL[7:(length(FVAL) - 1)],
pval_f = PLINES,
RMSE = gof$RMSE,
R2 = gof$R2),
bo_variance = vbo,
b1_variance = vb1),
class = "ge_reg")
if (verbose == TRUE) {
run_progress(pb,
actual = var,
text = paste("Evaluating trait", names(vars[var])))
}
listres[[paste(names(vars[var]))]] <- temp
}
return(structure(listres, class = "ge_reg"))
}
#' Plot an object of class ge_reg
#'
#' Plot the regression model generated by the function `ge_reg`.
#'
#'
#' @param x An object of class `ge_factanal`
#' @param var The variable to plot. Defaults to `var = 1` the first
#' variable of `x`.
#' @param type The type of plot to show. `type = 1` produces a plot with
#' the environmental index in the x axis and the genotype mean yield in the y
#' axis. `type = 2` produces a plot with the response variable in the x
#' axis and the slope/deviations of the regression in the y axis.
#' @param plot_theme The graphical theme of the plot. Default is
#' `plot_theme = theme_metan()`. For more details, see
#' [ggplot2::theme()].
#' @param x.lim The range of x-axis. Default is `NULL` (maximum and minimum
#' values of the data set). New arguments can be inserted as `x.lim =
#' c(x.min, x.max)`.
#' @param x.breaks The breaks to be plotted in the x-axis. Default is
#' `authomatic breaks`. New arguments can be inserted as `x.breaks =
#' c(breaks)`
#' @param x.lab The label of x-axis. Each plot has a default value. New
#' arguments can be inserted as `x.lab = "my label"`.
#' @param y.lim The range of x-axis. Default is `NULL`. The same arguments
#' than `x.lim` can be used.
#' @param y.breaks The breaks to be plotted in the x-axis. Default is
#' `authomatic breaks`. The same arguments than `x.breaks` can be
#' used.
#' @param y.lab The label of y-axis. Each plot has a default value. New
#' arguments can be inserted as `y.lab = "my label"`.
#' @param leg.position The position of the legend.
#' @param size.tex.lab The size of the text in the axes text and labels. Default
#' is `12`.
#' @param ... Currently not used..
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @seealso [ge_factanal()]
#' @method plot ge_reg
#' @return An object of class `gg, ggplot`.
#' @export
#' @examples
#' \donttest{
#' library(metan)
#' model <- ge_reg(data_ge2, ENV, GEN, REP, PH)
#' plot(model)
#' }
#'
plot.ge_reg <- function(x,
var = 1,
type = 1,
plot_theme = theme_metan(),
x.lim = NULL,
x.breaks = waiver(),
x.lab = NULL,
y.lim = NULL,
y.breaks = waiver(),
y.lab = NULL,
leg.position = "right",
size.tex.lab = 12,
...){
x <- x[[var]]
if(!type %in% c(1, 2)){
stop("Argument 'type' must be either 1 or 2", call. = FALSE)
}
if(type == 1){
y.lab <- ifelse(missing(y.lab), "Response variable", y.lab)
x.lab <- ifelse(missing(x.lab), "Environmental index", x.lab)
p <-
ggplot(x$data, aes(x = IndAmb, y = Y))+
geom_point(aes(colour = GEN), size = 1.5)+
geom_smooth(aes(colour = GEN), method = "lm", formula = y ~ x, se = FALSE)+
theme_bw()+
labs(x = x.lab, y = y.lab)+
plot_theme %+replace%
theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
axis.title = element_text(size = size.tex.lab, colour = "black"),
axis.ticks = element_line(color = "black"),
axis.ticks.length = unit(.2, "cm"),
legend.position = leg.position)
return(p)
}
if(type == 2){
y.lab <- ifelse(missing(y.lab), "Slope of the regression", y.lab)
x.lab <- ifelse(missing(x.lab), "Response variable", x.lab)
p <-
ggplot(x$regression, aes(x = b0, y = b1))+
geom_point(size = 1.5)+
geom_hline(yintercept = mean(x$regression$b1))+
geom_text_repel(aes(label = GEN))+
labs(x = x.lab, y = y.lab) +
plot_theme
p2 <-
ggplot(x$regression, aes(x = b0, y = s2di))+
geom_point(size = 1.5)+
geom_hline(yintercept = mean(x$regression$s2di))+
geom_text_repel(aes(label = GEN))+
labs(x = x.lab, y = "Deviations from the regression") +
plot_theme
p + p2
}
}
#' Print an object of class ge_reg
#'
#' Print the `ge_reg` object in two ways. By default, the results are shown
#' in the R console. The results can also be exported to the directory into a
#' *.txt file.
#'
#' @param x An object of class `ge_reg`.
#' @param export A logical argument. If `TRUE`, a *.txt file is exported to
#' the working directory.
#' @param file.name The name of the file if `export = TRUE`
#' @param digits The significant digits to be shown.
#' @param ... Options used by the tibble package to format the output. See
#' [`tibble::print()`][tibble::formatting] for more details.
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @method print ge_reg
#' @export
#' @examples
#' \donttest{
#'
#' library(metan)
#' model <- ge_reg(data_ge2, ENV, GEN, REP, PH)
#' print(model)
#' }
print.ge_reg <- function(x, export = FALSE, file.name = NULL, digits = 3, ...) {
opar <- options(pillar.sigfig = digits)
on.exit(options(opar))
if (export == TRUE) {
file.name <- ifelse(is.null(file.name) == TRUE, "ge_reg print", file.name)
sink(paste0(file.name, ".txt"))
}
for (i in 1:length(x)) {
var <- x[[i]]
cat("Variable", names(x)[i], "\n")
cat("---------------------------------------------------------------------------\n")
cat("Joint-regression Analysis of variance\n")
cat("---------------------------------------------------------------------------\n")
print(var$anova)
cat("---------------------------------------------------------------------------\n")
cat("Regression parameters\n")
cat("---------------------------------------------------------------------------\n")
print(var$regression)
cat("---------------------------------------------------------------------------\n")
cat("Variance of b0:", var[["bo_variance"]], "\n")
cat("Variance of b1:", var[["b1_variance"]], "\n")
cat("\n\n\n")
}
if (export == TRUE) {
sink()
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.