Rでt検定をするときに使用する、t.test()は集計前のデータのベクトルを与えることでt検定を行ってくれる（{stats}の）。
ところで、t検定はサンプルの平均値と分散、サンプル数がわかればt値と自由度が分かるので、p値等々ほしい値が分かる。
pythonだとscipy.statsのttest_ind_from_statsとかで記述統計量から計算できるが、ざっと見た感じRで記述統計量からt検定できるものが見つからなかった（多分ちゃんと探せばありそう）。
{stats}のt.test.defulat()を改造すると簡単にできそうだったので作ってみた。
そもそもRで統計処理しているときに記述統計量しかわからないという状況はそんなに多くないと思うけど。
t_test_stats <- function(mean_x, std_x, n_x, mean_y, std_y, n_y, alternative = c("two.sided", "less", "greater"), mu = 0, var.equal = FALSE, conf.level = 0.95){ alternative <- match.arg(alternative) if(length(mean_x) != 1 || length(std_x) != 1 || length(n_x) != 1 || !is.numeric(mean_x) || !is.numeric(std_x) || !is.numeric(n_x)){ stop("'mean_x' and 'std_x', 'n_x' must be a single number") } if(!missing(mean_y) || !missing(std_y) || !missing(n_y)){ if(length(mean_y) != 1 || length(std_y) != 1 || length(n_y) != 1 || !is.numeric(mean_y) || !is.numeric(std_y) || !is.numeric(n_y)){ stop("'mean_y' and 'std_y', 'n_y' must need a single number") } } if(!missing(mu) && (length(mu) != 1 || is.na(mu))) stop("'mu' must be a single number") if(!missing(conf.level) && (length(conf.level) != 1 || !is.finite(conf.level) || conf.level < 0 || conf.level > 1)) stop("'conf.level' must be a single number between 0 and 1") if(n_x < 2) stop("not enough 'x' observations") if(!missing(mean_y) && n_y < 2) stop("not enough 'y' observations") if(std_x < 0) stop("'std_x' must be a possiteve single number") if(!missing(std_y) && std_y < 0) stop("'std_y' must be a possitive single number") if(!missing(mean_y)){ dname <- paste("mean of", deparse(substitute(mean_x)),"and", "mean of", deparse(substitute(mean_y))) }else{ dname <- paste("mean of", deparse(substitute(mean_x))) } var_x <- std_x ^ 2 if(missing(mean_y)){ df <- n_x - 1 stderr <- sqrt(var_x/n_x) if(stderr < 10 *.Machine$double.eps * abs(mean_x)) stop("data are essentially constant") tstat <- (mean_x - mu)/stderr method <- "One Sample t-test" estimate <- mean_x names(estimate) <- "mean of x" }else{ mean_y <- mean_y var_y <- std_y ^ 2 method <- paste(if(!var.equal) "Welch", "Two Sample t-test") estimate <- c(mean_x, mean_y) names(estimate) <- c("mean of x", "mean of y") if (var.equal){ df <- n_x + n_y - 2 v <- 0 if (n_x > 1) v <- v + (n_x - 1) * var_x if (n_y > 1) v <- v + (n_y - 1) * var_y v <- v/df stderr <- sqrt(v * (1/n_x + 1/n_y)) }else{ stderrx <- sqrt(var_x/n_x) stderry <- sqrt(var_y/n_y) stderr <- sqrt(stderrx^2 + stderry^2) df <- stderr^4/(stderrx^4/(n_x - 1) + stderry^4/(n_y - 1)) } if(stderr < 10 *.Machine$double.eps * max(abs(mean_x), abs(mean_y))) stop("data are essentially constant") tstat <- (mean_x - mean_y - mu)/stderr } if(alternative == "less"){ pval <- stats::pt(tstat, df) cint <- c(-Inf, tstat + stats::qt(conf.level, df)) } else if(alternative == "greater"){ pval <- stats::pt(tstat, df, lower.tail = FALSE) cint <- c(tstat - stats::qt(conf.level, df), Inf) }else{ pval <- 2 * stats::pt(-abs(tstat), df) alpha <- 1 - conf.level cint <- stats::qt(1 - alpha/2, df) cint <- tstat + c(-cint, cint) } cint <- mu + cint * stderr names(tstat) <- "t" names(df) <- "df" names(mu) <- if(!missing(mean_y)) "difference in means" else "mean" attr(cint, "conf.level") <- conf.level rval <- list(statistic = tstat, parameter = df, p.value = pval, conf.int = cint, estimate = estimate, null.value = mu, stderr = stderr, alternative = alternative, method = method, data.name = dname) class(rval) <- "htest" return(rval) }
これは私用パッケージとしてGithubに{infunx}として入れておいたので
install.packages("remotes") remotes::install_github("indenkun/infunx") library(infunx)
でもロードできる。
これで、scipy.statsのttest_ind_from_statsのようにt検定ができる。
# aは平均2、標準偏差が1、サンプル数が3のベクトル # bは平均0、標準偏差が1、サンプル数が3のベクトル a <- c(1, 2, 3) b <- c(-1, 0, 1) # 普通のt.test()で2群のt検定をやってみる t.test(a, b)
##
## Welch Two Sample t-test
##
## data: a and b
## t = 2.4495, df = 4, p-value = 0.07048
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2669579 4.2669579
## sample estimates:
## mean of x mean of y
## 2 0
# 同じことをt_test_stats()で記述統計量をつかってやる t_test_stats(mean_x = 2, std_x = 1, n_x = 3, mean_y = 0, std_y = 1, n_y = 3)
##
## Welch Two Sample t-test
##
## data: mean of 2 and mean of 0
## t = 2.4495, df = 4, p-value = 0.07048
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2669579 4.2669579
## sample estimates:
## mean of x mean of y
## 2 0
と同じ結果が出る。
1サンプルのt検定もできる。
t.test(a, mu = 1)
##
## One Sample t-test
##
## data: a
## t = 1.7321, df = 2, p-value = 0.2254
## alternative hypothesis: true mean is not equal to 1
## 95 percent confidence interval:
## -0.4841377 4.4841377
## sample estimates:
## mean of x
## 2
t_test_stats(mean_x = 2, std_x = 1, n_x = 3, mu = 1)
##
## One Sample t-test
##
## data: mean of 2
## t = 1.7321, df = 2, p-value = 0.2254
## alternative hypothesis: true mean is not equal to 1
## 95 percent confidence interval:
## -0.4841377 4.4841377
## sample estimates:
## mean of x
## 2
記述統計量から簡単に検算したいときにつかえるかもしれない。
あとで自分で使うようにメモしておく。
