Nothing
## quadeAllPairsTest.R
##
## Copyright (C) 2015-2019 Thorsten Pohlert
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## A copy of the GNU General Public License is available at
## http://www.r-project.org/Licenses/
#' @name quadeAllPairsTest
#' @title All-Pairs Comparisons for
#' Unreplicated Blocked Data (Quade's All-Pairs Test)
#'
#' @description
#' Performs Quade multiple-comparison test for unreplicated
#' blocked data.
#'
#' @details
#' For all-pairs comparisons of unreplicated blocked data
#' Quade's test can be applied.
#' A total of \eqn{m = k(k-1)/2}
#' hypotheses can be tested. The null hypothesis
#' H\eqn{_{ij}: \theta_i = \theta_j} is tested in the two-tailed test
#' against the alternative
#' A\eqn{_{ij}: \theta_i \ne \theta_j, ~~ i \ne j}.
#'
#' The function has included two methods for approximate p-value estimation:
#' \describe{
#' \item{TDist}{p-values are computed from the t distribution}
#' \item{Normal}{p-values are computed from the standard normal distribution}
#' }
#'
#' If no p-value adjustment is performed (\code{p.adjust.method = "none"}),
#' than a simple protected test is recommended, i.e.
#' all-pairs comparisons should only be applied after a significant
#' \code{\link{quade.test}}. However, any method as implemented in
#' \code{\link{p.adjust.methods}} can be selected by the user.
#'
#' @references
#' W. J. Conover (1999), \emph{Practical nonparametric Statistics},
#' 3rd. Edition, Wiley.
#'
#' N. A. Heckert and J. J. Filliben (2003). NIST Handbook 148:
#' Dataplot Reference Manual, Volume 2: Let Subcommands and Library Functions.
#' National Institute of Standards and Technology Handbook Series, June 2003.
#'
#' D. Quade (1979), Using weighted rankings in the analysis of complete
#' blocks with additive block effects. \emph{Journal of the American
#' Statistical Association}, 74, 680-683.
#'
#' @template class-PMCMR
#'
#' @examples
#' ## Sachs, 1997, p. 675
#' ## Six persons (block) received six different diuretics
#' ## (A to F, treatment).
#' ## The responses are the Na-concentration (mval)
#' ## in the urine measured 2 hours after each treatment.
#' ##
#' y <- matrix(c(
#' 3.88, 5.64, 5.76, 4.25, 5.91, 4.33, 30.58, 30.14, 16.92,
#' 23.19, 26.74, 10.91, 25.24, 33.52, 25.45, 18.85, 20.45,
#' 26.67, 4.44, 7.94, 4.04, 4.4, 4.23, 4.36, 29.41, 30.72,
#' 32.92, 28.23, 23.35, 12, 38.87, 33.12, 39.15, 28.06, 38.23,
#' 26.65),nrow=6, ncol=6,
#' dimnames=list(1:6, LETTERS[1:6]))
#' print(y)
#'
#' ## Global test
#' quade.test(y)
#'
#' ## All-pairs comparisons
#' quadeAllPairsTest(y, dist="TDist", p.adjust.method="holm")
#'
#' @keywords htest nonparametric
#'
#' @seealso
#' \code{\link{quade.test}}, \code{\link{friedmanTest}}
#' @export
quadeAllPairsTest <- function(y, ...)
UseMethod("quadeAllPairsTest")
#' @rdname quadeAllPairsTest
#' @aliases quadeAllPairsTest.default
#' @method quadeAllPairsTest default
#' @template two-way-parms
#' @param dist the test distribution. Defaults to \code{"TDist"}.
#' @param p.adjust.method method for adjusting p values
#' (see \code{\link{p.adjust}}).
#' @importFrom stats pt
#' @importFrom stats pnorm
#' @importFrom stats complete.cases
#' @importFrom stats pairwise.table
#' @importFrom stats p.adjust
#' @importFrom stats p.adjust.methods
#' @export
quadeAllPairsTest.default <-
function(y,
groups,
blocks,
dist = c("TDist", "Normal"),
p.adjust.method = p.adjust.methods,
...)
{
dist <- match.arg(dist)
p.adjust.method <- match.arg(p.adjust.method)
DNAME <- deparse(substitute(y))
if(is.matrix(y)) {
groups <- factor(c(col(y)))
blocks <- factor(c(row(y)))
GRPNAMES <- colnames(y)
}
else {
if(anyNA(groups) || anyNA(blocks))
stop("NA's are not allowed in 'groups' or 'blocks'")
if(any(diff(c(length(y), length(groups), length(blocks))) != 0L))
stop("'y', 'groups' and 'blocks' must have the same length")
DNAME <- paste0(DNAME, ", ",
deparse(substitute(groups)), " and ",
deparse(substitute(blocks)))
if(any(table(groups, blocks) != 1))
stop("not an unreplicated complete block design")
groups <- factor(groups)
blocks <- factor(blocks)
o <- order(groups, blocks)
y <- y[o]
groups <- groups[o]
blocks <- blocks[o]
GRPNAMES <- levels(groups)
}
k <- nlevels(groups)
b <- nlevels(blocks)
## <FIXME split.matrix>
y <- matrix(unlist(split(c(y), blocks)), ncol = k, byrow = TRUE)
y <- y[complete.cases(y), ]
# n <- nrow(y)
r <- t(apply(y, 1L, rank))
q <- rank(apply(y, 1, function(u) max(u) - min(u)))
s <- q * (r - (k + 1) / 2)
w <- q * r
## s is a matrix of ranks within blocks (minus the average rank)
## multiplied by the ranked ranges of the blocks
A <- sum(s ^ 2)
B <- sum(colSums(s) ^ 2) / b
S <- colSums(s)
W <- colSums(w)
if (dist == "TDist") {
METHOD <- "Quade's test with TDist approximation"
DIST <- "t"
denom <- sqrt((2 * b * (A - B)) /
((b - 1) * (k - 1)))
compare.stats <- function(i, j) {
dif <- abs(S[i] - S[j])
tval <- dif / denom
return(tval)
}
PSTAT <- pairwise.table(compare.stats, levels(groups),
p.adjust.method = "none")
compare.levels <- function(i, j) {
dif <- abs(S[i] - S[j])
tval <- dif / denom
pval <- pval <- 2 * pt(abs(tval),
df = (b - 1) * (k - 1),
lower.tail = FALSE)
return(pval)
}
PVAL <- pairwise.table(compare.levels, levels(groups),
p.adjust.method = p.adjust.method)
} else {
METHOD <- paste("Quade's test",
"with standard-normal approximation",
sep = "\t")
DIST <- "z"
n <- b * k
denom <- sqrt((k * (k + 1) * (2 * n + 1) * (k - 1)) /
(18 * n * (n + 1)))
nn <- length(w[, 1])
ff <- 1 / (nn * (nn + 1) / 2)
compare.stats <- function(i, j) {
dif <- abs(W[i] * ff - W[j] * ff)
zval <- dif / denom
return(zval)
}
PSTAT <- pairwise.table(compare.stats, levels(groups),
p.adjust.method = "none")
compare.levels <- function(i, j) {
dif <- abs(W[i] * ff - W[j] * ff)
zval <- dif / denom
pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
return(pval)
}
PVAL <- pairwise.table(compare.levels, levels(groups),
p.adjust.method = p.adjust.method)
}
colnames(PSTAT) <- GRPNAMES[1:(k - 1)]
rownames(PSTAT) <- GRPNAMES[2:k]
colnames(PVAL) <- GRPNAMES[1:(k - 1)]
rownames(PVAL) <- GRPNAMES[2:k]
MODEL <- data.frame(x = y, groups, blocks)
ans <- list(
statistic = PSTAT,
p.value = PVAL,
method = METHOD,
p.adjust.method = p.adjust.method,
data.name = DNAME,
model = MODEL,
dist = DIST
)
class(ans) <- "PMCMR"
ans
}
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.