R実装と解説対応のない2標本の母平均の差の検定(母分散が異なる)

母分散が異なるの場合の対応のない2標本の母平均の差の検定とは、2つの母集団が正規分布に従い、ともに母分散が異なるとき、一方の母平均が他方の母平均と「異なる」または「大きい」、「小さい」かどうかを、検定統計量がt分布に従うことを利用して検定します。

統計的検定の流れ

検定の大まかな流れを確認しておきます。

帰無仮説H₀と対立仮設H₁をたてます
有意水準を決めて、棄却域を決定します
検定統計量を計算します
検定統計量が棄却域にあれば帰無仮説H₀を棄却し、対立仮設H₁を採択します。そうでなければ、帰無仮説H₀を採択します

母集団Aの母平均と母集団Bの母平均が異なることを検定したい場合

母集団Aの母平均を $\mu_1$ 、母集団Bの母平均を $\mu_2$ とします。
母集団Aの母平均が母集団Bの母平均と異なることを検定したい場合は、両側検定として次のような仮説をたてます。

帰無仮説 $H_0: \mu_1 = \mu_2$
対立仮説 $H_1: \mu_1 \neq \mu_2$

母集団Aの母平均と母集団Bの母平均より大きいことを検定したい場合

母集団Aの母平均を $\mu_1$ 、母集団Bの母平均を $\mu_2$ とします。
母集団Aの母平均が母集団Bの母平均より大きいことを検定したい場合は、片側検定として次のような仮説をたてます。

帰無仮説 $H_0: \mu_1 = \mu_2$
対立仮説 $H_1: \mu_1 > \mu_2$

母集団Aの母平均と母集団Bの母平均より小さいことを検定したい場合

母集団Aの母平均を $\mu_1$ 、母集団Bの母平均を $\mu_2$ とします。
母集団Aの母平均が母集団Bの母平均より小さいことを検定したい場合は、片側検定として次のような仮説をたてます。

帰無仮説 $H_0: \mu_1 = \mu_2$
対立仮説 $H_1: \mu_1 < \mu_2$

理論

標本xの標本平均を $\overline{x}_1$ 、標本の大きさを $n_1$ 、母集団の分布を正規分布 $N(\mu_1, \sigma_1^2)$ 、分散の推定値を $\hat{\sigma}_1^2$ とします。
同様に、標本yの標本平均を $\overline{x}_2$ 、標本の大きさを $n_2$ 、母集団の分布を正規分布 $N(\mu_2, \sigma_2^2)$ 、分散の推定値を $\hat{\sigma}_2^2$ とします。

検定統計量

母分散が異なる場合の対応のない2標本の母平均の差の検定における検定統計量 $t$ は、次式で表され、自由度 $\nu$ のt分布に従います。

$t = \frac{\overline{x}_1-\overline{x}_2-(\mu_1-\mu_2)}{\sqrt{\frac{\hat{\sigma}_1^2}{n_1}+\frac{\hat{\sigma}_2^2}{n_2}}}$

$\nu \approx \frac{(\frac{\hat{\sigma}_1^2}{n_1}+\frac{\hat{\sigma}_2^2}{n_2})^2}{\frac{\hat{\sigma}_1^4}{n_1^2(n_1-1)}+\frac{\hat{\sigma}_2^4}{n_2^2(n_2-1)}}$

信頼区間

信頼区間は次のように表すことができます。

$P(-t_{\frac{\alpha}{2}} \leq \frac{\overline{x}_1-\overline{x}_2-(\mu_1-\mu_2)}{\sqrt{\frac{\hat{\sigma}_1^2}{n_1}+\frac{\hat{\sigma}_2^2}{n_2}}} \leq t_{\frac{\alpha}{2}}) = 1-\alpha$

これを変形すると次のようになります。

$P(\overline{x}_1-\overline{x}_2 - t_{\frac{\alpha}{2}}\sqrt{\frac{\hat{\sigma}_1^2}{n_1}+\frac{\hat{\sigma}_2^2}{n_2}} \leq \mu_1-\mu_2 \leq \overline{x}_1-\overline{x}_2 + t_{\frac{\alpha}{2}}\sqrt{\frac{\hat{\sigma}_1^2}{n_1}+\frac{\hat{\sigma}_2^2}{n_2}}) = 1-\alpha$

よって、信頼係数 $1-\alpha$ の信頼区間は次のようになります。

$[\overline{x}_1-\overline{x}_2 - t_{\frac{\alpha}{2}}\sqrt{\frac{\hat{\sigma}_1^2}{n_1}+\frac{\hat{\sigma}_2^2}{n_2}}, \overline{x}_1-\overline{x}_2 + t_{\frac{\alpha}{2}}\sqrt{\frac{\hat{\sigma}_1^2}{n_1}+\frac{\hat{\sigma}_2^2}{n_2}}]$

R実装

Rで母分散が異なるの場合の対応のない2標本の母平均の差の検定を行う場合は、基本パッケージstatsに実装されているt.test関数を用います。
getAnywhere関数の引数にt.test.defualtを渡すと


getAnywhere(t.test.default)

t.test.defualt関数の実装を見ることができます。


A single object matching ‘t.test.default’ was found
It was found in the following places
  registered S3 method for t.test from namespace stats
  namespace:stats
with value

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(deparse(substitute(x)), "and", deparse(substitute(y)))
        if (paired) 
            xok <- yok <- complete.cases(x, y)
        else {
            yok <- !is.na(y)
            xok <- !is.na(x)
        }
        y <- y[yok]
    }
    else {
        dname <- deparse(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 of the differences"
        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 || !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, 
        alternative = alternative, method = method, data.name = dname)
    class(rval) <- "htest"
    return(rval)
}
<bytecode: 0x69251f8>
<environment: namespace:stats>

検定統計量

母分散が異なる場合の対応のない2標本の母平均の差の検定における検定統計量を計算するR実装のうち、本質的な箇所を抜き出すと以下のようになります。
自由度の計算が少しわかりにくいですが、次のようになっております。

$\nu \approx \frac{(\frac{\hat{\sigma}_1^2}{n_1}+\frac{\hat{\sigma}_2^2}{n_2})^2}{\frac{\hat{\sigma}_1^4}{n_1^2(n_1-1)}+\frac{\hat{\sigma}_2^4}{n_2^2(n_2-1)}} = \frac{\Biggl(\sqrt{\biggl(\sqrt{\frac{\hat{\sigma_1^2}}{n_1}}\biggr)^2+\biggl(\sqrt{\frac{\hat{\sigma_2^2}}{n_2}}\biggr)^2}\Biggr)^4}{\biggl(\sqrt{\frac{\hat{\sigma_1^2}}{n_1}}\biggr)^4\frac{1}{n_1-1}+\biggl(\sqrt{\frac{\hat{\sigma_2^2}}{n_2}}\biggr)^4\frac{1}{n_2-1}} = \frac{(\sqrt{\mbox{stderrx}^2+\mbox{stderry}^2})^4}{\frac{\mbox{stderrx}^4}{\mbox{nx}-1}+\frac{\mbox{stderry}^4}{\mbox{ny}-1}}$


nx <- length(x)                                                  # xの標本の大きさ
mx <- mean(x)                                                    # xの標本平均
vx <- var(x)                                                     # xの分散の推定値(不偏分散)
if (is.null(y)) {
}
else {
    ny <- length(y)                                              # yの標本の大きさ
    my <- mean(y)                                                # yの標本平均
    vy <- var(y)                                                 # yの分散の推定値(不偏分散)
    if (var.equal) {
    }
    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))
    }
    tstat <- (mx - my - mu)/stderr                              # 検定統計量
}

p値と信頼区間

p値と信頼区間に関するR実装のうち、本質的な箇所を抜き出すと以下のようになります。
pt(t₀, df)関数は、自由度dfのt分布における統計量t₀に対する下側確率P(t < t₀)を計算します。
qt(p, df)関数は、自由度dfのt分布における確率pに対する統計量t₀を計算します。


if (alternative == "less") {
    pval <- pt(tstat, df)                       # 対立仮説x<y: p値
    cint <- c(-Inf, tstat + qt(conf.level, df)) # 対立仮説x<y: 信頼係数conf.levelの限界値
}
else if (alternative == "greater") {
    pval <- pt(tstat, df, lower.tail = FALSE)   # 対立仮説x>y: p値, lower.tail=FALSEで上側確率を計算
    cint <- c(tstat - qt(conf.level, df), Inf)  # 対立仮説x>y: 信頼係数conf.levelの限界値
}
else {
    pval <- 2 * pt(-abs(tstat), df)             # 対立仮説x!=y: p値, 両側検定なので下側確率の2倍を計算
    alpha <- 1 - conf.level                     # 有意水準
    cint <- qt(1 - alpha/2, df)                 # 両側検定なので、1から有意水準の1/2を引いた確率に対する統計量を計算
    cint <- tstat + c(-cint, cint)              # 対立仮説x!=y: 信頼係数conf.levelの限界値
}
cint <- mu + cint * stderr                      # 信頼区間

例

あるクラスAとBの身長を計測したところ、次のようになりました。
身長の分布が正規分布をすると仮定すると、クラスAとBの平均身長は差があるかを検定してみます。
ただし、有意水準5%とします。

クラスA：166, 168, 168, 169, 170, 171, 172, 173, 174, 175
クラスB：166, 166, 167, 167, 168, 168, 168


t.test(c(166, 168, 168, 169, 170, 171, 172, 173, 174, 175),
         c(166, 166, 167, 167, 168, 168, 168),
         alternative = "two.sided", mu = 0, conf.level = 0.95)

    Welch Two Sample t-test

data:  c(166, 168, 168, 169, 170, 171, 172, 173, 174, 175) and c(166, 166, 167, 167, 168, 168, 168)
t = 3.5201, df = 11.305, p-value = 0.004607
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.302622 5.611664
sample estimates:
mean of x mean of y 
 170.6000  167.1429

p値0.004607 < 有意水準0.05なので、棄却。
よって、クラスAとBの平均身長は差があるといえる。

R実装と解説対応のない2標本の母平均の差の検定(母分散が異なる)

統計的検定の流れ

母集団Aの母平均と母集団Bの母平均が異なることを検定したい場合

母集団Aの母平均と母集団Bの母平均より大きいことを検定したい場合

母集団Aの母平均と母集団Bの母平均より小さいことを検定したい場合

理論

検定統計量

信頼区間

R実装

検定統計量

p値と信頼区間

例

こちらもおすすめ

R 関数の引数が省略記号（…）を持つかどうかを判定する方法

R K-means法のクラスタ数を機械的に決定する方法

R×qgraph ネットワーク・グラフ一覧