# File src/library/stats/R/oneway.test.R # Part of the R package, https://www.R-project.org # # Copyright (C) 1995-2015 The R Core Team # # 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 2 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 # https://www.R-project.org/Licenses/ oneway.test <- function(formula, data, subset, na.action, var.equal = FALSE) { if(missing(formula) || (length(formula) != 3L)) stop("'formula' missing or incorrect") dp <- as.character(formula) if(length(dp) != 3L) stop("a two-sided formula is required") DNAME <- paste(dp[[2L]], "and", dp[[3L]]) m <- match.call(expand.dots = FALSE) if(is.matrix(eval(m$data, parent.frame()))) m$data <- as.data.frame(data) m$var.equal <- NULL ## need stats:: for non-standard evaluation m[[1L]] <- quote(stats::model.frame) mf <- eval(m, parent.frame()) response <- attr(attr(mf, "terms"), "response") y <- mf[[response]] if(length(mf[-response]) > 1L) g <- factor(do.call("interaction", mf[-response])) else g <- factor(mf[[-response]]) k <- nlevels(g) if(k < 2L) stop("not enough groups") n.i <- tapply(y, g, length) if(any(n.i < 2)) stop("not enough observations") m.i <- tapply(y, g, mean) v.i <- tapply(y, g, var) w.i <- n.i / v.i sum.w.i <- sum(w.i) tmp <- sum((1 - w.i / sum.w.i)^2 / (n.i - 1)) / (k^2 - 1) METHOD <- "One-way analysis of means" if(var.equal) { n <- sum(n.i) STATISTIC <- ((sum(n.i * (m.i - mean(y))^2) / (k - 1)) / (sum((n.i - 1) * v.i) / (n - k))) PARAMETER <- c(k - 1, n - k) PVAL <- pf(STATISTIC, k - 1, n - k, lower.tail = FALSE) } else { ## STATISTIC <- sum(w.i * (m.i - mean(y))^2) / ## ((k - 1) * (1 + 2 * (k - 2) * tmp)) m <- sum(w.i * m.i) / sum.w.i STATISTIC <- sum(w.i * (m.i - m)^2) / ((k - 1) * (1 + 2 * (k - 2) * tmp)) PARAMETER <- c(k - 1, 1 / (3 * tmp)) PVAL <- pf(STATISTIC, k - 1, 1 / (3 * tmp), lower.tail = FALSE) METHOD <- paste(METHOD, "(not assuming equal variances)") } names(STATISTIC) <- "F" names(PARAMETER) <- c("num df", "denom df") RVAL <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = DNAME) class(RVAL) <- "htest" RVAL }