ベイズ統計モデリング～単回帰～
今回は単回帰をベイズ統計モデリングによって検証していきます。今回も想像上の研究条件を考え疑似データを作成しました。今回の疑似データでは、「身体活動量とテスト成績の関係」という点に着目しました。自分にとって運動に絡めたほうが理解が進むのでこのような形にしています。イメージとしては↓の図のような感じです。
データセットは正規分布による乱数発生によって作成しました。データセット作成に関するコードは↓となっています。プログラミング言語はRを使っています。
#データセット
ta=rnorm(20,50,10)
tb=rnorm(20,65,10)
tc=rnorm(20,80,10)
ea=rnorm(20,20,5)
eb=rnorm(20,40,5)
ec=rnorm(20,60,5)
data1=data.frame(ta,ea)
data2=data.frame(tb,eb)
data3=data.frame(tc,ec)
names(data1) <- c("score", "physical_activity")
names(data2) <- c("score", "physical_activity")
names(data3) <- c("score", "physical_activity")
data12=rbind(data1,data2)
data_all=rbind(data12,data3)このコードでどのようなデータセットが作成できたか確かめるために散布図を作成しました(回帰直線も加えています)。
想定していた身体活動量(physical_activity)が多いほどテストの成績(score)が高くなっているという疑似データが作成できました。ちなみに図の作成に関するコードは↓になっています。Rのライブラリであるggplot2を用いて図を作成しています。
library(ggplot2)
g = ggplot(data_all, aes(x =physical_activity, y = score))
g = g + geom_point()
g = g + geom_smooth(method = "lm")
g = g + scale_color_nejm()
plot(g)このデータセットを用いてベイズ統計モデリングにおける単回帰を実行してみました。詳しいコードは、【馬場真哉 (2019)．「実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門」講談社】という書籍を参考にしています。そのため、ここではコードの記載は省略します。
単回帰ということで、モデルの構造としては、
ui=β0+βi*xiとして、yiは平均がuiで分散がσ^2に従うと想定しています。yiがテストのスコア、xiが身体活動量となっています。式にすると、
scorei ～ Normal(Intercept+beta*physical activity, sigma^2)となりますね。
実際にMCMCによる単回帰のベイズ統計モデリングを行ってみると↓のような結果になりました。
Rhatは1となっており、上手く収束できていますね。モデル構造の切片(Intercept)は36.65となっており、Betaは0.73となっています。このようなことから、score(テストの成績)=36.65+0.73*physical activity(身体活動量)という関係があることが判明しました。
最後にパラメータ(Intercept、beta、sigma)の事後分布とトレースプロットを確認しておきます。
今回は正規分布の乱数発生により作成した疑似データということもあり、事後分布もきれいにモデリングされていますね。収束も問題ないようです。
今回はここまでとします。最後まで読んでくださりありがとうございました。
ベイズ統計に関してはまだまだ勉強中の身ですので、データの記載や解釈等に間違いがありましたら、ぜひコメント等で教えてください。
参考文献
・馬場真哉 (2019)．「実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門」講談社