Rで統計的仮説検定:ウェルチのt検定

Rでデータサイエンス

ウェルチのt検定

  • 2群の平均値の差の検定
  • 2群に対応なし。
  • 2群の分散に等分散が仮定できない。
  • 帰無仮説:2群の平均値は同一である。

サンプル

seed <- 20230502
sample_n_small <- 10
sample_n_big <- 40
sd_x <- 5
sd_y <- 15
# 平均値が異なる
set.seed(seed = seed)
x <- rnorm(n = sample_n_big, mean = 10, sd = sd_x)
y <- rnorm(n = sample_n_big, mean = 20, sd = sd_y)
t.test(x = x, y = y, alternative = "two.sided", paired = F, var.equal = F, conf.level = 0.95)
# 帰無仮説は棄却される。

    Welch Two Sample t-test

data:  x and y
t = -5.0664, df = 48.943, p-value = 0.000006172
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -16.542457  -7.146214
sample estimates:
mean of x mean of y 
 9.304418 21.148754 
# 平均値が同一
set.seed(seed = seed)
x <- rnorm(n = sample_n_big, mean = 10, sd = sd_x)
y <- rnorm(n = sample_n_big, mean = 10, sd = sd_y)
t.test(x = x, y = y, alternative = "two.sided", paired = F, var.equal = F, conf.level = 0.95)
# 帰無仮説は棄却されない。

    Welch Two Sample t-test

data:  x and y
t = -0.78892, df = 48.943, p-value = 0.434
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -6.542457  2.853786
sample estimates:
mean of x mean of y 
 9.304418 11.148754 
# 平均値が異なる。但しサンプルサイズを10とする。
set.seed(seed = seed)
x <- rnorm(n = sample_n_small, mean = 10, sd = sd_x)
y <- rnorm(n = sample_n_small, mean = 20, sd = sd_y)
t.test(x = x, y = y, alternative = "two.sided", paired = F, var.equal = F, conf.level = 0.95)
# 本サンプルの場合、平均値は異なるが帰無仮説が棄却されない。

    Welch Two Sample t-test

data:  x and y
t = -1.5542, df = 11.576, p-value = 0.147
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -21.738594   3.680391
sample estimates:
mean of x mean of y 
 9.715514 18.744615 

シミュレーション

# 関数
fun_sim_non_normal <- function(n_x, n_y, location_x, location_y, scale_x, scale_y, iter) {
    t_test_p_value <- vector()
    for (iii in seq(iter)) {
        x <- rlogis(n = n_x, location = location_x, scale = scale_x)
        y <- rlogis(n = n_y, location = location_y, scale = scale_y)
        t_test_p_value[iii] <- t.test(x = x, y = y, alternative = "two.sided", paired = F, var.equal = F, conf.level = 0.95)$p.value
    }
    return(sum(t_test_p_value >= 0.05)/iter)
}
fun_sim_normal <- function(n_x, n_y, mean_x, mean_y, sd_x, sd_y, iter) {
    t_test_p_value <- vector()
    for (iii in seq(iter)) {
        x <- rnorm(n = n_x, mean = mean_x, sd = sd_x)
        y <- rnorm(n = n_y, mean = mean_y, sd = sd_y)
        t_test_p_value[iii] <- t.test(x = x, y = y, alternative = "two.sided", paired = F, var.equal = F, conf.level = 0.95)$p.value
    }
    return(sum(t_test_p_value >= 0.05)/iter)
}
iter <- 5000

正規分布に従う場合

サンプルサイズが10の場合
平均値は異なる場合
set.seed(seed = seed)
n_x <- sample_n_small
n_y <- sample_n_small
mean_x <- 10
mean_y <- 20
fun_sim_normal(n_x = n_x, n_y = n_y, mean_x = mean_x, mean_y = mean_y, sd_x = sd_x, sd_y = sd_y, iter = iter)
[1] 0.5516
平均値が同一の場合
set.seed(seed = seed)
n_x <- sample_n_small
n_y <- sample_n_small
mean_x <- 10
mean_y <- 10
1 - fun_sim_normal(n_x = n_x, n_y = n_y, mean_x = mean_x, mean_y = mean_y, sd_x = sd_x, sd_y = sd_y, iter = iter)
[1] 0.0496
サンプルサイズが40の場合
平均値は異なる場合
set.seed(seed = seed)
n_x <- sample_n_big
n_y <- sample_n_big
mean_x <- 10
mean_y <- 20
fun_sim_normal(n_x = n_x, n_y = n_y, mean_x = mean_x, mean_y = mean_y, sd_x = sd_x, sd_y = sd_y, iter = iter)
[1] 0.0248
平均値が同一の場合
set.seed(seed = seed)
n_x <- sample_n_big
n_y <- sample_n_big
mean_x <- 10
mean_y <- 10
1 - fun_sim_normal(n_x = n_x, n_y = n_y, mean_x = mean_x, mean_y = mean_y, sd_x = sd_x, sd_y = sd_y, iter = iter)
[1] 0.0532

正規分布に従わない場合

  • 本ページではロジスティック分布を例とする。
  • 平均:位置パラメータ(location)
  • 分散:\(\dfrac{\pi^2\,\textrm{scale}^2}{3},\quad\)scale:尺度(scale)パラメータ
scale_x <- 5
scale_y <- 8
サンプルサイズが10の場合
平均値は異なる場合
set.seed(seed = seed)
n_x <- sample_n_small
n_y <- sample_n_small
location_x <- 10
location_y <- 20
fun_sim_non_normal(n_x = n_x, n_y = n_y, location_x = location_x, location_y = location_y, scale_x = scale_x, scale_y = scale_y, iter = iter)
[1] 0.5728
平均値が同一の場合
set.seed(seed = seed)
n_x <- sample_n_small
n_y <- sample_n_small
location_x <- 10
location_y <- 10
1 - fun_sim_non_normal(n_x = n_x, n_y = n_y, location_x = location_x, location_y = location_y, scale_x = scale_x, scale_y = scale_y, iter = iter)
[1] 0.045
サンプルサイズが40の場合
平均値は異なる場合
set.seed(seed = seed)
n_x <- sample_n_big
n_y <- sample_n_big
location_x <- 10
location_y <- 20
fun_sim_non_normal(n_x = n_x, n_y = n_y, location_x = location_x, location_y = location_y, scale_x = scale_x, scale_y = scale_y, iter = iter)
[1] 0.0498
平均値が同一の場合
set.seed(seed = seed)
n_x <- sample_n_big
n_y <- sample_n_big
location_x <- 10
location_y <- 10
1 - fun_sim_non_normal(n_x = n_x, n_y = n_y, location_x = location_x, location_y = location_y, scale_x = scale_x, scale_y = scale_y, iter = iter)
[1] 0.0508

コード

methods("t.test")
[1] t.test.default* t.test.formula*
see '?methods' for accessing help and source code
getS3method(f = "t.test", class = "default")
function (x, y = NULL, alternative = c("two.sided", "less", "greater"), 
    mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, 
    ...) 
{
    alternative <- match.arg(alternative)
    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 (!is.null(y)) {
        dname <- paste(deparse1(substitute(x)), "and", deparse1(substitute(y)))
        if (paired) 
            xok <- yok <- complete.cases(x, y)
        else {
            yok <- !is.na(y)
            xok <- !is.na(x)
        }
        y <- y[yok]
    }
    else {
        dname <- deparse1(substitute(x))
        if (paired) 
            stop("'y' is missing for paired test")
        xok <- !is.na(x)
        yok <- NULL
    }
    x <- x[xok]
    if (paired) {
        x <- x - y
        y <- NULL
    }
    nx <- length(x)
    mx <- mean(x)
    vx <- var(x)
    if (is.null(y)) {
        if (nx < 2) 
            stop("not enough 'x' observations")
        df <- nx - 1
        stderr <- sqrt(vx/nx)
        if (stderr < 10 * .Machine$double.eps * abs(mx)) 
            stop("data are essentially constant")
        tstat <- (mx - mu)/stderr
        method <- if (paired) 
            "Paired t-test"
        else "One Sample t-test"
        estimate <- setNames(mx, if (paired) 
            "mean difference"
        else "mean of x")
    }
    else {
        ny <- length(y)
        if (nx < 1 || (!var.equal && nx < 2)) 
            stop("not enough 'x' observations")
        if (ny < 1 || (!var.equal && ny < 2)) 
            stop("not enough 'y' observations")
        if (var.equal && nx + ny < 3) 
            stop("not enough observations")
        my <- mean(y)
        vy <- var(y)
        method <- paste(if (!var.equal) 
            "Welch", "Two Sample t-test")
        estimate <- c(mx, my)
        names(estimate) <- c("mean of x", "mean of y")
        if (var.equal) {
            df <- nx + ny - 2
            v <- 0
            if (nx > 1) 
                v <- v + (nx - 1) * vx
            if (ny > 1) 
                v <- v + (ny - 1) * vy
            v <- v/df
            stderr <- sqrt(v * (1/nx + 1/ny))
        }
        else {
            stderrx <- sqrt(vx/nx)
            stderry <- sqrt(vy/ny)
            stderr <- sqrt(stderrx^2 + stderry^2)
            df <- stderr^4/(stderrx^4/(nx - 1) + stderry^4/(ny - 
                1))
        }
        if (stderr < 10 * .Machine$double.eps * max(abs(mx), 
            abs(my))) 
            stop("data are essentially constant")
        tstat <- (mx - my - mu)/stderr
    }
    if (alternative == "less") {
        pval <- pt(tstat, df)
        cint <- c(-Inf, tstat + qt(conf.level, df))
    }
    else if (alternative == "greater") {
        pval <- pt(tstat, df, lower.tail = FALSE)
        cint <- c(tstat - qt(conf.level, df), Inf)
    }
    else {
        pval <- 2 * pt(-abs(tstat), df)
        alpha <- 1 - conf.level
        cint <- qt(1 - alpha/2, df)
        cint <- tstat + c(-cint, cint)
    }
    cint <- mu + cint * stderr
    names(tstat) <- "t"
    names(df) <- "df"
    names(mu) <- if (paired) 
        "mean difference"
    else if (!is.null(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"
    rval
}
<bytecode: 0x000002a242c42160>
<environment: namespace:stats>

参考引用資料

  1. https://ja.wikipedia.org/wiki/%E3%82%A6%E3%82%A7%E3%83%AB%E3%83%81%E3%81%AEt%E6%A4%9C%E5%AE%9A
  2. https://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B8%E3%82%B9%E3%83%86%E3%82%A3%E3%83%83%E3%82%AF%E5%88%86%E5%B8%83

最終更新

Sys.time()
[1] "2024-04-06 08:19:05 JST"

R、Quarto、Package

R.Version()$version.string
[1] "R version 4.3.3 (2024-02-29 ucrt)"
quarto::quarto_version()
[1] '1.4.542'
packageVersion(pkg = "tidyverse")
[1] '2.0.0'

著者