1#  File src/library/stats/R/quade.test.R
2#  Part of the R package, https://www.R-project.org
3#
4#  Copyright (C) 1995-2019 The R Core Team
5#
6#  This program is free software; you can redistribute it and/or modify
7#  it under the terms of the GNU General Public License as published by
8#  the Free Software Foundation; either version 2 of the License, or
9#  (at your option) any later version.
10#
11#  This program is distributed in the hope that it will be useful,
12#  but WITHOUT ANY WARRANTY; without even the implied warranty of
13#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14#  GNU General Public License for more details.
15#
16#  A copy of the GNU General Public License is available at
17#  https://www.R-project.org/Licenses/
18
19quade.test <- function(y, ...) UseMethod("quade.test")
20
21quade.test.default <-
22function(y, groups, blocks, ...)
23{
24    DNAME <- deparse1(substitute(y))
25    if(is.matrix(y)) {
26        d <- dim(y)
27        groups <- factor(.col(d))
28        blocks <- factor(.row(d))
29    }
30    else {
31        if(anyNA(groups) || anyNA(blocks))
32            stop("NA's are not allowed in 'groups' or 'blocks'")
33        if(any(diff(c(length(y), length(groups), length(blocks))) != 0L))
34            stop("'y', 'groups' and 'blocks' must have the same length")
35        DNAME <- paste0(DNAME, ", ",
36                        deparse1(substitute(groups)), " and ",
37                        deparse1(substitute(blocks)))
38        if(any(table(groups, blocks) != 1))
39            stop("not an unreplicated complete block design")
40        ord <- order(groups)
41        y <- y[ord]
42        groups <- factor(groups[ord])
43        blocks <- factor(blocks[ord])
44    }
45    k <- nlevels(groups)
46    b <- nlevels(blocks)
47    ## <FIXME split.matrix>
48    y <- matrix(unlist(split(c(y), blocks)), ncol = k, byrow = TRUE)
49    y <- y[complete.cases(y), ]
50#    n <- nrow(y)
51    r <- t(apply(y, 1L, rank))
52    q <- rank(apply(y, 1, function(u) max(u) - min(u)))
53    s <- q * (r - (k+1)/2)
54    ## S is a matrix of ranks within blocks (minus the average rank)
55    ## multiplied by the ranked ranges of the blocks
56    A <- sum(s^2)
57    B <- sum(colSums(s)^2) / b
58    if(A == B) {
59        ## Treat zero denominator case as suggested by Conover (1999),
60        ## p.374.
61        STATISTIC <- NaN
62        PARAMETER <- c(NA, NA)
63        PVAL <- (gamma(k+1))^(1-b)
64    } else {
65        STATISTIC <- (b - 1) * B / (A - B)
66        ## The same as 2-way ANOVA on the scores S.
67        PARAMETER <- c(k - 1, (b-1) * (k-1))
68        PVAL <- pf(STATISTIC, PARAMETER[1L], PARAMETER[2L], lower.tail = FALSE)
69    }
70    names(STATISTIC) <- "Quade F"
71    names(PARAMETER) <- c("num df", "denom df")
72
73    structure(list(statistic = STATISTIC,
74                   parameter = PARAMETER,
75                   p.value = PVAL,
76                   method = "Quade test",
77                   data.name = DNAME),
78              class = "htest")
79}
80
81quade.test.formula <-
82function(formula, data, subset, na.action, ...)
83{
84    if(missing(formula))
85        stop("'formula' missing")
86    ## <FIXME>
87    ## Maybe put this into an internal rewriteTwoWayFormula() when
88    ## adding support for strata()
89    if((length(formula) != 3L)
90       || (length(formula[[3L]]) != 3L)
91       || (formula[[3L]][[1L]] != as.name("|"))
92       || (length(formula[[3L]][[2L]]) != 1L)
93       || (length(formula[[3L]][[3L]]) != 1L))
94        stop("incorrect specification for 'formula'")
95    formula[[3L]][[1L]] <- as.name("+")
96    ## </FIXME>
97    m <- match.call(expand.dots = FALSE)
98    m$formula <- formula
99    if(is.matrix(eval(m$data, parent.frame())))
100        m$data <- as.data.frame(data)
101    ## need stats:: for non-standard evaluation
102    m[[1L]] <- quote(stats::model.frame)
103    mf <- eval(m, parent.frame())
104    DNAME <- paste(names(mf), collapse = " and ")
105    names(mf) <- NULL
106    y <- do.call("quade.test", as.list(mf))
107    y$data.name <- DNAME
108    y
109}
110