東京国際大学大学院

心理データ解析

東京大学大学院農学生命科学研究科 大森宏


R の基本的な使い方を学ぶ

1 + 1      #数値演算 
x <- c(2, 3, 5)      #ベクトルの定義 
y <- c(10, 14, -7)         
x + y      #ベクトルの和 
2x      #ベクトルのスカラー倍 
sum(x)      #ベクトル要素の和 
1:10      #連続した値 
x <- 1:10      # 
y <- c(1,3,2.8,2.5,3,3.5,3.8,4,3.5,3.6)      # 
length(y)      # y の長さ 
plot(x, y)      #グラフ表示 
plot(x, y, type="l")      #グラフ表示 

データ分布と基本統計量

データ

データ:母集団の特性を調べるために抽出した標本の観測値.標本は通常,無作為標本(ランダムサンプル)
大きく,質的データと量的データに分けられる.

量的データ

長さ,身長,重さ,売り上げ,得点など.
データ(観測値):x1x2,…,xn, サンプルサイズ(データ数)は n.
統計量:データから得られる数値.サンプルサイズ,平均,最大値…など.

データの視覚化

位置情報

データの中心的な位置を示す統計量

ちらばりの情報

データのばらつきの程度を示す統計量

英語得点データ

英語得点データを用いて,データの基本統計量の計算演習を行う.
#英語の得点
eigo <- c(
36,70,56,68,76,60,50,63,62,42,64,60,50,68,71,67,
50,65,67,57,72,64,61,66,46,80,46,51,59,32,55,65,
65,52,57,64,23,57,53,54,38,71,57,69,77,61,51,64,
63,43,65,61,51,69,72,68,53,66,68,58,73,65,62,67,
47,81,47,52,59,33,56,66,67,52,58,65,24,58,54,55)
length(eigo)      #データ数 
mean(eigo)      #標本平均 
var(eigo)      #標本分散 
sd(eigo)      #標本標準偏差 
boxplot(eigo, main="英語得点の箱ヒゲ図")      #箱ひげ図 
boxplot.stats(eigo)      #箱ひげ図用統計量 
hist(eigo, breaks=seq(0, 100, by=5), xlab="英語得点", ylab="頻度", main="")
title(main = "英語得点のヒストグラム")      #グラフタイトル 
stem(eigo, scale=2)      #幹葉表示 
#国語の得点
kokugo <- c(
68,56,54,15,39,30,78,56,45,55,41,44,89,57,34,60,66,
40,62,76,67,70,68,53,46,82,61,47,56,75,52,59,21,73,
76,25,63,54,81,57,38,62,70,53,47,29,66,61,23,72,41,
65,58,54,40,49,39,50,36,52,46,71,34,64,50)
 
boxplot(eigo,kokugo, xaxt="n")
axis(1, 1:2, labels=c("英語","国語"))
title(main = "英語と国語得点の箱ひげ図")
幹葉表示の R の出力
> stem(eigo,scale=2)

The decimal point is 1 digit(s) to the right of the |     

2 | 34
2 |
3 | 23
3 | 68
4 | 23
4 | 6677
5 | 0001112223344
5 | 5566777788899
6 | 0011122334444
6 | 5555556667777888899
7 | 011223
7 | 67
8 | 01

eigobox eigobox

データ操作

課題

  1. 国語データ得点を用いて,英語データと同様の解析をせよ.
  2. A君の得点は英語,国語とも65点であった.英語と国語,どちらの方が成績が良かった と言えるか.

確率分布

離散確率分布

 ある量の集まり P = { p1,…,pn } の中で,
fig1
という性質をもつもの.
 また,離散的な変数 X = { x1,…,xn } のおのおのの値に対し,それが 生起する確率 pi が与えられているとき,X を離散確率変数という. これは,

表 1 :離散確率変数(n = 5)
 変数 X     x1   x2   x3   x4   x5 
 確率 P     p1  p2  p3  p4  p5

と表せる.

連続型確率分布

 関数 f(x) が,
fig2
という性質を持つとき,これを連続型確率分布という.

累積分布関数

Pr[A] を事象 A が生起する確率とする.連続型確率変数 X に対し,X が x 以下である確率,
fig3
で定義される関数 F(x) を累積分布関数という. この関数を用いると,確率変数 X が区間 (a, b) に落ちる確率は,
fig4
で表せる.なお,離散型確率変数でも積分を和に変えることにより,同様に階段型の累積分布関数が定義できる.

分布の代表値

平均

 平均(mean)μ は,分布の中心的な位置座標を表す.確率変数 X に対しては, X の期待値 E[X ] とも表記し,
fig5
と定義される.離散確率変数では,積分を総和に変えることにより括弧内のように定義できる.

分散

 分散(variance) σ2 は,分布の拡がりの程度を表す. 確率変数 X に対しては,Var[X ] と表記する. 確率変数 X の関数 (X - μ)2 の期待値でもあり,
fig6
と定義される.
 また,σ を標準偏差(SD, standard deviation)といい,平均と同じ次元で分布の拡がりの大きさ を表す量である.

1次元の分布

離散分布の例

離散一様分布

サイコロの目の分布

離散一様分布
 変数 X     1   2   3   4   5   6 
 確率 P     1/6  1/6  1/6  1/6  1/6  1/6

# 離散一様分布の平均と分散の計算
x <- 1:6      # 1 から 6 まで並べたベクトル 
p <- rep(1/6, 6)      # 1/6 を6個並べたベクトル 
m <- sum(x*p)      #平均 
s2 <- sum(p*(x-m)^2)      #分散 

二項分布

 成功確率 p の事象を n 回繰り返したときの成功回数 x の分布
fig13
平均: E[X ] = np,分散: Var[X ] = npq
X 〜 B(n,p) と書くこともある.

ポアソン(Poisson)分布

 正のパラメータ λ と,0 以上の整数 X に対し,確率密度が
fig16
である分布
平均: E[X ] = λ,分散: Var[X ] = λ

連続型分布の例

一様分布(uniform distribution)

 2 つのパラメータ a,b(a<b)をもつ確率密度関数が
dunif
で表される分布.

平均: E[X ] = (a + b)/2,分散: Var[X ] = (b - a)2/12

課題:π の近似の精度を上げて π の近似値を再計算せよ.

正規分布(normal distribution)

 平均 μ,分散 σ2 の2つのパラメータをもつ確率密度関数が
dnorm
で表される分布で,N(μ,σ2)と表記する.μ は位置パラメータで, スケールパラメータ σ を標準偏差(standard deviation)という.

χ2(カイ2乗)分布

 正の自由度パラメータ n をもち,正の値しか持たない分布.
平均:E[X ] = n,分散:Var[X ] = 2n.

F 分布

 正の2つの自由度パラメータ m,n をもち,正の値しか取らない分布

t 分布

 正の自由度パラメータ n をもつ分布. 標準正規分布より裾が重く(x が 0 より離れてもなかなか確率密度が 0 に近づかない), 自由度が小さいほど裾が重くなる. 自由度 n = 1 のときはコーシー分布になり,n = ∞ のときは標準正規分布となる.

正規母集団からの標本に基づく推論

独立な正規分布の合成分布

平均 μ1, 分散 σ12,の正規分布 からの標本 x 〜 N( μ1,σ12 ) と, 平均 μ2, 分散 σ22,の正規分布 からの標本 y 〜 N( μ2,σ22 ) があり,両者が互いに独立であるとする.(y の値は x の値の影響を受けない.)

正規分布に基づく母数の区間推定

 正規分布は,平均 μ と分散 σ2 の2つの母数を持つ. 2つの母数とも未知であるのが普通であるが,片方が既知であるときは母数に関する推論は 簡単に行える.このため,多少非現実的な設定であるが,まず,既知の場合を考え,その後,より 一般的である2つの母数とも未知である場合を扱う.

分散既知の場合の母平均 μ の区間推定

 正規分布する母集団で母分散 σ2 がわかっている場合は,未知の母平均 μ に関する区間推定は以下のように行える.
いま,正規分布 N( μ,σ2 ) において,大きさ n の標本 x1x2,…,xn を抽出したとき,母平均は標本平均で推定される.標本平均 x- の分布は,
meanCI
となる.標準正規分布の 97.5%分位点を z0.975(= 1.96)とすると, 標準正規分布する確率変数 z が -z0.975 から z0.975 に入る 確率は 0.95 となる.つまり,
meanCI
となる.最後の式を母集団平均 μ の 95% 信頼区間と言う.
このように,母数の信頼区間を標本から推定することを区間推定という.区間推定においては, 信頼区間の幅 2d が小さい程よい.すなわち,母分散が小さい母集団で, 標本の大きさ(サンプルサイズ)が大きい程,精度の高い推定が行える.

95% の意味

 同じ正規母集団から標本抽出を繰り返すと,毎回標本平均として異なる値がえられ,それに 対応して信頼区間も異なる.この信頼区間の 95% が真の平均 μ を含む,という意味である.
つまり,100回の標本抽出により,100 個の信頼区間を作ったら平均的にみて,95 個の信頼区間が 真の平均 μ を含むことが期待できる.
下の図は,平均 0 分散 2 の正規分布 N( 0, 2 ) から大きさ 10 の標本を取りだし,分散が既知であるとして, 母平均に対する信頼区間を 100 個生成したものである."×" が標本平均を示す.左の "*" は,信頼区間 が母平均の真値 0 を含まなかった場合である.
normsamp
# 平均 μ の 95 %信頼区間 100 回生成の R スクリプト
v <- 2; n <- 10 # 母分散と標本の大きさ
d <- qnorm(0.975)*sqrt(v/n) # 95%信頼区間の幅  
x <- c(-2.5, -2.5, 2.5, 2.5) # グラフの x の範囲  
y <- c(0, 100, 100, 0) # グラフの y の範囲  
plot(x, y, type="n", xlab="", ylab="") # グラフ領域確保  
segments(0, 0, 0, 100, col="red") # 母平均のライン  
for(i in 0:100){ # 100回の繰り返し  
m <- mean(rnorm(n, mean=0, sd=sqrt(v))) # N(0, v)からの大きさ n の標本平均 
segments(m-d, i, m+d, i) # 信頼区間の表示 
points(m, i, pch=4, col="red", cex=0.8) # 標本平均の赤×表示  
if(m-d>0 || m+d<0) text(-2.5, i, "*") # 信頼区間が母平均を含まない(失敗)した場合  
} #  
title(main="N( 0, 2 )からの大きさ 10 の標本から得られた  
平均 μ の 95 %信頼区間を 100 回作成", cex.main=1.0) #  

平均既知の場合の母分散 σ2 の区間推定

 正規母集団で母平均 μ がわかっているとき,大きさ n の標本 x1x2,…,xn を抽出したとき,母分散は,
varCI
で推定される.ところで,標本は
varCI
と分布するので,自由度 n の χ2 の 2.5%分位点と 97.5%分位点をそれぞれ, χ2(n)0.025,χ2(n)0.975 とすると,
varCI
が成り立つ.下の式の区間を母分散 σ2 の 95% 信頼区間と言う.

平均未知の場合の母分散 σ2 の区間推定

 正規母集団では,母数が未知であるのが普通であろう.このとき, 大きさ n の標本 x1x2,…,xn を抽出したとき,母平均 μ と母分散 σ2 は,そえぞれ 標本平均 x- と標本分散 s2
sampmv
で推定される.母平均 μ の信頼区間を述べる前に母分散 σ2 の 信頼区間の構成法を述べる.
 ところで,標本や標本平均は,
sampmv
と分布する.一方,
sampmv
と計算されるので,(n - 1)s22 という量は,
sampmv
と,自由度 n - 1 の χ2 分布,χ2(n - 1),に従うことがわかる.
 自由度 n - 1 の χ2 分布の 2.5%分位点と 97.5%分位点をそれぞれ, χ2(n - 1)0.025,χ2(n - 1)0.975 とすると,
varCI
が成り立つ.下の式の区間を母分散 σ2 の 95% 信頼区間と言う.

分散未知の場合の母平均 μ の区間推定

 前節で考えたように正規母集団の母数が未知のときは, 大きさ n の無作為標本 x1x2,…,xn から,母平均 μ と母分散 σ2 は,それぞれ 標本平均 x- と標本分散 s2
sampmv
で推定される.
 標本平均 x- の分布は標準化すると,
tCI
のように標準正規分布となり,標本分散に関係する量は,
tCI
のように自由度 n - 1 の χ2 分布する.これより,zU をその自由度 n - 1 で 割った量の平方根との比は,
tCI
のように自由度 n - 1 の t 分布に従う.
 自由度 n - 1 の t 分布の 97.5%分位点を t(n - 1)0.975 とすると, t 分布する確率変数 t が -t(n - 1)0.975 から t(n - 1)0.975 に入る 確率は 0.95 となる.つまり,
tCI
となる.最後の式を母分散未知のときの母集団平均 μ の 95% 信頼区間と言う.

仮説検定

帰無仮説(H0)と対立仮説(H1

統計学で扱う仮説とは,母集団に対する断定や推測.たとえば,

などである.

統計的仮説検定で用いられる仮説は,まず,帰無仮説という形式で与えられる.
帰無仮説棄却されることに意味がある仮説である.
帰無仮説と反対の仮説を対立仮説という.

上の3番目の例でみると,

帰無仮説: 母集団 A と母集団 B の平均は等しい. (H0: μA = μB
対立仮説: 母集団 A と母集団 B の平均は等しくない. (H0: μA ≠ μB

母集団 A と母集団 B は異なる処理(薬の投与など)をしているので,実験の目的 は,母集団 A と母集団 B の平均は異なる(処理効果がある)ことを言いたい (対立仮説が正しいことを望む)のだが,まずは 「等しい(処理効果無し)」と仮定してみようという考え方.

仮説検定

検定統計量

 標本から算出される量で,検定に用いられるもので,t 値,F 値などがある. この値から帰無仮説を受託(採択)するか 棄却(対立仮説の採択)するかを判定する.

有意水準

 統計的仮説検定では,たとえば2つの母集団平均が等しいという帰無仮説を考えると, この帰無仮説のもとで,検定統計量(標本平均の差に基づく t 値など)以上 (もしくは未満)の値が得られる確率を求める.
くだけた言い方をすれば,帰無仮説が正しいとしたときに,標本のようなデータが得られる確率 を求める.
これが十分小さい(ほとんどありえない)ときは,平均が等しいと仮定したことが誤りであったと判断して 帰無仮説を棄却し,2つの母集団平均には差があると結論づける.
この確率がそれほど小さく ない場合は,このような統計量が得られることもありえると考え,帰無仮説を採択し,平均が等しいと考え てもよいとする.
棄却か採択かの判断の基準となる確率を有意水準といい, 5 %1 % がよく用いられる.

片側検定と両側検定

実験状況によっては,薬投与などの処理を行った集団(処理群)平均 μA が,薬を投与しない 集団(対照群)の平均 μB より小さくなることはないことが事前に わかっているような場合が ある.このようなとき,

帰無仮説,H0: μA = μB
対立仮説,H0: μA > μB
となる.これは,事前情報より,μA < μB となる可能性 をまったく考えない場合である.
このため検定には,片側 5 %点や 1 %点を用いる.

両側検定と信頼区間

母集団平均に対する両側検定は,母集団平均に対する信頼区間と大きな関係がある. いま,帰無仮説(H0)と対立仮説(H1)が,

H0:μ = μ0
H1:μ ≠ μ0
であり,母分散 σ2 が既知のときを考える.
標本の大きさがnで,標本平均が x- であったとすると,母平均μに対する 95%信頼区間は,
Pr[ − 1.96 < √n(x- − μ )/σ < 1.96 ] = 0.95,
Pr[ x- − 1.96×σ/ √n < μ < x- + 1.96×σ/ √n ] = 0.95
となる.
 一方,この検定の検定統計量は,標本平均の標準化値の絶対値
|z| = √n|x- − μ0 |/σ
で, 有意水準 5 %で帰無仮説を受諾するのは,検定統計量 |z| が両側 5 %点である 1.96 以下のときである.つまり,
帰無仮説を受諾 ⇔ − 1.96 < √n(x- − μ0 )/σ < 1.96
である.この両者の関係より,
帰無仮説を受諾 ⇔ 母平均の信頼区間に μ0 が含まれる.
帰無仮説を棄却 ⇔ 母平均の信頼区間に μ0 が含まれない.
が成り立つ

検定における2種類の過誤

検定は,仮説を棄却するか採択するかのいずれかであるが, 統計量は分布をもつので,この判定には間違いが起こることがある.
以下のように,この過誤には 2 種類がある.

統計的検定における2種類の過誤
  仮説の棄却 仮説の採択
仮説が真のとき 第1種の過誤 正解
仮説が偽のとき 正解 第2種の過誤

第1種の過誤が有意水準である.また,第2種の過誤の確率を β としたとき, 仮説が偽のとき正しく仮説を棄却する確率,1 - β,を検出力という. よい検定は,第1種の過誤を固定したもとで検出力の高い検定方式である.

正規母集団の母平均に対する t 検定

1つの母集団に対する検定

 平均 μ,分散 σ2 がともに未知である正規母集団に対して,
帰無仮説 H0: μ = μ0
対立仮説 H1: μ ≠ μ0
の両側検定を考える.  いま,母集団から大きさ n の無作為標本 x1x2,…,xn を抽出したところ,標本平均が x-,標本分散が s2 であったとする. 帰無仮説のもとでは,標本平均は,
ttest
と分布するので,これを,標本平均の標準誤差 s/√n で標準化した t は,
ttest
のように自由度 n - 1 の t 分布に従う.この分布の97.5%分位点を t(n - 1)0.975 とすると, 有意水準 5 %の検定は,
|t| > t(n - 1)0.975
のとき帰無仮説を棄却する.|t| が検定統計量で,この値を |t| 値という.

なお,この検定は,対のある標本に適用できる.対のある標本とは,n 組のペアー標本,

(x1y1 ),(x2y2 ), …,(xnyn )
からなる.正規性の仮定のもとでは,
xi 〜 N( μi,σx2 ), yi 〜 N( μi + δ,σy2 )
ここで興味ある母数は δ であり,μ1,…,μn は攪乱母数 である.yixi の差を取ると,
ziyixizi 〜 N( δ,σz2 )
となるので,1つの母集団に対する検定に帰着する.なおこの問題は, 反復のない 2×n の2元配置と考えて解くこともできる.

2つの母集団に対する検定

 2つの母集団 A,B があり,それぞれが平均を μA,μB, 分散を σA2,σB2 の正規分布に従って いるが,その値は未知であるとする.いま,両集団の分散の値が等しく, σA2=σB2=σ2,と仮定 できるとしよう.このとき,

帰無仮説,H0: μA = μB
対立仮説,H1: μA ≠ μB

の検定は t 分布を用いて行える.

母集団 A から大きさ nA,母集団 B から大きさ nB の標本を抽出した. 母集団 A からの標本の標本平均が x-A, 標本分散が sA2 であり,母集団 B の 標本平均が x-B, 標本分散が sB2 であった.母集団 A,B が共通の 分散 σ2 をもつとすると,その推定値 s2 は 以下のように推定できる.

母集団 A からの標本の偏差平方和:  SA=(nA−1)sA2
母集団 B からの標本の偏差平方和:  SB=(nB−1)sB2
母集団 A,B 全体での偏差平方和:  SSASB =(nA−1)sA2+ (nB−1)sB2
母集団 A,B 共通の標本分散: bunsan

また,母集団Aの標本分布は,N(μA,σ2)であり,母集団Bでは, N(μB,σ2)であることから,それぞれの標本平均は,
x-A 〜 N(μA,σ2/nA), x-B 〜 N(μB,σ2/nB
と分布する.これより,標本平均の差x-Ax-Bは,
diffmean
と分布する.

帰無仮説(H0: μA = μB)のもとでは,μA−μB=0,なので, 標本平均の差は,

diffmean
と分布する.これを標準化した z 値,
diffz
において,標準偏差 σ の代わりに標本標準偏差 s を代入した t 値,
difft
が自由度 nA+nB−2 の t 分布に従うことを利用して検定ができる. なお,母集団 A,B からの標本の大きさがともに等しく, nA=nB=n であるときは,式がずっと簡単になる.

母集団A,Bで分散の同等性が疑われるときは,ウェルチ(Welch)の検定を用いる.

ノンパラメトリック検定

 ノンパラメトリック検定(nonparametric test)とは,母集団分布に関して,正規分布などのある特定の分布 を仮定しないで統計的検定を行う方法である.この手法の利点は,多少の制約がある場合もあるが,どのような 母集団分布からのデータであっても適用可能なことである.
 このため, 標本中に他の観測値から飛び離れた値と思われる 異常値(outlier)が含まれているような場合でも正しい検定を与えることができる.すなわち, 頑健(robust)な検定法である.
一方,弱点としては,分布に関する情報を用いないので,特定の分布の元での最良の検定(正規母集団での小標本 に対する t 検定など)に比べ検定(検出)力(power)が低下することである.  ノンパラメトリック検定においてはデータの値を直接使わず,これを大きさの順に並べてその順位(rank)を 用いることが多い.このことは,データのもつ情報を全部使い切っていない,情報の損失があることを意味する. 他方,異常値の影響は]それだけ受けにくくなっている.  ノンパラメトリック検定の以前のやり方は,n 個の標本が与えられたときに,なんらかの考え方に基づいて 検定統計量(test statistic)T を求める.n が小さいときには,帰無仮説のもとでの T の分布 が直接計算され,その分布に基づいた数表を用いて検定を行った.n が大きいときは,帰無仮説のもとでの T の平均 E[T ],分散 Var[T ] を求め,中心極限定理から,
が標準正規分布 N(0,1) に近似的に従うことを利用して,標準正規分布表を用いて近似的な検定を行った. たとえば,両側検定の場合,|Z | の値が 1.96 以上となれば有意水準 0.05 で帰無仮説は棄却される. なお,T は離散分布となるので近似の精度を上げたいときには分子に 1/2 をいれて,連続性の補正,
を行う.
 現在では,オープンソースの統計環境ソフト R にノンパラメトリック検定を行う関数がかなり実装されている ので,数表をみる必要がないことが多くなった.
 使用法としては,標本の正規性に少しでも疑問がもたれるときは"保険"の意味でノンパラメトリック検定を 用いるとよい.この検定で有意と判断されれば非のうちどころはまったくなくなる.

1標本問題(one sample problem)

概要

 1つの母集団に関する検定問題である.ある連続分布をもつ母集団から大きさ n の無作為標本 X1,…,Xn が抽出されたとする.ここで,母集団分布の位置母数 (location parameter)に関する両側検定問題

H0 : ξ = ξ0,H1 : ξ ≠ ξ0

および,片側検定問題

H0 : ξ = ξ0,H1 : ξ > ξ0

を考えてみよう.正規性の仮定のもとでは,Xi 〜 N( ξ,σ2 ) とおき, ξ と σ2 の推定量として標本平均 x- と標本不偏分散 s2 を用い,t 検定を行う.
 ノンパラメトリック検定では,母集団の位置母数は平均でなくメディアン(中央値)で測られる.これは, コーシー分布のように平均が存在しない分布にも対応するためである.一方,メディアン M
と定義されるので,すべての連続分布で存在する.
 さてこれからは,一般性を失うことなく ξ0 = 0 とする.ξ0 ≠ 0 のときは, Xi のかわりに Xi − ξ0 を考えればよい.

符号検定(sign test)

 いま,X1,…,Xn のうち正の値(+の符号をもつ)の数を T とする.ただし,Xi = 0 となる標本が得られたときは,この標本はなかったものとして 取り除いておく.帰無仮説のもとでは式(5)により,T は2項分布 B(n,1/2) に従う.従って, 両側検定は
という棄却域をつくればよい.c の値は2項分布を計算すれば適当な有意水準に対して正確に 求められる.また,n が大きいときには,
となるので,正規近似(1)を行えば任意の有意水準に対して近似的な検定ができる.片側検定 に対しては,式(7)の棄却域を片側にとればよい.
 メディアンに対する符号検定は 1/2 分位点に対する検定でもある.このことに着目すれば,符号検定 は任意の分位点に対して検定できる.

ウィルコクソン(Wilcoxon)の符号順位和検定(signed rank sum test)

 符号検定は,観測値の大小はまったく無視してその符号だけを用いている.これは観測値の情報 をかなり無駄にしていると考えられる.この点を改良したのが Wilcoxon の符号順位和検定である. しかしながら,母集団分布が対称な分布であるという仮定があらたに加わる.
 標本の絶対値 |X | を大きさの順に順位をつけ,これを Ri とおく. たとえば X1 = 3.1,X2 = -4.5,X3 = 0.9 で あったならば R1 = 2,R2 = 3,R3 = 1 で ある.母集団は連続分布を仮定しているので,|Xi| は確率 1 ですべて互いに異なる値 をとり,Ri は一意的に定まる.
 しかしながら,実際のデータでは測定精度などや測定値の丸めにより,標本値に同じ値が含まれる ことがよくある.このような場合をタイ(tie)があるという.タイがあるときは,タイとなったデータ の順位となり得る値の平均を順位として与えるのが普通である.
 たとえば,X1 = 1.5,X2 = -0.3,X3 = 1.5, X4 = -2.8 であったとする.X1X3 は順位 として 2 か 3 を取ったはずである.そこで,この値の平均 2.5 を順位として与える.この結果, R1 = 2.5,R2 = 1,R3 = 2.5, R4 = 4 となる.このようにすればタイのあるなしにかかわらず, ΣRi = n(n+1)/2 となる.また,Xi = 0 というデータが 得られたときはこのデータはなかったものとして取り除いておくのは符号検定法と同じである.
 さて,Xi > 0 のときの順位和を R+Xi < 0 のときの順位和を R- とおき,この両者のうち小さい 方を T とおく.母集団分布が対称であるので,帰無仮説 H0 : ξ0 = 0 のもとでは R+R- の値はほぼ同じで, ΣRi/2 = n(n+1)/4 に近いことが期待される.よって,T が小さいときには 帰無仮説を棄却することができる.

表1.順位のあり得る組み合わせ(n=4)
 1   2   3   4     R+     R-      T   
1000
911
822
733
644
733
644
555
555
464
373
464
373
282
191
0100

 T の分布ようすを n = 4 のときに考えてみよう.順位の値は 1,2,3,4 の 4 通り ある.ある順位をとったデータが正であったか負であっかはまったくの偶然で決まる.そのあり得る 組み合わせは表1に示すように全部で 16 通りである.どの組み合わせも同じ確率で起こるので, これにより T の分布が決定できる.
 たとえば,T = 0 となるのは 2 通りであるので,その確率は 2/16 = 0.125 と計算される. もし,T = 0 のとき帰無仮説を棄却すると決めたならば,その有意水準は 0.125 となる. このように T の帰無仮説のもとでの分布を計算すれば,検定の有意確率が求まる.
 なお,n が大きいときは,
となることが知られているので,正規近似を用いて近似的な検定ができる.    前節で取り上げた園芸療法での MSE の結果についての解析を再び行う. データを以下の表にまとめた.元データに差と符号と順位も付け加えた.このデータを用いて,園芸療法 を施す前と後では,心理機能に差がないという帰無仮説,すなわち,

H0:差 = 0

の検定を行う.
符号検定
 差が 0 のデータを取り除くとデータ総数は16となった.16 例のうち,正の符号が 7,負の符号が 6 で あるので,2項分布を用いて検定が行える.表2のデータが csv ファイル形式 で R の作業フォルダに 'mse.csv' という名前であるとする.符号検定に 対する R のスクリプトは以下のようになる.
mse <- read.csv("mse.csv") #データ読み込み
d <- mse$after - mse$before #処方後得点−処方前得点
np <- length(d[d > 0])
nm <- length(d[d < 0])
binom.test(c(np,nm))
結果:'+'の確率 = 0.5625,p 値 = 0.8036 で帰無仮説は棄却されず, 園芸療法が心理機能に変化を引き起こすとは認められなかった.
Wilcoxon の順位和検定
 差が 0 のデータを取り除くとデータ総数は16となった.R+ = 86, R- = 50 となった.検定の R のスクリプトは以下のようになる.
wilcox.test(d) #Wilcoxon 検定を行う関数
結果:p 値 = 0.363 で帰無仮説は棄却されない.
t 検定
 比較のために t 検定も行った.
t.test(d)#t 検定を行う関数
結果:平均 = 0.895,t 値 = 1.1648(df = 18),p 値 = 0.2593 で帰無仮説は棄却されなかった.
データの吟味
 今回の調査では,クライアント全体では園芸療法の効果は認められなかった.また,MSE の満点が30点なので, 元々心理機能に問題がなく満点近い得点のクライアントでは,効果が認められないのも当然といえる. その上,園芸療法処方の前後で半年 程度のタイムラグがあるので,何もしなくても認知症の症状が進行する場合もあり,効果が検出しずらい こともある程度予測できる.そこで,データをよく調べるため,差と年齢のプロットを図1に行った.


図1.園芸療法の”効果”と年齢との関係

 これをみると,年齢が高いほど効果が減少しているようにみえる.そこで,多少「後づけ」ではあるが クライアントのうち若年層( 85 歳以下)を取り出すことにした.若年層は 10 例となり,サンプルサイズは かなり小さくなった.
符号検定
 表3のデータで帰無仮説の検定を行う.R のスクリプトは以下の通り.
young <- (1:19)[mse$age<86] #年齢86歳未満のデータ番号を取得
d2 <- d[young]
np <- length(d2[d2 > 0])
nm <- length(d2[d2 < 0])
binom.test(c(np,nm))
結果:'+'の確率 = 0.8889,p 値 = 0.03906 で,帰無仮説は 5%有意で棄却された.これより, 高齢でない層では,園芸療法が心理機能に変化を引き起こすと認められた.
Wilcoxon の順位和検定
R+ = 42.5,R- = 2.5 となった. R のスクリプトは以下の通り.
wilcox.test(d2)
結果:p 値 = 0.01955 で帰無仮説は 1% 有意ではないが,5% 有意で棄却された.
t 検定
t.test(d2)
結果:平均 = 2.8,t 値 = 3.3308(df = 9),p 値 = 0.008788 で帰無仮説は 1% 有意で棄却された.

2標本問題(two sample problem)

概要

 2つの母集団から無作為標本 X1,…,XmY1,…,Yn が抽出されたとする. 等分散の正規性の仮定のもとでは,

Xi 〜 N( μx,σ2 ), Yj 〜 N( μy,σ2 )

と仮定し,帰無仮説は,

H0 : μx = μy
となる.標本平均 X-Y- とこみにした標本不偏分散 s2 を用いて t 検定を行う.
 ノンパラメトリック検定においては,無作為標本 XiYj がそれぞれの 累積分布関数 Fx(・),Fy(・) から抽出されたと考える.帰無仮説は, すべての z に対して

H0Fx(z ) = Fy(z ) )
となる.

Mann-Whitney の U 検定

 順位和検定(rank-sum test)とも呼ばれる.2 つの母集団分布の平均的な位置の違いを検出する方法である. ある程度大きな標本数では t 検定の約 95% の検出効率である.  標本 X1,…,XmY1,…,Yn をまとめたものを Zi,i=1,…,m+n と 表す.すなわち,
ZiXi,i=1,…,m, Zi+mYi,i=1,…,n
である.ここで,Z を小さいものから順に順位をつけこれを Ri とおく.もし, タイがあったときは前節と同様に平均順位を割り付ける.ここで,
とおく.すなわち,RxXi の順位の和であり, RyYi の順位の和である. RxRy = (m+n+1)(m+n)/2 である.
 さて,UxXY より大きくなる個数とする.いま, Xi の順序統計量(大きさの順に並べたもの)を, X1',…,Xm' とし,これに対応する順位を R1',…,Rm' とおく.明らかに X1' は R1'−1 個の Y より大きく, X2' は R2'−2 個の Y より大きい.同様に考えて, Xm' は Rm'−m 個の Y より大きい.よって,
となる.同様に,
であり,UxUy = mn が成立する.検定統計量 UUxUy のうち小さい方とする.X の分布と Y の分布が離れていれば U の値は小さくなるはずである.
 m,n の値が大きいときは,式(12)の帰無仮説のもとで,
または,
となることが知られているので,正規近似(1)を行って検定することができる.なお,この正規近似は m,n が 7 より大きければかなり正確であることも示されている.

Kolmogorov-Smirnov の 2 標本検定

 累積分布関数 F(・),G(・) をもつ 2 つの母集団からそれぞれ大きさ m,n の無作為 標本 X1,…,XmY1,…,Yn が抽出されたとする.このとき,
を経験(標本)累積分布関数という.標本 Yi に対する経験累積分布関数を Gn(x ) としたとき,2 つの母集団からの経験累積分布の最大偏差
を求める.次に,
を求め,これを検定統計量として用いる.
 2 つの母集団の累積分布関数が等しいという帰無仮説のもとでの K の漸近分布 (m,n → ∞ としたときの分布)から以下で定義される臨海値が求められている.すなわち,
であり,この値は表3に示されている.

表3.Kolmogorov-Smirnov 2 標本両側検定の臨界値
   確率:γ       0.99      0.95      0.90       0.85      0.80   
   臨界値:kγ       1.63      1.36      1.22       1.14      1.07   

3-4.解析例

 前節では,心理機能検査(MSE)データを対のある 2 標本データとして解析した.この節では, このデータの園芸療法前データを対照とし,園芸療法後データを処理とみなして解析してみる.この場合, 実際とは異なりこのデータの解析法としては正しくないが, 対照 19 名,処理 19 名のデータとみなすわけである.
データ分布
 対照区と処理区とのデータ分布の違いを箱ひげ図でみる.R のスクリプトは以下のようになる.これを 見ると,対照区と処理区で分布に大きな違いがないようにみえる.検定で確認してみよう.
boxplot(mse$before, mse$after, names= c("対照区","処理区"), ylab="MSE 得点", cex.axis=0.8)
Mann-Whitney の U 検定
 R のスクリプトは以下の通り.
wilcox.test(mse$after, mse$before)
結果:p 値 = 0.6812 で帰無仮説は棄却されない.
t 検定
 R のスクリプトは以下の通り.
t.test(mse$after, mse$before, var.equal=T)
結果:t 値 = 0.3725(df = 36),p 値 = 0.7117 で帰無仮説は棄却されない.
Kolmogorov-Smirnov 検定
 まず,対照区と処理区の経験累積分布の重ね合わせのグラフを書き,両者の乖離の程度をみて検定を行う. R のスクリプトは以下の通り.
plot(ecdf(mse$before), do.points=F, verticals=T,xlab="MSE 得点", ylab="累積確率", xlim=range(mse$after,mse$before), main="")
plot(ecdf(mse$after), do.points=F, verticals=T, add=T, col.h='red', col.v='red')
title(main="MSE 得点の経験累積分布")
legend(locator(1), legend=c("処理","対照"), lty=c(1,1),col=c("red","black"))
ks.test(mse$after, mse$before)
結果:D = 0.1053,p 値 = 1 で帰無仮説は棄却されない.経験累積分布をみても両者はほとんど 重なっている.
「若年層」のデータ分布
 前節の解析では,「若年層」のデータでは,園芸療法の効果の有意性が示された.対照区と処理区とみなした場合 はどうなるか試してみる.R のスクリプトは以下の通り.これをみると,処理区の方が若干上の方にシフトしている ようにみえる.検定で確認する.
boxplot(mse$before[young], mse$after[young], names= c("対照区","処理区"), ylab="MSE 得点", cex.axis=0.8)
title(main="高齢でない層の MSE の得点分布", cex.main=1.0)
Mann-Whitney の U 検定
 R のスクリプトは以下の通り.
wilcox.test(mse$after[young], mse$before[young])
結果:p 値 = 0.4452 で帰無仮説は棄却されない.
t 検定
 R のスクリプトは以下の通り.
t.test(mse$after[young], mse$before[young], var.equal=T)
結果:t 値 = 0.7368(df = 18),p 値 = 0.4707 で帰無仮説は棄却されない.
Kolmogorov-Smirnov 検定
 まず,対照区と処理区の経験累積分布の重ね合わせのグラフを書き,両者の乖離の程度をみて検定を行う. R のスクリプトは以下の通り.
plot(ecdf(mse$before[young]), do.points=F, verticals=T,xlab="MSE 得点", ylab="累積確率", xlim=range(mse$after[young],mse$before[young]), main="")
plot(ecdf(mse$after[young]), do.points=F, verticals=T, add=T, col.h='red', col.v='red')
title(main="高齢でない層の MSE 得点の経験累積分布", cex.main=0.8)
legend(locator(1), legend=c("処理","対照"), lty=c(1,1),col=c("red","black"))
ks.test(mse$after[young], mse$before[young])
結果:D = 0.3,p 値 = 0.7591 で帰無仮説は棄却されない. 経験累積分布をみると,両者に多少の「ずれ」が認められるが,有意な「ずれ」ではなかった.

正規母集団の母分散に対する検定

母集団分散の検定

  平均 μ,分散 σ2 がともに未知である正規母集団に対して,
帰無仮説 H0: σ2 = σ02
対立仮説 H1: σ2 ≠ σ02
の検定を考える.  いま,母集団から大きさ n の無作為標本 x1x2,…,xn を抽出したところ,標本平均が x-,標本分散が s2 であったとする.すると, 帰無仮説ももとで(under H0),標本分散に関係した量が,
vartest
と自由度 n - 1 の χ2 分布に従うので,U を検定統計量にして検定が行える.

 有意水準 5 %の検定は,自由度 n - 1 の χ2 分布の 2.5%点と 97.5%点をそれぞれ χ2(n - 1)0.025, χ2(n - 1)0.975 とすると,

U < χ2(n - 1)0.025U > χ2(n - 1)0.975
のいずれかの不等式を満たしたとき帰無仮説を棄却し,母分散は σ02 と有意に異なる と結論づける.

2つの母集団分散の同等性の検定

 2つの母集団 A,B があり,それぞれが平均を μA,μB, 分散を σA2,σB2 の正規分布に従って いるが,その値は未知であるとする.このとき,2つの母分散の同等性の検定,

帰無仮説,H0: σA2 = σB2
対立仮説,H1: σA2 ≠ σB2

の検定を考える.

 母集団 A から大きさ nA,母集団 B から大きさ nB の標本を抽出した. 母集団 A からの標本の標本平均が x-A, 標本分散が sA2 であり,母集団 B の 標本平均が x-B, 標本分散が sB2 であるとする.すると, 標本分散に関係した量はそれぞれ

vartest
と χ2 分布に従い,それぞれが独立である.これらの量の比は,
vartest
のように,自由度 nA - 1,nB - 1 の F 分布に従う.

 ところで,帰無仮説が正しいとする と,σA2 = σB2 とおけるので, 母集団の分散比は,γ0 = σA2B2 = 1, となる.このとき,標本分散の分散比の統計量 γ が,

vartest
と,自由度 nA - 1,nB - 1 の F 分布に従うので,この γ 値を検定統計量にして 2つの母分散が等しいという帰無仮説の検定が行える.
 すなわち,有意水準 5 %の検定を行うには,自由度 nA - 1,nB - 1 の F 分布 の 2.5%点と 97.5%点をそれぞれ F(nA - 1,nB - 1)0.025, F(nA - 1,nB - 1)0.975 とすると,検定統計量 γ が,
γ < F(nA - 1,nB - 1)0.025γ > F(nA - 1,nB - 1)0.975
のいずれかの不等式を満たしたとき帰無仮説を棄却し,2つの母集団の分散は有意に異なると結論づける.

2つの母集団の分散比の信頼区間

 2つの母集団 A,B の分散 σA2,σB2 の分散比, γ0 = σA2B2,の 95%信頼区間は, 上記の考えから簡単に求めることができる.すなわち,互いに独立に χ2 分布する 変量の比が,標本分散の分散比 γ と母集団分散比 γ0 の比となり,
vartest
と分布する.これより,母集団分散比の 95%信頼区間は,
vartest
となる.

課題

  • 問題4−2: 淡水性ウナギの汽水域での生理活性の違いのデータ、ウナギデータ, を2標本データとして違いを検定せよ.ただし,2母集団で分散の違いがあるようなので,まず, 分散の同等性の検定を行い,違いが認められるときは,データの対数変換を行うなどして,分散安定化 変換を行って t 検定を行ってみよ.また,データの対数変換などを行わず.ウェルチ検定やノンパラメトリック 検定を行った場合との比較検討を行え.

    成功確率(比率)に関する検定

    標準正規分布による近似検定(大標本理論)

     成功確率 p のベルヌイ試行を n 回行ったときの成功回数 X は, X 〜 B(n, p),のように2項分布に従う.X の平均と分散はそれぞれ, E[X ] = np,Var[X ] = np(1 - p),である.
     ここで,成功確率が p0 であるという帰無仮説,
    H0: p = p0
    の検定を考える.帰無仮説のもとでは,成功回数 X は,X 〜 B(n, p0), と分布するので,X をその平均と標準偏差で標準化すると,中心極限定理から,
    binotest
    のように標準正規分布に漸近的に従う.
     これより,近似的な 5%両側検定は,標準正規分布の 97.5%分位点の z0 = 1.96 より 検定統計量 T = |z| の値が大きくなったとき帰無仮説を棄却することで得られる. なお,二項分布は離散的なので,Yates(イエーツ)の連続性の補正を行った検定統計量を用い,
    binotest
    のとき帰無仮説を棄却する方が近似の精度がよいと言われている.

     このように,中心極限定理を利用して,標準正規近似を行って検定を行うやり方を大標本(large sample)理論 といい,コンピュータが発達する以前はもっぱら大標本理論に基づいた検定を行っていたが,現在では正確な確率 が短時間で計算できるので,大標本理論の重要性は大きく低下したといえる.

    比率の正規近似に基づく信頼区間

     成功確率 p のベルヌイ試行を n 回行ったとき x 回成功したとすると,成功確率は, p^ = x/n,と推定される.この推定値は最尤推定値である.
     成功回数 x は二項分布し,その平均は E[x ] = np,分散は Var[x ] = np(1 - p),で あるので,成功確率推定量 p^ の平均は E[p^ ] = E[x/n] = p, 分散は Var[p^ ] = Var[x/n] = Var[x ]/n2 = p(1 - p)/n, となる.これより,
    binotest
    と漸近的に分布するので,標準正規分布の 97.5%点の z0 = 1.96 を用いると, 近似的に
    binotest
    という不等式が成り立つ.これを整理すると,
    binotest
    という p の2次不等式を解くことに帰着する.いま,p の2次方程式の根を
    binotest
    とすると,この根を用い,p の 95%信頼区間は近似的に
    binotest
    となる.
     また,連続性の補正を行うには,成功確率の推定値 p を,信頼区間の下限と上限でそれぞれ
    binotest
    というように変えて,信頼区間が少し広くなるようにする. R では,これらの式を用いて信頼区間を構成しているようである.

     ここで,さらに近似を加えて,z02 の項を消去すると, p の2次方程式の根は,

    binotest
    となるので,p の 近似的な 95%信頼区間は,
    binotest
    と簡略化される.
     なお,この信頼区間は, 成功確率推定量 p^ の分散において,真の成功確率 p の 代わりにその推定量 p^ に置き換えて,Var[p^ ] = p^(1 - p^)/n,とみなした場合と同じで, この信頼区間は教科書等でよく出てくる.
     簡略化された信頼区間で連続性の補正を入れるには,
    binotest
    として,信頼区間の幅を拡げる.

     ところで,正規近似による信頼区間の構成では,場合により信頼区間が負になったり 1 を超えることがあるが, このときは,0 と 1 で切り詰める.

    二項確率の計算による正確な検定

     現在では二項確率が R などのコンピュータソフトにより直接計算できるので, 正規分布による近似検定を行う意味はあまりないといえる. 直接計算のやり方を簡単な例で考えてみる.
     成功確率の正確な検定の例として,n = 10 回の試行で x = 7 回の成功が観察されたベルヌイ試行で,
    帰無仮説:H0: p = 0.4
    対立仮説:H1: p ≠ 0.4
    の検定を行ってみる.
     まず上側確率を求める.帰無仮説のもとで,x = 7 回の成功が得られる確率を q7 とすると, これ以上の成功回数 8 回,9 回,10 回が得られる確率,q8,q9,q10, を加え合わせる.これらは,
    binotest
    となり,上側確率は,
    binotest
    となる.
     下側確率は,帰無仮説のもとで,x = 0,1,…,回の成功が得られる確率で q7 より 小さいものを, q0,q1,…,と求める.これらは,
    binotest
    となり,q2 > q7 なので,求める下側確率は,
    binotest
    となる.
     これより p 値は,
    Pupper + Plower = 0.05476188 + 0.0463574 = 0.1011193
    となる.

    適合度検定

    Pearson χ2 検定

     前節の比率の検定は,χ2 分布を用いる適合度検定と大きな関係がある. ここでは n 回のベルヌイ試行で X 回成功したときに,成功確率が p0 であるという,
    帰無仮説,H0: p = p0, 対立仮説,H1: p ≠ p0
    の検定を考えた.そこでは,X を標準化して標準正規分布にもって行ったが,これを2乗して χ2 分布を用いることもできる.すなわち,
    fit
    という関係がある.
     ところで,n 回のベルヌイ試行の結果と帰無仮説のもとでの期待値を表にすると,

       成 功   失 敗 
     観測度数  X n - X
     期待度数  np0 n(1 - p0)

    となる.ここで,ピアソン(Pearson)のχ2 値,

    fit
    を計算すると,
    fit
    となる.つまり,χ2 値は,試行回数 n が大きくなるにつれて 帰無仮説のもとで自由度 1 の χ2 分布に漸近的に従う.よって,これより検定が行える.

    確率分布との適合度

    確率分布が既知のとき

     データが想定している確率分布に適合しているかは,ピアソン(Peason)の χ2 適合度検定で行う ことができる.いま,離散分布の,たとえば m = 5 のセルに対して,観測されたカウントデータと対応する 想定確率が,

      セル1  セル2 セル3  セル4  セル5    計   
     観測度数  n1 n2 n3 n4 n5 n
     想定確率分布  p1 p2 p3 p4 p5 1

    のようになっていたとする.このとき,ピアソン(Peason)の χ2 値は,
    tekigo
    のように近似的に自由度 m - 1 の χ2 分布に従う.これにより,データが想定確率分布に 適合しているかの検定が行える.検定の帰無仮説は,

    H0:データは想定確率分布に従う.

    である.  この近似は n が大きく,各セルの度数 ni が 5 以上である ことが望ましい,とされている.  一方,正規分布などの連続分布では,適当に階級分けして離散化すればこの検定が行える.ただし,階級分け は任意なので,階級分けのやり方によっては結果が異なる恐れがある.

    確率分布のパラメータをデータから推定する場合

     確率分布のタイプ(二項分布やポアソン分布など)は想定できるが,パラメーターはデータから推定することが 普通であろう.このときは,推定されたパラメーターのもとでの推定確率分布を用いて,セル数が m = 5 のときは,

      セル1  セル2 セル3  セル4  セル5    計   
     観測度数  n1 n2 n3 n4 n5 n
     推定確率分布  p^1 p^2 p^3 p^4 p^5 1

    のような表ができる.推定したパラメーターの数が k であったとすると, このとき,ピアソン(Peason)の χ2 値は,
    tekigo
    のように近似的に自由度 m - k - 1 の χ2 分布に従う.これにより,データが推定確率分布に 適合しているかの検定が行える.

    参考文献

    1. 心理・教育のための統計法(第 2 版),山内光哉,1998,サイエンス社
    2. 工学のためのデータサイエンス入門−フリーな統計環境Rを用いたデータ解析−,間瀬茂ら,2004, 数理工学社
    3. 実践生物統計学−分子から生態まで−(第 1 章,第 2 章), 東京大学生物測定学研究室編(大森宏ら), 2004,朝倉書店
    4. R で学ぶデータマインニング I −データ解析の視点から−,熊谷悦生・船尾暢男,2007,九天社
    5. R で学ぶデータマインニング II −シミュレーションの視点から−,熊谷悦生・船尾暢男,2007,九天社

    Copyright (C) 2008, Hiroshi Omori. 最終更新:2008年 6月 4日