ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2017年12月09日

回帰分析の悩みどころ (「アヒル本」7.1-7.5) [スライド紹介]

 こんにちは。mutopsyです。この記事は,『StanとRでベイズ統計モデリング』読書会 (Osaka.Stan #4) で使用したスライドの紹介記事です。

 『StanとRでベイズ統計モデリング』,通称「アヒル本」のChapter 7「回帰分析の悩みどころ」の前半部分を解説しているスライドです。回帰分析を用いる際に留意するべき点,あるいは工夫できる点として,「交互作用」「対数をとるか否か」「非線形の関係」「多重共線性」「交絡」の5つに焦点を当てています。内容に関してはスライドをご覧頂くか,「アヒル本」を読んで頂ければ良いかと思いますので,この記事ではスライド内で紹介しているStanコードを(コピペしやすいように)改めて紹介したいと思います。


 以下のStanコードは,2つの説明変数(イケメンか否か・年収)で1つの応答変数(魅力)を予測する重回帰モデルです(この例は清水先生の記事を参考にしました)。まずは交互作用項を含まないモデル(すなわち,交互作用に関わるパラメータが0であると仮定するモデル)から。

Stanコード01: 交互作用を含まないモデル (スライド11枚目)
data {
  int N;  //データの数
  int<lower=1, upper=10> Miryoku[N];  //魅力の評定値データ
  int<lower=0, upper=1> Ikemen[N];  //イケメンか否かの2値データ
  real<lower=0> Nenshu[N];  //年収データ
}

parameters{
  real b[3];  //回帰係数 b[1], b[2], b[3]
  real sigma;  //残差の標準偏差
}

model{
  for (n in 1:N){
    Miryoku[n] ~ normal(b[1] + b[2]*Ikemen[n] + b[3]*Nenshu[n],  //重回帰モデル
		sigma); 
  }
}

 このコードのモデル部分を以下のように書き換えて,2つの説明変数を掛けたものを新たな項として加えることで,交互作用項を含むモデルになります。

Stanコード02: 交互作用を含むモデル (スライド13枚目)
model{
  for (n in 1:N){
    Miryoku[n] ~ normal(b[1] + b[2]*Ikemen[n] + b[3]*Nenshu[n]
		+ b[4]*Ikemen[n]*Nenshu[n],  //この部分が交互作用項を表す
		sigma);
  }

 また,以下のtransformed paramatersブロックをmodelブロックの前に挿入することで,回帰係数の和についても推定することができます。これは,年収が魅力に与える効果の大きさ(傾き)を,イケメンと非イケメンのそれぞれについて推定することに相当します。

Stanコード03: パラメータの和を推定 (スライド16枚目)
transformed parameters{
  real Intrcpt_I0;		
  real Slope_I0;
  real Intrcpt_I1;
  real Slope_I1;

  Intrcpt_I0 = b[1];
  Slope_I0 = b[3];
  Intrcpt_I1 = b[1] + b[2];
  Slope_I1 = b[3] + b[4];
}


 非線形な関係を表すモデルも,モデル式を変えてやれば簡単に表現できます。

Stanコード04: 二次曲線の当てはめ (スライド33枚目)
data {
  int N;  //データの数
  real X[N];  //変数X(説明変数)
  real Y[N];  //変数Y(応答変数)
}

parameters{
  real a;
  real b;
  real x0;
  real<lower=0> s_Y;
}

model{
  for (n in 1:N){
    Y[n] ~ normal(a + b*(X[n]-x0)^2, s_Y);
  }
}

Stanコード05: 時系列データへの指数曲線の当てはめ (スライド35枚目)
data {
  int N;
  real X[N];
  real Y[N];
}

parameters {
  real<lower=0, upper=100> a;
  real<lower=0, upper=5> b;
  real<lower=0> s_Y;
}

model {
  for (n in 1:N)
    Y[n] ~ normal(a*(1 - exp(-b*X[n])), s_Y);
}


 最後に,多重共線性がある場合(=説明変数同士が強く相関している場合)に何が起こるかをシミュレーションするためのRコードを紹介します。ここではベイズではなく普通にlm()を使って得られる結果を用いています。パラメータをいろいろ変えてみながら遊んでみると楽しいです。

Rコード01: 多重共線性のあるデータに重回帰分析を適用するとどうなるか (スライド42枚目)
##パラメータ等の指定##
rAB <- 0.9   #説明変数間の母相関係数。この値をいろいろ変えてみる。
rAY <- 0.5  #AとYの母相関係数
rBY <- 0.6  #BとYの母相関係数
n <- 100  #サンプルサイズ
Rep <- 1000  #サンプリング回数

##分散・共分散行列の作成##
Mat <- matrix(c(1, rAB, rAY, rAB, 1, rBY, rAY, rBY, 1), ncol=3) 

##N = nのデータセットをRep回生成し,それぞれに対して重回帰分析##
Res_beta <- data.frame(b2 = 1:n, b3 = 1:n)
for (i in 1:Rep){
  d <- as.data.frame(mvrnorm(n= n, mu= c(0, 0, 0), Sigma= Mat, empirical= FALSE))
  colnames(d) <- c("A","B","Y")  
  reg <- lm(Y ~ A + B, data = d)
 Res_beta[i,] <- reg$coefficients[2:3]    
}
Res_beta  #Rep回分の偏回帰係数 (b2とb3) が格納されたデータフレーム

 思い返してみると,このOsaka.Stanでの発表が,僕にとってのベイズ(とStan)の始まりだったように思います。Osaka.Stanに参加するまでまともにStanを動かしたことのなかった僕ですが,なんやかんやでそれなりに自由にモデルが書けるようになりました。モデリングの感覚をつかむには,自分で手を動かしてみるのが一番だと思います。その点,アヒル本ではStanコードの各行の意味までしっかりと説明がなされているので,実際にコードを走らせながら読めば相当身になると思います。

posted by mutopsy at 00:00 | 統計

2017年12月08日

ベイズ推定の初歩──二項分布を例に── [スライド紹介]

 こんにちは。mutopsyです。この記事では,約2年前に公開した「ベイズ推定の初歩──二項分布を例に──」というスライドを紹介します。



 この資料は,2015年7月18日に関西学院大学で開催された読書会で用いたものです。文献は“Bayesian Cognitive Modeling: A Practical Course”,通称「コワイ本」と呼ばれる本です。今年の夏,日本語訳が出版されて再び話題に上がりましたね(『ベイズ統計で実践モデリング: 認知モデルのトレーニング』)。
 このスライドでは,コワイ本の第3章「Inferences with binomials(二項分布を使った推論)」について解説をしています。内容としては,ベイズ推定の基本的な考え方と,確率的プログラミング言語の1つであるWinBUGSを使って実際にベイズ推定を行う方法について,二項分布に従う変数を例に説明しています。二項分布は整数値を返す離散分布で,コインを投げて表が出る回数であったり,n問中k問正解する回数であったりと,直感的に理解しやすい分布なので,導入としてよく用いられます。しかし,ここで紹介している基本的な内容は,正規分布のような連続分布に対しても適用可能です。そういった理由もあって,スライドのタイトルを「ベイズ推定の初歩──二項分布を例に──」としました。
 当時の僕は博士前期課程2年で,まだベイズの「ベ」の字も知らないような初心者でしたが,この本を読んでこの資料を作って,なんとなくベイズの「ヘ」ぐらいまでは分かったかなあといったところです。今であれば,『ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版』,通称「犬4匹本」や,『StanとRでベイズ統計モデリング (Wonderful R) 』といった,日本語で読める良著も増えており,ベイズを学ぶ敷居はだいぶ低くなりつつあるように思います。インターネット上にも有志の方々による分かりやすい資料がたくさん公開されています。とてもいい時代ですね。
 ちなみに,コワイ本で用いられているWinBUGSは現在開発が終了していますが,その後継であるJAGSが利用できます。また,より効率的なアルゴリズムを採用しているStanもよく使われています(筆者もStanを愛用しています)。さらに,特定の目的に特化した,簡単に使えるソフトウェアやパッケージも増えつつあります。今やベイズはコワくない。Let's Do Bayesian Data Analysis!
posted by mutopsy at 00:00 | 統計

2017年12月07日

統計ブログを作りました。

こんにちは。mutopsyこと武藤拓之です。
大阪大学で認知心理学の研究をしています。
心理学と統計が好きな,どこにでもいる普通のベイジアンです。

勢いあまって統計ブログというものを作ってしまいました。
RやStanを使って統計に関する記事を書いたり,あるいは書かなかったりする予定です。

このブログの名前はベイズマックス。あなたの知的好奇心をくすぐります。さういふブログに私はしたい。
何卒,よろしくお願いいたします。
posted by mutopsy at 16:08 | 雑記