狠狠撸

RとStanで
一元配置分散分析
ベイズ推定を用いて3群の平均値差を比較
千葉県東上総児童相談所
児童心理司佐名隆徳

自己紹介
? 関西大学卒業（単位ギリギリ）→県の心理職（公僕）
? 大学院には正直行きたかったけどお金がないので断念した
? 延長因子分析とSEMの小包化について興味があるが資料がない
? 好きな女性有名人：多部未華子と黒田エイミ

分散分析とは
? 質的な独立変数（ｘ）の値によって、従属変数（ｙ）の平均がどのように異
なるかを分析するための方法
? 質的な独立変数の値は、その要因の水準とよばれる（ビールの種類、温
度、など）
? 対応のない要因：要因の各水準に異なる被験者が割り当てられる
?異なる水準に含まれる従属変数の値は互いに独立
? 対応のある要因：異なる水準に含まれる従属変数の値に相関を想定

分散分析なのに平均値差なの？
? 分散を用いて、多群の平均値を分析するから分散分析
? 分散分析は、F値（＝要因変動/誤差変動）を求めることで、その要因の
効果があるかどうか（有意かどうか）を判断する
? F値は大きい（要因変動の割合の方が大きい）ほど、p値が小さくなる
?「要因の効果が○％水準で有意」といえる値になっていく

分散分析表（２要因）例

なぜRとStan？
? MCMCソフトウェアであるStan(無料)を、統計ソフトであるR(無料)上で用いて分析する
? MCMCは、ベイズ推定において事後分布等のサンプリングを行うアルゴリズムの総称
? ベイズ推定の主目的：①母数の事前情報を推定に用いる、②サンプルサイズが少な
い場合にもそれに直接影響されない母数推定を行う
? ベイズ推定では、母数についての情報を事前分布として設定し、それを尤度関数と
合わせて各母数に関する事後分布に基づいて推定値を得る

余談
? 実は分散分析って、大学2年かそこらの実験実習（まじめに受けてなかっ
た）でやったくらいで、正直よく分かってない
? しかも実験系の統計は苦手
? でも最近よく分散分析をベイズでという質問を受けるので、仕事のストレ
スでアパシー状態になりながらも今回作成してみた

注意点？
? そもそも分散分析をベイズでやるのってどうなの？
ベイズの利点として、自由なモデリングが可能なため分析における制約
が少ないことが挙げられる。
だけど分散分析は強い仮定の下で分析を行うため、ぼんやりとベイズ
のメリットの１つを殺している気がしなくもない。
まあでも正直勉強不足過ぎてよくわかんない。

分析方法
? 対応あり?一元配置分散分析をベイズ推定で実施
? 対応のある場合は設定する事前分布に多変量正規分布を用いるため、
Stanコードがちょっと分かりづらい（文系の限界）
? 堅く言うと、各水準の差の事後分布を確認することが目的になる
? Stanコードで各水準の差について計算する指示を行うため、多重比較的
な方法を最初から分析に組み込める

使用データ（仮想データ）
? N=10(スーパードライ、ギネス、よな
よな)
? 変数は「各ビールを飲んでみてうま
いと思った得点」（高いほどうまい）
? 目的：各ビールのうまい得点の平
均値に差があるかどうかを分析する
（ちなみに僕は、よなよなエールうま
いと思っています）
id sup gi yonayona
1 3 6 1
2 3 9 1
3 3 8 2
4 4 5 1
5 8 9 1
6 7 3 2
7 8 5 3
8 7 6 4
9 7 8 3
10 9 9 5

使用パッケージ
library(ggplot2)
library(StanHeaders)
library(rstan)
library(lordif)
library(psych)
（下２つは不要かも???）

データの取り込みや命名、他
dat3<-read.csv("beer3.csv")
x1<-dat3$sup #x1はスーパードライ
x2<-dat3$gi #x2はギネス
x3<-dat3$yonayona #x3はよなよなエール
x <- structure(
.Data = c(x1, x2, x3), .Dim = c(10,3))
#ベクトル型に渡すために、structuer()関数を用いてデータを整形＋.Dim因数でデータの構造
を指定(10×3)＋list()でデータリスト型に変換

（xの中身はこんな感じ）
> x
[,1] [,2] [,3]
[1,] 3 6 1
[2,] 3 9 1
[3,] 3 8 2
（略）
[9,] 7 8 3
[10,] 9 9 5

データリスト作成など
n <- nrow(x) #nはxの行数、という命令
> n
[1] 10
data <-list(N = n, X = x)
#Stanで分析に使用するため、list()関数を用いてNとXのデータリストに格納

（dataの中身はこんな感じ）
> data
$N
[1] 10
$X
[,1] [,2] [,3]
[1,] 3 6 1
[2,] 3 9 1
（略）
[9,] 7 8 3
[10,] 9 9 5

Stan関連コード dataブロック
stancode <- "
data{
int<lower=0> N;#「N」っていう整数があるぞ。Nは最低値が０だぞ。
vector[3] X[N];#要素が３つのベクトル「X」があるぞ。Xは０～Nだぞ。
}

parametersブロック
parameters{
vector[3] mu;#muは３つのベクトルがあるぞ。
real<lower=0> sig[3];#sigは３つのベクトルがある、０以上の実数だぞ
real<lower=-1,upper=1> rho;#群内の要因が加わる混合計画では相関を
加える
}

transformed parametersブロック
transformed parameters{
matrix[3,3] Sigma;#Sigmaは[3,3]行列です。あとで定義。
Sigma[1,1] = sig[1]*sig[1];
Sigma[1,2] = sig[1]*sig[2]*rho;#共分散行列の非対角成分に相関を掛ける
Sigma[1,3] = sig[1]*sig[3]*rho;
Sigma[2,3] = sig[2]*sig[3]*rho;
Sigma[2,1] = Sigma[1,2];
}

modelブロック
model{
X ~ multi_normal(mu, Sigma);#Xは多変量正規分布に従う
rho ~ uniform(-1,1);
sig ~ cauchy(0,2.5);# sigは平均0,SD2.5のコーシー分布に従う。
}
#多変量正規分布のパラメータは、K：正の整数、μ：長さKのベクトル、σ：K*Kの分散共分
散行列になる
# 事前分布の影響がありそうな小さなサンプルに対しては，分散の事前分布として半コー
シー分布を選ぶと良いらしい（※半コーシー分布については他参照）

generated quantitiesブロック
generated quantities{
real delta21;#deltaという新たなパラメータを定義
real delta13;
real delta23;
delta21 = mu[2] - mu[1];#deltaで平均値差を算出する、と定義
delta13 = mu[1] - mu[3];
delta23 = mu[2] - mu[3];
}

MCMCを実行
fit <- stan(model_code = stancode, iter = 10000, chains = 1, data = data)
#<-stan()：()内で指定したMCMCを実行するというコマンド
#iter：サンプリング数。各chainsでサンプリングを何度行うか
#chains：マルコフ連鎖数。ホントは４～５が望ましい
#data：この場合はスライド１４で作成したデータリスト”data”を利用する

print(fit, digits_summary = 2)
>
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[1] 5.93 0.01 0.85 4.23 5.38 5.92 6.47 7.66 4037 1
mu[2] 6.84 0.01 0.81 5.25 6.33 6.83 7.35 8.41 3905 1
mu[3] 2.32 0.01 0.51 1.32 1.99 2.31 2.63 3.37 3664 1
sig[1] 2.64 0.01 0.68 1.68 2.16 2.52 2.96 4.36 3180 1
sig[2] 2.47 0.01 0.65 1.54 2.01 2.36 2.80 4.00 3623 1
sig[3] 1.61 0.01 0.44 1.01 1.30 1.52 1.83 2.67 3146 1
rho 0.26 0.00 0.20 -0.12 0.12 0.25 0.40 0.64 3897 1
（共分散は省略）
delta21 0.91 0.01 1.00 -1.07 0.26 0.92 1.57 2.81 5000 1
delta13 3.61 0.01 0.89 1.86 3.03 3.61 4.17 5.41 5000 1
delta23 4.52 0.01 0.84 2.87 3.98 4.53 5.06 6.16 5000 1
lp__ -36.19 0.05 2.00 -40.93 -37.29 -35.84 -34.72 -33.32 1817 1

解釈の仕方
#print()で表示された事後分布を確認。
#delta21は、平均で0.91の差があり、95％確信区間は[-1.07~2.81]
#delta13は、平均で3.61の差があり、95％確信区間は[1.86~5.41]
#delta23は、平均で4.52の差があり、95％確信区間は[2.87~6.16]

今後
? 交互作用にも言及した２要因の分散分析
（温度と種類など）をベイズに適用するケー
スについて書くつもり。（→みたいなやつ）
? ケースワーク技法のエビデンス検証にベイ
ズを適用する、なんてプロジェクトを遂行中
なので、いい結果が出るといいなと思う。
? 自由にモデリングしてゴリゴリ分析、みたい
な日はやってくるのだろうか…。
id condition sup gi
1 hot 3 6
2 hot 3 9
3 hot 3 8
4 hot 4 5
5 hot 8 9
6 cool 7 3
7 cool 8 5
8 cool 7 6
9 cool 7 8
10 cool 9 9

狠狠撸

搁と厂迟补苍で分散分析

More Related Content

搁と厂迟补苍で分散分析