#' Run spatial row and column model
#'
#' @description
#' \code{model_spatial} runs spatial row and column model based on the SpATS package
#'
#' @param data The data frame on which the model is run.
#' It should come from \code{\link{format_data_PPBstats.data_agro}}
#'
#' @param variable variable to analyse
#'
#'
#' @details
#' The model is run with the SpATS function of package SpATS.
#' See ?SpATS for more information.
#'
#' More information can be found in the book \href{https://priviere.github.io/PPBstats_book/family-1.html#spatial-analysis}{here}.
#'
#' @return
#' The function returns a list with two elements :
#' \itemize{
#' \item info : a list with variable
#' \item model :
#' \itemize{
#' \item model: the output from SpATS function
#' \item summary: summary of the model
#' \item df_residuals: effective degree of freedom of the residuals
#' \item MSerror: Mean Square error
#' }
#' }
#'
#' @author
#' Pierre Riviere
#'
#' @seealso
#' \itemize{
#' \item \code{\link{check_model}}
#' \item \code{\link{check_model.fit_model_spatial}}
#' }
#'
#' @importFrom SpATS SpATS
#' @importFrom SpATS summary.SpATS
#' @importFrom methods is
#'
#' @export
#'
model_spatial = function(data, variable){
# 1. Error messages, update arg ----------
if(!is(data, "data_agro")){ stop(substitute(data), " must be formated with type = \"data_agro\", see PPBstats::format_data_PPBstats().") }
check_data_vec_variables(data, variable)
data_tmp = data
colnames(data_tmp)[which(colnames(data_tmp) == variable)] = "variable"
# 2. Set up data set ----------
data_tmp$col = as.numeric(data_tmp$X)
data_tmp$row = as.numeric(data_tmp$Y)
data_tmp$col_f = as.factor(data_tmp$col)
data_tmp$row_f = as.factor(data_tmp$row)
# 3. Run the model ----------
m = suppressMessages(
SpATS::SpATS(
response = "variable",
genotype = "germplasm",
genotype.as.random = TRUE,
spatial = ~ PSANOVA(col, row, nseg = c(nlevels(data_tmp$X), nlevels(data_tmp$Y))),
random = ~ col_f + row_f,
data = data_tmp)
)
# 4. MSerror for mean comparisons ----------
s = summary.SpATS(m, which = "dimensions")
deviance_model = deviance(m)
df_residual = as.numeric(s$p.table.dim["Residual", "Effective"])
MSerror = deviance_model / df_residual
# 5. Return results ----------
out = list(
"info" = list("variable" = variable, "data" = data),
"model" = list(
"model" = m,
"summary" = summary(m),
"df_residual" = df_residual,
"MSerror" = MSerror
)
)
class(out) <- c("PPBstats", "fit_model_spatial")
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.