1# File src/library/stats/R/t.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 19t.test <- function(x, ...) UseMethod("t.test") 20 21t.test.default <- 22function(x, y = NULL, alternative = c("two.sided", "less", "greater"), 23 mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, 24 ...) 25{ 26 alternative <- match.arg(alternative) 27 28 if(!missing(mu) && (length(mu) != 1 || is.na(mu))) 29 stop("'mu' must be a single number") 30 if(!missing(conf.level) && 31 (length(conf.level) != 1 || !is.finite(conf.level) || 32 conf.level < 0 || conf.level > 1)) 33 stop("'conf.level' must be a single number between 0 and 1") 34 if( !is.null(y) ) { 35 dname <- paste(deparse1(substitute(x)),"and", 36 deparse1(substitute(y))) 37 if(paired) 38 xok <- yok <- complete.cases(x,y) 39 else { 40 yok <- !is.na(y) 41 xok <- !is.na(x) 42 } 43 y <- y[yok] 44 } 45 else { 46 dname <- deparse1(substitute(x)) 47 if (paired) stop("'y' is missing for paired test") 48 xok <- !is.na(x) 49 yok <- NULL 50 } 51 x <- x[xok] 52 if (paired) { 53 x <- x-y 54 y <- NULL 55 } 56 nx <- length(x) 57 mx <- mean(x) 58 vx <- var(x) 59 if(is.null(y)) { 60 if(nx < 2) stop("not enough 'x' observations") 61 df <- nx-1 62 stderr <- sqrt(vx/nx) 63 if(stderr < 10 *.Machine$double.eps * abs(mx)) 64 stop("data are essentially constant") 65 tstat <- (mx-mu)/stderr 66 method <- if(paired) "Paired t-test" else "One Sample t-test" 67 estimate <- 68 setNames(mx, if(paired)"mean of the differences" else "mean of x") 69 } else { 70 ny <- length(y) 71 if(nx < 1 || (!var.equal && nx < 2)) 72 stop("not enough 'x' observations") 73 if(ny < 1 || (!var.equal && ny < 2)) 74 stop("not enough 'y' observations") 75 if(var.equal && nx+ny < 3) stop("not enough observations") 76 my <- mean(y) 77 vy <- var(y) 78 method <- paste(if(!var.equal)"Welch", "Two Sample t-test") 79 estimate <- c(mx,my) 80 names(estimate) <- c("mean of x","mean of y") 81 if(var.equal) { 82 df <- nx+ny-2 83 v <- 0 84 if(nx > 1) v <- v + (nx-1)*vx 85 if(ny > 1) v <- v + (ny-1)*vy 86 v <- v/df 87 stderr <- sqrt(v*(1/nx+1/ny)) 88 } else { 89 stderrx <- sqrt(vx/nx) 90 stderry <- sqrt(vy/ny) 91 stderr <- sqrt(stderrx^2 + stderry^2) 92 df <- stderr^4/(stderrx^4/(nx-1) + stderry^4/(ny-1)) 93 } 94 if(stderr < 10 *.Machine$double.eps * max(abs(mx), abs(my))) 95 stop("data are essentially constant") 96 tstat <- (mx - my - mu)/stderr 97 } 98 if (alternative == "less") { 99 pval <- pt(tstat, df) 100 cint <- c(-Inf, tstat + qt(conf.level, df) ) 101 } 102 else if (alternative == "greater") { 103 pval <- pt(tstat, df, lower.tail = FALSE) 104 cint <- c(tstat - qt(conf.level, df), Inf) 105 } 106 else { 107 pval <- 2 * pt(-abs(tstat), df) 108 alpha <- 1 - conf.level 109 cint <- qt(1 - alpha/2, df) 110 cint <- tstat + c(-cint, cint) 111 } 112 cint <- mu + cint * stderr 113 names(tstat) <- "t" 114 names(df) <- "df" 115 names(mu) <- if(paired || !is.null(y)) "difference in means" else "mean" 116 attr(cint,"conf.level") <- conf.level 117 rval <- list(statistic = tstat, parameter = df, p.value = pval, 118 conf.int = cint, estimate = estimate, null.value = mu, 119 stderr = stderr, 120 alternative = alternative, 121 method = method, data.name = dname) 122 class(rval) <- "htest" 123 rval 124} 125 126t.test.formula <- 127function (formula, data, subset, na.action, ...) 128{ 129 if (missing(formula) || (length(formula) != 3L)) 130 stop("'formula' missing or incorrect") 131 oneSampleOrPaired <- FALSE 132 if (length(attr(terms(formula[-2L]), "term.labels")) != 1L) 133 if (formula[[3]] == 1L) 134 oneSampleOrPaired <- TRUE 135 else 136 stop("'formula' missing or incorrect") 137 m <- match.call(expand.dots = FALSE) 138 if (is.matrix(eval(m$data, parent.frame()))) 139 m$data <- as.data.frame(data) 140 ## need stats:: for non-standard evaluation 141 m[[1L]] <- quote(stats::model.frame) 142 m$... <- NULL 143 mf <- eval(m, parent.frame()) 144 DNAME <- paste(names(mf), collapse = " by ") # works in all cases 145 names(mf) <- NULL 146 response <- attr(attr(mf, "terms"), "response") 147 if (! oneSampleOrPaired) { 148 g <- factor(mf[[-response]]) 149 if (nlevels(g) != 2L) 150 stop("grouping factor must have exactly 2 levels") 151 DATA <- setNames(split(mf[[response]], g), c("x", "y")) 152 y <- do.call("t.test", c(DATA, list(...))) 153 if (length(y$estimate) == 2L) { 154 names(y$estimate) <- paste("mean in group", levels(g)) 155 names(y$null.value) <- 156 paste("difference in means between", 157 paste("group", levels(g), collapse = " and ")) 158 } 159 } 160 else { # 1-sample and paired tests 161 respVar <- mf[[response]] 162 if (inherits(respVar, "Pair")){ 163 DATA <- list(x = respVar[,1], y = respVar[,2], paired=TRUE) 164 y <- do.call("t.test", c(DATA, list(...))) 165 } 166 else { 167 DATA <- list(x = respVar) 168 y <- do.call("t.test", c(DATA, list(...))) 169 } 170 } 171 y$data.name <- DNAME 172 y 173} 174