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