僕が作った架空の例題で説明します。独立した２群の平均値には差があるのか・・・？
データプレパレーション
例題 回復期病棟で３週間入院した脳卒中患者を対象とした。A病院の患者群をA群、B病院の患者群をB群としてFIM利得を比較した。両群の患者属性、入院時FIM、治療内容については傾向スコアにて調整した。
同じ分散をもつ正規分布から抽出したサンプルを使用します
A <- round(rnorm(28, 25, 10),0)
B <- round(rnorm(31, 30, 10),0)
乱数なので出るたびに数値は異なりますので注意してください
ここでは下記のcsvファイルを使用します
作業ディレクトリを指定して取り込んでください
dat <- read.csv("t-test.csv", header=T, fileEncoding="Shift-JIS")
head(dat)
それぞれの群のサンプルサイズが異なり、また独立した群の比較なので、縦のデータに変換しましょう
#縦のデータセットに変更
dat2 <- tidyr::gather(dat, key=group, value=data)
#欠損値を削除
tate <- na.omit(dat2)
head(tate)
tidyr:: これはライブラリtidyrの関数を使ってます、という記号です。::がついてたらパッケージを起動「library(tidyr)」してください。
グループ毎のデータ要約
boxplot(data~group, data=tate)
psych::describeBy(tate$data, tate$group)
t検定
今回は同じ分散を設定した正規分布からの乱数ですので等分散を仮定（var.equal=T）した検定となります
t.test(data~group, var.equal=T, data=tate)
A群平均は25.11、B群平均は30.71、p値は0.022（＜0.05）となりました
有意水準5％であれば有意な差と言えます
必要なのが検定の結果だけであればここまででOKです！
＊これから先はデータサイエンス入門になります
合併した分散（併合分散, プールした分散）
対応のある2群の平均値は1標本の問題となるので、シンプルな分散でした
もし両群の分散（\(\sigma^2\)）が既知であったならば、分散の加法性よりZ値は以下のようになります
A群の平均を\(\bar{A}\)、B群の平均を\(\bar{B}\)
A群の母平均を\(\mu_A\)、B群の母平均を\(\mu_B\)
A群のサンプルサイズを\(n_A\)、B群のサンプルサイズを\(n_B\)
\(Z=\dfrac{(\bar{A} – \bar{B})-(\mu_A – \mu_B)}{\sqrt{\dfrac{\sigma^2}{n_A}+\dfrac{\sigma^2}{n_B}}} \, \sim N(0, \, 1)\)
ただし、\(\sigma^2\) は未知なので、それぞれの不偏分散\(s_A^2\)、\(s_B^2\)に置き換えなければなりません
t分布より
\(t=\dfrac{Z}{\sqrt{\dfrac{W}{k}}} \, \sim t(k)\)
この \(W\) は自由度 \(k\) の \(\chi^2\) 分布である
A, Bの母集団は正規分布であるため、カイ二乗分布の性質より以下のことが成り立ちます（ここがポイント）
\(\chi_A^2=\dfrac{\Sigma{(A_i-\bar{A})^2}}{\sigma^2}\,\sim\chi(n_A\,-\,1)\)
\(\chi_B^2=\dfrac{\Sigma{(B_i-\bar{B})^2}}{\sigma^2}\,\sim\chi(n_B\,-\,1)\)
\(w=\chi_A^2 \, + \, \chi_B^2\)とするとｶｲ二乗分布の再生性より
\(w = \dfrac{(n_A-1)s_A^2+(n_B-1)s_B^2}{\sigma^2} \, \sim\chi(n_A + n_B-2)\)
これで\(\sigma^2\) を不偏分散\(s_A^2\)、\(s_B^2\)に置き換えることができます
\(\sigma^2=\dfrac{(n_A-1)s_A^2+(n_B-1)s_B^2}{w}\)
これまでの結果から \(T\) は自由度\(n_A+n_B-2\)のｔ分布に従います
\(T=\dfrac{Z}{\sqrt{(\dfrac{w}{n_A+n_B-2})}}\)
\(Z\) と \(w\) を上式に当てはめると以下のようになります
\(T=\dfrac{(\bar{A} – \bar{B})-(\mu_A – \mu_B)}{\sqrt{\dfrac{\sigma^2}{n_A}+\dfrac{\sigma^2}{n_B}}}\sqrt{\dfrac{\sigma^2(n_A+n_B-2)}{(n_A-1)s_A^2+(n_B-1)s_B^2}}\)
\(=\dfrac{(\bar{A} – \bar{B})-(\mu_A – \mu_B)}{s^2\sqrt{\dfrac{1}{n_A}+\dfrac{1}{n_B}}}\)
\(s^2=\dfrac{(n_A-1)s_A^2+(n_B-1)s_B^2}{n_A+n_B-2}\)・・・これが合併した分散
これでｔ値を求めることができます。サンプルデータを当てはめてみましょう！
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))))
Rと全く同じt値を求める場合は・・・
sA <- sd(tate[tate$group=="A","data"])
sB <- sd(tate[tate$group=="B","data"])
mA <- mean(tate[tate$group=="A","data"])
mB <- mean(tate[tate$group=="B","data"])
s2 <- ((28-1)*sA^2 + (31-1)*sB^2)/(28+31-2)
(t <- (mA-mB)/(sqrt(s2*(1/28+1/31))))
検定方法
t値による検定方法は対応ありをご確認ください
コメント