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