t検定（1群） - 丸井綜研

「二項検定と二種類の過誤 - 丸井綜研」で、Common Lispで二項検定ができるようになりました。今日は1群のt検定を実装したいと思います。

t分布

t検定をするには、まずt分布の確率密度関数と累積分布関数を作る必要があります（検定だけなら累積分布関数だけでよいです）。確率密度関数の計算にはガンマ関数が、累積分布関数の計算には正則不完全ベータ関数が必要ですが、これらはすでに作ってあるのでそのまま使います。

Wikipedia日本語版の記事に書いてあった数式をそのままCommon Lispのコードにしたのが以下のものです。いくつかの場合でRでの計算結果と比べてみましたが、今のところ小数点以下6～7桁目くらいまでは一致しているようです。

;; PDF: Student's t Probability Density Function
;; (dt x df)
;; x: a value
;; df: degrees of freedom
;;
(defun dt (tval df)
  (* (/ (gamma (/ (+ df 1) 2)) (gamma (/ df 2)))
     (/ 1 (sqrt (* df pi)))
     (/ 1 (expt (+ 1 (/ (* tval tval) df)) (/ (+ df 1) 2)))))

;; CDF: Student's t Cumulative Distribution Function
;; (pt x df)
;; x: a value
;; df: degrees of freedom
;;
(defun pt (tval df)
  (ribeta (/ df 2.0) (/ df 2.0)
          (/ (+ tval (sqrt (+ (* tval tval) df)))
             (* 2 (sqrt (+ (* tval tval) df))))))

t検定のコード

お待ちかねのt検定です。これもWikipediaの数式を参考にしつつ、愚直に実装するだけです。出力としては両側検定のときのp値が欲しかったので、t値の値で場合分けをして計算してあります。

;;; one-sample t-test of mean with known mean (returns p-value of
;;; two-sided test).
;;;
;;; H0: mean of x = mu
;;;
(defun ttest1 (x &optional (mu 0))
  (let ((tvalue (/ (- (mean x) mu)
                   (sqrt (/ (vars x) (length x)))))
        (df (1- (length x))))
    (if (>= tvalue 0)
        (+ (- 1 (pt tvalue df)) (pt (- tvalue) df))
        (+ (pt tvalue df) (- 1 (pt (- tvalue) df))))))

使用例

例えば、とあるダイエット食品を使用した10名について、開始6か月後の体重減少ぶんを測定したとします（仮想データ）。それが有効なのかを検定してみましょう。

CL-USER> (setq x '(-3.5 -4.3 -6.4 -7.5 +4.0 -10.0 -2.4 +0.4 +0.9 -2.0))

CL-USER> (mean x)
-3.0800002

平均で3 kgほど体重減少しています。中には7.5 kgや10 kgも落とした人がいる反面、4.0 kg増えてしまっているケースもあります。帰無仮説を「体重変化は0.0 kgだ」とした上で、両側検定をしてみましょう。両側なので「体重変化があったか（増えたほう・減ったほうどちらであっても）」を見ることになります。有意水準αは0.05とし、それを下回ったら帰無仮説を棄却して、統計的有意な体重変化があったとします。

CL-USER> (ttest1 x)
0.04613601467151503d0

得られたp値が0.046なので、帰無仮説は棄却しましょう。……とは言え、統計的有意差よりも実質差のほうがはるかに大切です。実際に増減した体重がわずかであっても人数を集めれば有意差は出せてしまうので、検定の使い方は注意しないとですね。

Rでやると

同じことをRでやってみると以下のようになりました。

> x <- c(-3.5, -4.3, -6.4, -7.5, +4.0, -10.0, -2.4, +0.4, +0.9, -2.0)
> t.test(x)

    One Sample t-test

data:  x
t = -2.3113, df = 9, p-value = 0.04614
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -6.09447473 -0.06552527
sample estimates:
mean of x
    -3.08