僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。
例題）回復期病棟で３週間入院した脳卒中患者を対象とした。X病院の患者群をX群、Y病院の患者群をY群としてFIM利得を比較した。両群の患者属性、入院時FIM、治療内容については傾向スコアにて調整した
データの準備と要約
使用するパッケージ（パッケージのダウンロード）
library(tidyr)
library(ggplot2)
同じ分散をもつ正規分布から抽出したサンプル t-test を使用します
X <- round(rnorm(28, 25, 10),0)
Y <- round(rnorm(31, 30, 10),0)
乱数なので出るたびに数値は異なりますので注意してください
ここでは下記のcsvファイルを使用して学習します
データセットt-testをdatに格納します（ファイルの読み込み）
dat <- read.csv("t-test.csv", header=T, fileEncoding="Shift-JIS")
head(dat)
それぞれの群のサンプルサイズが異なり、また独立した群の比較なので、縦のデータに変換しましょう（縦のデータと横のデータ）
dat_long <- gather(dat, key = "group", value = "data", na.rm = TRUE)
# 箱ひげ図を描画
ggplot(dat_long, aes(x = group, y = data)) +
geom_boxplot() +
theme_minimal()
tidyr:: これはでパッケージtidyrを起動できます。つまり「::」を使用することで、library関数の記述は不要になります。
各グループの記述統計
psych::describeBy(dat_long$data, dat_long$group)
t検定
今回は同じ分散を設定した正規分布からの乱数ですので等分散を仮定（var.equal=TRUE）した検定となります
t.test(data~group, var.equal=TRUE, data=dat_long)
A群平均は25.11、B群平均は30.71、p値は0.022（＜0.05）となりました
有意水準5％であれば有意な差と言えます
必要なのが検定の結果だけであればここまででOKです。ここからは、もう少し詳しく見ていきましょう。
合併した分散（併合分散, プールした分散）
２群の分散が同じだと仮定した場合（等分散）
対応のある2群の平均値は1標本の問題となるので、シンプルな分散でした（対応のあるｔ検定）
もし両群の分散（\(\sigma^2\)）が既知であったならば、分散の加法性よりZ値は以下のようになります
X群の平均を\(\bar{X}\)、Y群の平均を\(\bar{Y}\)
X群の母平均を\(\mu_x\)、Y群の母平均を\(\mu_y\)
X群のサンプルサイズを\(n_x\)、Y群のサンプルサイズを\(n_y\)
\(Z=\dfrac{(\bar{X} – \bar{Y})-(\mu_x – \mu_y)}{\sqrt{\dfrac{\sigma^2}{n_x}+\dfrac{\sigma^2}{n_y}}} \, \sim N(0, \, 1)\)
ただし、\(\sigma^2\) は未知なので、それぞれの不偏分散 \(\hat{s_x}^2\)、\(\hat{s_y}^2\)に置き換えなければなりません
t分布より
\(t=\dfrac{Z}{\sqrt{\dfrac{W}{k}}} \, \sim t(k)\)
この \(W\) は自由度 \(k\) の \(\chi^2\) 分布です
X, Yの母集団は正規分布であるため、カイ二乗分布の性質より以下のことが成り立ちます（ここがポイント）
\(\chi_X^2=\dfrac{\Sigma{(X_i-\bar{X})^2}}{\sigma^2}\,\sim\chi(n_X\,-\,1)\)
\(\chi_Y^2=\dfrac{\Sigma{(Y_i-\bar{Y})^2}}{\sigma^2}\,\sim\chi(n_Y\,-\,1)\)
\(w=\chi_X^2 \, + \, \chi_Y^2\)とするとｶｲ二乗分布の再生性より
\(w = \dfrac{(n_x-1)s_X^2+(n_y-1)s_Y^2}{\sigma^2} \, \sim\chi(n_x + n_y-2)\)
これで\(\sigma^2\) を不偏分散$\hat{s_x}^2、\hat{s_y}^2$に置き換えることができます
\(\sigma^2=\dfrac{(n_x-1)\hat{s_x}^2+(n_y-1)\hat{s_y}^2}{w}\)
これまでの結果から \(T\) は自由度\(n_x+n_y-2\)のｔ分布に従います
\(T=\dfrac{Z}{\sqrt{(\dfrac{w}{n_x+n_y-2})}}\)
\(Z\) と \(w\) を上式に当てはめると以下のようになります
\(T=\dfrac{(\bar{X} – \bar{Y})-(\mu_x – \mu_y)}{\sqrt{\dfrac{\sigma^2}{n_x}+\dfrac{\sigma^2}{n_y}}}\sqrt{\dfrac{\sigma^2(n_x+n_y-2)}{(n_x-1)\hat{s_x}^2+(n_y-1)\hat{s_y}^2}}\)
\(=\dfrac{(\bar{X} – \bar{Y})-(\mu_x – \mu_y)}{s^2\sqrt{\dfrac{1}{n_x}+\dfrac{1}{n_y}}}\)
\(s^2=\dfrac{(n_x-1)\hat{s_x}^2+(n_y-1)\hat{s_x}^2}{n_x+n_y-2}\)・・・これが合併した分散
これでｔ値を求めることができます。サンプルデータを当てはめてみましょう！
t値
s2 <- ((28-1)*9.7^2 + (31-1)*8.57^2)/(28+31-2)
(t <- (25.1-30.7)/(sqrt(s2*(1/28+1/31))))
データセットから標準偏差と平均を求める方法
sx <- sd(dat_long[dat_long$group=="X", "data"])
sy <- sd(dat_long[dat_long$group=="Y","data"])
mx <- mean(dat_long[dat_long$group=="X", "data"])
my <- mean(dat_long[dat_long$group=="Y", "data"])
s2 <- ((28 - 1)*sx^2 + (31 - 1)*sy^2)/(28 + 31 - 2)
(t <- (mx - my)/(sqrt(s2*(1/28 + 1/31))))
２群の分散が異なると仮定した場合は、Welchの検定の適用になります
RでWelchの検定を実行
t.test(data~group, var.equal=FALSE, data=dat_long)
◆ ｔ値＝-2.3408の求め方
$\hat{s_X}^2$ と $\hat{s_Y}^2$ が等しいと想定できな場合は、次のようにｔ統計量を定義します
$T=\dfrac{\bar{x}-\bar{y}}{\sqrt{\dfrac{\hat{s_x}^2}{n_x}+\dfrac{\hat{s_x}^2}{n_y}}}$
自由度は以下のようになります（整数になるとは限りません）
$d_1=\dfrac{\hat{s_x}^2}{n_x}$、$d_1=\dfrac{\hat{s_y}^2}{n_y}$
$df = \dfrac{(d_1+d_2)^2}{\dfrac{d_1^2}{n_x-1}+\dfrac{d_2^2}{n_y-1}}$
◆ 実際の数値を代入してみましょう
自由度
d1 <- 9.696656^2/28
d2 <- 8.572023^2/31
df <- (d1+d2)^2/((d1^2/(28-1))+(d2^2/(31-1)))
print(df)
ｔ値
t <- (25.10714-30.70968)/sqrt(d1+d2)
print(t)
p値
pt(-2.340832, 54.24476)*2
検定方法
結果の解釈については対応のあるｔ検定をご確認ください
コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください