Nothing
#' Multiple comparison procedures to the means.
#'
#' \code{MCPtest} applies several tests of multiple comparisons present
#' in the literature. The tests chosen are based on the evaluation
#' of the researcher to decide the choice of test for analysis
#' in the experiment.
#' @param y Model (aov or lm), numeric vector containing the
#' response variable or the mean of the treatments.
#' @param trt Constant (y = model) or a vector containing
#' the treatments.
#' @param dferror Degrees of freedom of the Mean Square Error.
#' @param mserror Mean Square Error.
#' @param replication Number de repetitions of the treatments in the experiment.
#' For unbalanced data should be informed the harmonic mean of repetitions.
#' This argument should be informed only if \code{ismean = TRUE}.
#' @param alpha Significant level. The default is \code{alpha = 0.05}.
#' @param main Title of the analysis.
#' @param MCP Allows choosing the multiple comparison test;
#' the \emph{defaut} is "all". This option will go perform all tests.
#' However, the options are: the Skott-Knott midrange test ("MGM"),
#' the Skott-Knott Range test ("MGR"), the Tukey midrange test ("TM"),
#' the Scott-Knott's test ("SK").
#' @param ismean Logic. If \code{FALSE} (default), the \code{y} argument represents
#' a model (aov or lm) or a numeric vector containing the response variable.
#' If \code{TRUE} the \code{y} argument represents the mean of treatments.
#' @param parallel Logic. If \code{FALSE} (default), the \code{MCPtest} function
#' runs without parallelization. Otherwise, the code is executed with
#' parallelization. Note that the parallelization is not always more efficient.
#' @return \code{MCPtest} returns the print of a list of results. First,
#' the summary of \code{y}. Second, the statistics
#' of the test chosen. And finally, the mean group results for each test.
#' If \code{MCPtest} function is stored
#' in an object, the results will be printed and
#' also stored in the object.
#' @details The \code{MCP} argument allows you to choose various tests
#' of multiple comparisons at once. For example,
#' \code{MCP = c("MGM", "MGR")}, and so on.
#' @references
#' BATISTA, Ben Deivide de Oliveira. Testes de comparacoes multiplas
#' baseados na distribuicao da \emph{midrange} \emph{estudentizada}
#' externamente. 2016. 194f. Tese (Doutorado em Estatistica e
#' Experimentacao Agropecuaria) - Universidade Federal de Lavras, 2016.
#' (Portuguese, Brazil)
#'
#' SCOTT, A. J.; KNOTT, M. A cluster analysis method for grouping means in
#' the analysis of variance. Biometrics, International Biometric Society,
#' v. 30, n. 3, p. 507-512, 1974.
#'
#' DUNCAN, D. B. Multiple range and multiple F Tests. Biometrics, International
#' Biometric Society, v. 11, n. 1,p. 1-42, 1955.
#' @examples
#' # Simulated data (completely randomized design)
#'
#' # Response variable
#' rv <- c(100.08, 105.66, 97.64, 100.11, 102.60, 121.29, 100.80,
#' 99.11, 104.43, 122.18, 119.49, 124.37, 123.19, 134.16,
#' 125.67, 128.88, 148.07, 134.27, 151.53, 127.31)
#'
#' # Treatments
#' treat <- factor(rep(LETTERS[1:5], each = 4))
#'
#' # Anova
#' res <- anova(aov(rv~treat))
#' DFerror <- res$Df[2]
#' MSerror <- res$`Mean Sq`[2]
#'
#' # Loading the MCPtests package
#' library(MCPtests)
#'
#' # applying the tests
#' results <- MCPtest(y = rv,
#' trt = treat,
#' dferror = DFerror,
#' mserror = MSerror,
#' alpha = 0.05,
#' main = "Multiple Comparison Procedure: MGM test",
#' MCP = c("MGM"))
#'
#' # Other option for the MCP argument is "all". All tests are used.
#'
#' results$Groups # Results of the tests
#' results$Statistics # Main arguments of the tests
#' results$Summary # Summary of the response variable
#'
#' # Using the y argument as aov or lm model
#' res <- aov(rv~treat)
#'
#' MCPtest(y = res, trt = "treat", alpha = 0.05,
#' main = "Multiple Comparison Procedure: MGM test",
#' MCP = c("MGM"))
#'
#' # For unbalanced data: It will be used the harmonic mean of
#' # the number of experiment replicates
#'
#' # Using the previous example
#' rv <- rv[-1]
#' treat <- treat[-1]
#'
#' res <- lm(rv~treat) # Linear model
#'
#' # Multiple comparison procedure: MGR test
#' MCPtest(y = res, trt = "treat", alpha = 0.05,
#' main = "Multiple Comparison Procedure: MGR test",
#' MCP = c("MGR"))
#'
#' # Assuming that the available data are the averages
#' # of the treatments and the analysis of variance
#'
#' # Analysis of Variance Table
#'
#' # Response: rv
#' # Df Sum Sq Mean Sq F value Pr(>F)
#' # treat 4 4135.2 1033.80 14.669 4.562e-05 ***
#' # Residuals 15 1057.1 70.47
#'
#' mean.treat <- c(100.87, 105.95, 117.62, 127.97, 140.30)
#' treat <- factor(LETTERS[1:5])
#' DFerror <- 15
#' MSerror <- 70.47488
#' replic <- 4
#'
#' MCPtest(y = mean.treat,
#' trt = treat,
#' dferror = DFerror,
#' mserror = MSerror,
#' replication = replic,
#' alpha = 0.05,
#' main = "Multiple Comparison Procedure: MGM test",
#' MCP = c("MGM"),
#' ismean = TRUE)
#'
#' @import "utils" "graphics" "SMR" "foreach"
#' @importFrom "stats" "deviance" "df.residual" "qtukey" "aov" "pchisq" "sd"
#' @importFrom "parallel" "makeCluster" "detectCores" "stopCluster"
# @importFrom "doParallel" "registerDoParallel"
#' @export
MCPtest <- function(y, trt = NULL, dferror = NULL, mserror = NULL, replication = NULL, alpha = 0.05, main = NULL,
MCP = "all", ismean = FALSE, parallel = FALSE){
#####################################################
#Defensive programming
if (ismean == TRUE) {
if (is.null(replication)) {
stop("The replication argument must be informed", call. = FALSE)
}
}
if (is.numeric(y)) {
if (!is.numeric(y)) {
stop("The y argument must be numeric", call. = FALSE)
}
if (!is.factor(trt)) {
stop("The trt argument must be factor", call. = FALSE)
}
if (length(y) != length(trt)) {
stop("The y and trt arguments must have same length", call. = FALSE)
}
if (is.null(dferror)) {
stop("The dferror argument must be informed", call. = FALSE)
}
if (is.null(mserror)) {
stop("The mserror argument must be informed", call. = FALSE)
}
}
#####################################################
if (all(MCP == "all")){
MCP = c("MGM", "MGR", "SNKM", "TM", "SK") # FAZER ALTERACOES
}
#####################################################
#Defensive programming
if (is.null(trt)) {
stop("The trt argument is required", call. = FALSE)
}
mcps <- c("MGM", "MGR", "SNKM", "TM", "SK") # FAZER ALTERACOES
nas <- pmatch(MCP, mcps)
if (any(is.na(nas))) {
stop("The options for the MCP argument are 'MGM', 'MGR', 'SNKM', 'TM' and 'SK'", call. = FALSE) # FAZER ALTERACOES
}
################################################
name.y <- paste(deparse(substitute(y)))
name.trt <- paste(deparse(substitute(trt)))
if (is.null(main)) {
main <- paste(name.y, "~", name.trt)
}
#y: "aov" or "lm" class
if ("aov" %in% class(y) | "lm" %in% class(y)) {
if (is.null(main)) {
main <- y$call
}
A <- y$model
dferror <- df.residual(y)
mserror <- deviance(y)/dferror
y <- A[, 1]
if (length(trt) > length(names(A)[-1]) + 1) {
stop("The length of the trt argument is invalid", call. = FALSE)
}
ipch <- pmatch(trt, names(A))
if (any(is.na(ipch))) {
ncolnam <- length(colnames(A))
optrt <- colnames(A)[2]
if (ncolnam > 2) {
for (i in 3:ncolnam) {
optrt <- paste(optrt, names(A)[i], sep = " ")
}
}
stop("The options for the trt argument are\n Options: ", optrt, "\n Note: Observe which of the options have signified practice", call. = FALSE)
}
ipch <- pmatch(trt, names(A))
nipch <- length(ipch)
ntrt <- length(names(A))
if (any(is.na(ipch))) {
optrt <- names(A)[2]
for (i in 3:ntrt) {
optrt <- paste(optrt, names(A)[i], sep = " ")
}
stop("Any of the options of trt argument is wrong \n Options: ", optrt, call. = FALSE)
}
name.t <- names(A)[ipch][1]
trt <- A[, ipch]
if (nipch > 1) {
trt <- A[, ipch[1]]
for (i in 2:nipch) {
name.t <- paste(name.t, names(A)[ipch][i], sep = ":")
trt <- paste(trt, A[, ipch[i]], sep = ":")
}
}
name.y <- names(A)[1]
}
#if (console) cat(gettext("Multiple Comparison Procedures\n", domain = "R-MCP"))
#if (console) cat(gettext("Study: ", domain = "R-MCP"), main, "\n\n")
# Mean standard error
ep <- function(x){
epm <- sd(x)/sqrt(length(x))
return(epm)
}
if (ismean == TRUE) {
data.na <- subset(data.frame(y, trt), is.na(y) == FALSE)
means <- round(y, 2)
n <- length(y)
rn <- rep(replication, n)
rh <- 1/mean(1/rn)
Mean <- mean(y)
CV <- sqrt(mserror) * 100/Mean
} else {
data.na <- subset(data.frame(y, trt), is.na(y) == FALSE)
means <- round(tapply(data.na[, 1], data.na[, 2], "mean"), 2)
std <- round(tapply(data.na[, 1], data.na[, 2], "ep"), 2)
rn <- tapply(data.na[, 1], data.na[, 2], "length")
mi <- round(tapply(data.na[, 1], data.na[, 2], "min"), 2)
ma <- round(tapply(data.na[, 1], data.na[, 2], "max"), 2)
n <- nrow(means)
rh <- 1/mean(1/rn)
Mean <- mean(data.na[, 1])
CV <- sqrt(mserror) * 100/Mean
}
################################################
#Observation for unbalanced data
if (length(unique(rn)) != 1) {
#if (console) cat(gettext("Unbalanced data: It will be used the harmonic mean of \n the number of experiment replicates \n", domain = "R-MCP"))
}
################################################
if (ismean == TRUE) {
summarydata <- data.frame(Means = means,
r = rh)
#if (console) cat(gettext("Summary:\n", domain = "R-MCP"))
colnames(summarydata) <- c(gettext("Means", domain = "R-MCP"),
"r"
)
#if (console) cat(summarydata)
} else {
summarydata <- data.frame(Means = means,
std = std,
r = rh,
Min = mi,
Max = ma)
#if (console) cat(gettext("Summary:\n", domain = "R-MCP"))
colnames(summarydata) <- c(gettext("Means", domain = "R-MCP"),
gettext("std", domain = "R-MCP"),
"r", "Min", "Max"
)
#if (console) print(summarydata)
}
if (length(unique(rn)) != 1) {
#if (console) cat(gettext("\n Harmonic mean of the number of experiment replicates",
# domain = "R-MCP"), rh, "\n")
}
#DMS midrange
if (any(MCP == "SNKM")) {
nnn <- 2:n
aaa <- rep(1 - alpha / 2, times = (n - 1))
MRq.snkm <- SMR::qSMR(aaa, nnn, dferror) # Studentized midrange of the SNKM test
MRq <- MRq.snkm[n - 1] # Studentized midrange of the MGM/TM test
dms.int <- MRq * sqrt(mserror / rh) # Internal DMS of the MGM/TM test
dms <- dms.int + sqrt(0.5) * sqrt(mserror / rh) / n^0.5 # add of mean square error/nmed^0.5
dms <- c(dms, dms.int) # External DMS of the MGM/TM test
dmsint <- MRq.snkm * sqrt(mserror / rh) # Internal DMS of the SNKM test
dms.snkm <- dmsint + sqrt(0.5) * sqrt(mserror / rh) / n^0.5 # External DMS of the SNKM test
} else{
MRq <- SMR::qSMR(1 - alpha / 2, n, dferror) # Studentized midrange of the MGM/TM test
dms.int <- MRq * sqrt(mserror / rh) # Internal DMS of the MGM/TM test
dms <- dms.int + sqrt(0.5) * sqrt(mserror/rh) / n^0.5 # External DMS of the MGM/TM test
dms <- c(dms, dms.int) # DMS of the MGM/TM test
}
#DMS range:
Rq <- qtukey(1 - alpha, n, dferror) # Studentized range of the MGR test
dms.range <- Rq * sqrt(mserror / rh) # DMS of the MGR test
#Initial statistics:
statistics.MGM <- NA
statistics.MGR <- NA
statistics.SNKM <- NA
statistics.TM <- NA
statistics.SK <- NA
#Initial groups:
test.MGM <- NA
test.MGR <- NA
test.SNKM <- NA
test.TM <- NA
test.SK <- NA
#Defensive programming
if (!any(parallel == c(TRUE,FALSE))) {
stop("The options for the parallel argument are TRUE or FALSE!", call. = FALSE)
}
# Means Grouping Midrange Test (BATISTA, 2016)
if (any(MCP == "MGM")){
#if (console) cat(gettext("\nMeans Grouping Midrange Test (BATISTA, 2016)\n\n", domain = "R-MCP"))
statistics <- data.frame(Exp.Mean = Mean,
CV = CV,
MSerror = mserror,
DF = dferror,
n = n,
Stud.Midrange = MRq,
Ext.DMS = dms[1],
Int.DMS = dms[2])
#if(console) cat(gettext("Statistics: \n", domain = "R-MCP"))
rownames(statistics) <- " "
colnames(statistics) <- c(gettext("Exp.Mean", domain = "R-MCP"),
"CV",
gettext("MSerror", domain = "R-MCP"),
gettext("DF", domain = "R-MCP"),
"n",
gettext("Stud.Midrange", domain = "R-MCP"),
gettext("Ext.DMS", domain = "R-MCP"),
gettext("Int.DMS", domain = "R-MCP")
)
statistics.MGM <- statistics
#if (console) print(statistics)
test <- MGMtest(y, trt, n, dferror, mserror, alpha, dms)
test[1] <- round(test[1], 2)
test.MGM <- test
#if(console) cat(gettext("\nGroups: \n", domain = "R-MCP"))
#if (console) print(test)
}
# Means Grouping based on the Range Test (BATISTA, 2016)
if (any(MCP == "MGR")){
#if (console) cat(gettext("\nMean Grouping Range Test (BATISTA, 2016)\n\n", domain = "R-MCP"))
statistics <- data.frame(Exp.Mean = Mean,
CV = CV,
MSerror = mserror,
DF = dferror,
n = n,
Stud.Range = Rq,
DMS = dms.range)
#if (console) cat(gettext("Statistics: \n", domain = "R-MCP"))
rownames(statistics) <- " "
colnames(statistics) <- c(gettext("Exp.Mean", domain = "R-MCP"),
"CV",
gettext("MSerror", domain = "R-MCP"),
gettext("DF", domain = "R-MCP"),
"n",
gettext("Stud.Range", domain = "R-MCP"),
"DMS"
)
statistics.MGR <- statistics
#if (console) print(statistics)
test <- MGRtest(y, trt, n, dferror, mserror, alpha, dms.range)
test[1] <- round(test[1], 2)
test.MGR <- test
#if (console) cat(gettext("\nGroups: \n", domain = "R-MCP"))
#if (console) print(test)
}
# Student-Newman-Keuls (SNK) Midrange Test (BATISTA, 2016)
if (any(MCP == "SNKM")){
#if (console) cat(gettext("\nSNK Midrange Test (BATISTA, 2016)\n\n", domain = "R-MCP"))
statistics <- data.frame(Exp.Mean = Mean,
CV = CV,
MSerror = mserror,
DF = dferror,
n = nnn,
Stud.Midrange = MRq.snkm,
DMS = dms.snkm)
statistics <- statistics[order(statistics[, 6], decreasing = FALSE),]
cont <- seq(1, n - 1, by = 1)
rownames(statistics) <- cont
rownames(statistics) <- rownames(MRq.snkm, do.NULL = FALSE, prefix = "comp")
#if (console) cat(gettext("Statistics: \n", domain = "R-MCP"))
colnames(statistics) <- c(gettext("Exp.Mean", domain = "R-MCP"),
"CV",
gettext("MSerror", domain = "R-MCP"),
gettext("DF", domain = "R-MCP"),
"n",
gettext("Stud.Midrange", domain = "R-MCP"),
"DMS"
)
statistics.SNKM <- round(statistics, 4)
#if (console) print(round(statistics, 4))
test <- SNKMtest(y, trt, n, dferror, mserror, alpha, dms.snkm)
test[1] <- round(test[1], 2)
test.SNKM <- test
#if (console) cat(gettext("\nGroups: \n", domain = "R-MCP"))
#if (console) print(test)
}
# Tukey Midrange Test (BATISTA, 2016)
if (any(MCP == "TM")) {
#if (console) cat(gettext("\nTukey Midrange Test (BATISTA, 2016)\n\n", domain = "R-MCP"))
statistics <- data.frame(Exp.Mean = Mean,
CV = CV,
MSerror = mserror,
DF = dferror,
n = n,
Stud.Midrange = MRq,
Ext.DMS = dms[1],
Int.DMS = dms[2])
#if (console) cat(gettext("Statistics: \n", domain = "R-MCP"))
rownames(statistics) <- " "
colnames(statistics) <- c(gettext("Exp.Mean", domain = "R-MCP"),
"CV",
gettext("MSerror", domain = "R-MCP"),
gettext("DF", domain = "R-MCP"),
"n",
gettext("Stud.Midrange", domain = "R-MCP"),
gettext("Ext.DMS", domain = "R-MCP"),
gettext("Int.DMS", domain = "R-MCP")
)
statistics.TM <- statistics
#if (console) print(statistics)
test <- TMtest(y, trt, n, dferror, mserror, alpha, dms)
test[1] <- round(test[1], 2)
test.TM <- test
#if (console) cat(gettext("\nGroups: \n", domain = "R-MCP"))
#if (console) print(test)
}
# Scott-Knott's test (SCOTT-KNOTT, 1974)
if (any(MCP == "SK")) {
#if (console) cat(gettext("\nScott-Knott's Test (SCOTT-KNOTT, 1974)\n\n", domain = "R-MCP"))
statistics <- data.frame(Exp.Mean = Mean,
CV = CV,
MSerror = mserror,
DF = dferror,
n = n
)
#if (console) cat(gettext("Statistics: \n", domain = "R-MCP"))
rownames(statistics) <- " "
colnames(statistics) <- c(gettext("Exp.Mean", domain = "R-MCP"),
"CV",
gettext("MSerror", domain = "R-MCP"),
gettext("DF", domain = "R-MCP"),
"n"
)
statistics.SK <- statistics
#if (console) print(statistics)
test <- sktest(y, trt, dferror, mserror, rh, alpha, parallel)
test.SK <- test
#if (console) cat(gettext("\nGroups: \n", domain = "R-MCP"))
#if (console) print(test)
}
#All statistics
stat.tests <- list(Statistics.MGM = statistics.MGM,
Statistics.MGR = statistics.MGR,
Statistics.SNKM = statistics.SNKM,
Statistics.TM = statistics.TM,
Statistics.SK = statistics.SK)
#All groups
group.tests <- list(group.MGM = test.MGM,
group.MGR = test.MGR,
group.SNKM = test.SNKM,
group.TM = test.TM,
group.SK = test.SK)
################
# Output results
################
alltest <- 1:5
nas <- nas[order(nas, na.last = NA)]
ntest <- alltest[-nas]
if (length(nas) == 5) {
statistics <- stat.tests
grouptest <- group.tests
} else {
statistics <- stat.tests[- ntest]
grouptest <- group.tests[- ntest]
}
# if(length(nas) == 4) {
# grouptest <- group.tests
# } else{
# grouptest <- group.tests[- ntest]
# }
output <- list(Summary = summarydata,
Groups = grouptest,
Statistics = statistics,
Tests = MCP)
#invisible(output)
return(output)
}
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.