2018の9日目の記事です。カイ2乗検定のようなことをStanでやってみようという内容です。何か致命的な間違いがあれば教えていただければと思います。パッケージの読み込みなどまずはパッケージの読み込みから。Mojaveのlocale問題の解決)今回,この記事を書くための解析を行なった環境を紹介します。ggplot2_3.1.0[後略]rstan2.18にしていないことがばれますね。それとlocaleが見慣れないものになっています。調べてみると,Rstudioのコミュニティにこんな情報が。日本語だけの問題でもなかったようです。みんな困ってたんですね。いうわけで,localeをUTF-8と日本語を読み込めるように呪文を走らせて変更しておきます。ggplot2_3.1.0[後略]ちゃんとlocaleが変わっていることが確認できます。これでいつも通りにデータの読み込みなどができるはず。それでは,よいRライフを。2018の9日目の記事を終わります。calenderなので,ちゃんと?Stanを使います。2018はこちらカイ2乗検定的なことがしたいさて本題です。カイ2乗検定的なことをしていきます。実は,カイ2乗検定をベイジアンモデリングでやってみるというものはすでにTJOさんが試されています。p値を計算したくなる検定の数々を試しにStanによるベイジアンモデリングで代替してみた今回は,TJOさんの記事をもとに,ほとんど同ことをやってみます。私のようなプログラミングに弱い人間は,Stanコードだけ見ていると,やっていることがわからなくなってしまうことが多いので,Stanコードの中身を極力シンプルにして,自分用の解説をしていこうかと思います。(1990)のデータでやってみます。デンマークの労働者の満足度(低,高)とスーパーバイザーの満足度(低,高)の人数の関係を調べたもののようです。データは以下のような感じです。dat1行目:低い満足度の労働者数。スーパーバイザーの満足度が低い場合162名,高い場合110名2行目:高い満足度の労働者数。スーパーバイザーの満足度が低い場合196名,高い場合247名3行目:スーパーバイザーの満足度の高低。4行目:このまま数値化すると,highが1,lowが2となる。lowを1とし,highを2として,参照点をlowにするためのコード。Stanに渡すデータはこんな感じにします。as.numeric(dat$sv_s))1行目:Nは分割表の行数。2)に変換カイ2乗検定的なことを行うStanコードTJOさんによると,こんな感じの表を分析する場合,二項ロジット分析の枠組みを使うとのことで,今回の場合はこんなモデルになるようです。低い満足度の労働者)のうち,確率pで生じる。確率pは,よくある回帰分析の予測値を確率に返すようにしている処理(たしか逆ロジット変換でよいはず)を行なっています。この辺りがロジット分析の枠組みな感じでしょうか。これをStanコードに書き下します。下記のStanコードを”chi2_binomial01.stan”として作業ディレクトリに保存します。25行目:推定結果を用いて新たにppをN個とdiff_ppという変数を算出する25行目:pp[1]としてスーパーバイザーの満足度が低い群での労働者の満足度が高い群となる事後確率を算出26行目:pp[2]としてスーパーバイザーの満足度が高い群での労働者の満足度が高い群となる事後確率を算出27行目:diff_ppとしてスーパーバイザーの満足度の高い群の事後確率(pp[2])と低い群の事後確率(pp[1])の差分を算出。スーパーバイザーの満足度によって満足度の高い労働者が増えるか検討するさあ推定だあとは推定するだけ。Rhat=1).パラメタaとbのトレースプロット無事に収束しているようですし,bの区間推定値も0を挟んでいないので,労働者の満足度高低で人数の違いはあるいってよさそうです。diff_ppも区間推定値に0を挟まないので,スーパーバイザーの満足度によっても満足度の高い労働者数も増えると言えそうです。一応bの推定値を可視化しておきましょう。Stan!反省文この記事は二番煎じに過ぎず,さらに本来自分でちゃんと確認すべきことが確認しきれておりません。記事にしたものの,なんかもやもやしていてところがあって,まだ理解できた!という感じがありません。factorを出してくれるので,それをStanで再現し,かつ,残差分析までやってみたり,ということにチャレンジしたかったので,分割表をJASPやBayesFactorパッケージでどのように求めているか書かれた論文を読んでいたのですが,サンプルコードがあったrjags試しにやったら推定途中でRstudioがアボーンするし,いろいろと迷子になってしまって今回は諦めました。明らかに鍛錬不足,,,がんばります。コメントをキャンセルメールアドレスが公開されることはありません。コメントデータの処理方法の詳細はこちらをご覧ください。「最強のM-1漫才師は誰だ」シリーズへの挑戦カテゴリ。