橘田いずみさんが紹介する「絶対外さない！お取り寄せ餃子5選」

こんにち餃子，ということで2016年1月22日放送「白熱ライブ　ビビット」より。
絶対外さない！お取り寄せ餃子5選として，餃子評論家の橘田いずみさんが餃子を紹介されていました。
今回のポイントは，“おうちで焼いてもお店の味をそのまま再現できる”ことだそうです。
わたくし，個人的に橘田いずみさんを応援しています。

橘田いずみのザ・餃子

1．餃子激戦区　宇都宮のお取り寄せ餃子

曰く，見た目がびっくりで若い人におすすめ，彩花というお店ののキング餃子（5個入り450円）。
普通の餃子の約3倍らしく，一時期健啖家を自負していた自分としては非常に胃袋を刺激されました。
ところで橘田さん，宇都宮のお取り寄せ餃子のすべてを食べ比べたとのことで，さすがの貫禄です。

2．餃子激戦区「浜松」　10個＿＿＿円のお取り寄せ餃子

浜松餃子栄福というお店の餃子（10個300円）。
材料を手で粗めに切ることでしっかりとした食感が残る餡となり，隠し味の山芋も合間ってとてもおいしいそうです。
この値段と数はいっぱい食べたいおひとりさまからわいわい食べたい複数人まで，幅広く楽しめて素敵です。

3．世田谷区梅丘　羽根つき餃子

自宅でも簡単につくれる羽根つき餃子として紹介されていたのは，世田谷区梅丘にある羽根木餃子というお店の羽根木餃子（5個入り600円）で，箱（小130円，大150円）に入れてくれるのでプレゼント用にもよいとのこと。
このお店ではトマト餃子，パクチー餃子といったものも販売しており，お取り寄せは半年待ちだそうです（！）。

4．鹿児島　1個で3種類の食感が楽しめるご当地餃子

丸っこくてお菓子のような見た目の，かごしま黒豚餅餃子（14個入り1080円）。
火が通った状態で送られてくるので，オーブントースターで焦げ目がつくまで焼くだけ，つまり誰でも焼ける餃子なんだとか。
スープに入れて水餃子，揚げて揚げ餅と，いろんな楽しみ方が紹介されていました。

5．香川県　2年待ちのクロワッサン餃子

たれ屋というお店のクロワッサン餃子（生餃子）（50個入り2650円）。
14年かけて完成した，すべてに徹底的にこだわりぬいた餃子らしいです。
2年待ちとはいえ，食べてみたい……。

番外編，手軽に作れる絶品たれ

番組内で餃子の裏技として紹介されていた，たれのつくり方についても備忘をかねて以下にまとめます。

• お酢に胡椒のみのシンプルなたれ

• 砂糖・はちみつ・酢・白味噌を混ぜるだけのお手軽な味噌だれ

2016年の餃子は「もちもち」らしい

完全なる私得でしかない備忘録です。

スーパーJチャンネルの特集，達人おススメ！「進化する餃子」④より。
2016年の餃子は「カリカリ」から「もちもち」へ，とのことです。

餃子の達人として，餃子評論家こと橘田いずみさんが「もちもち」をコンセプトに紹介されていたお店の情報をまとめておきます。
橘田いずみさんは個人的に応援しています。

橘田いずみのザ・餃子

渋谷区恵比寿の花林糖餃子

醤油（540円），味噌（650円），パクチー（650円）と3種のスープにもちもちの皮をした餃子の入った「スープ餃子」，橘田さんは味噌がお好きだそうで。
もっちりとした皮だと，うどんのようにしっかりと食べ応えがありスープの味もしみる，そんな一品になりそうでおいしそうですね。
tabelog.com

新宿区飯田橋のPAIRON

台湾の水餃子をさらに進化させて，中にいろいろなものを入れた色鮮やかな餃子が目を惹くお店でした。
番組で橘田さんとレポーターの方が召し上がっていたのは白龍餃子，青龍餃子，黒龍餃子，レッドドラゴンという4種類で各390円。
個人的にはシナモンの入った白龍餃子が気になります。
tabelog.com

世田谷区成城学園前のギョーザブラザーズ

ミニ中華まんのようにも見える（放送では「白玉？」なんて言われていました）かわいい餃子，もちもちまるチーズ／キムチ餃子（421円）が紹介されていました。
白いチーズ入りの方の餃子は求肥のようだと，言っているのを聞いてうまい表現だなと感心した次第です。
定番の焼き餃子もカリカリ・ジューシーとのことで食べてみたいです。
tabelog.com

豊島区池袋の中国家庭料理楊2号店

四川料理のお店とのことで辛いものが食べられない私には縁がないかもしれない，と思っていましたが，水餃子が絶品らしく非常に気になります！
紹介されていたものは真っ赤なスープが見ているだけで痛い紅油水餃子（920円）でした。
食べログのレビューで知ったのですが，孤独のグルメにも登場したお店なんですね。
tabelog.com

もちもちとした皮の餃子，とても気になりますね。
ごはんと一緒に餃子を食べるのと同じよう感覚になる部分もあるのかもしれません。
餃子と言えばカリカリ派でしたが，年も明けたことですし新たな境地を開拓してみようかしら。

この記事を書くために深夜3時に特集を見直していてしみじみ思いましたが，夜中に見るものじゃないですね，本当に。

mvrnorm関数で「'Sigma' is not positive definite」というエラーが出たときの対処法

シミュレーションなどのために，多変量正規分布からのサンプリングを行う関数として{MASS}パッケージのmvrnorm関数があります。
この関数は次のようなプログラム例を通してサンプリングを行います。

library(MASS)

#3変量正規分布の例
#平均値として0が3つのベクトルを作る
mux=numeric(3)

#3*3でそれぞれ0.5の相関係数をとる相関行列
sigma<-matrix(0.5, 3, 3)
diag(sigma)<-1

#サンプルサイズ=100，3変数で多変量正規分布からのサンプリング
data<-mvrnorm(n=100, mu=mux, Sigma=sigma, empirical=FALSE)

empiricalというのは与えた相関行列を寸分たがわず再現するかしないか，というロジカルな引数で，TRUEだとcor(data)として相関係数を求めたときに与えた相関行列と完全に一致するため母相関行列のように使えるかと思います。
今回はFALSEとしているので，標本から推定した相関係数のように使うことができます。

さて，ここから本題ということで次のような状況を考えてみましょう。
x1-x10までの10個の変数があるときに，x1-x9までは相互に0.1で相関しており，x10だけは他の変数と0.5で相関しているという状況です。
これを前提として，10変量正規分布からのサンプリングを行うプログラムを次のように書きます。

library(MASS)

mux<-numeric(10)
sigma<-matrix(0.1, 10, 10)
sigma[10,]=0.5
sigma[,10]=0.5
diag(sigma)<-1

data<-mvrnorm(n=100, mu=mux, Sigma=sigma, empirical=FALSE)

そうすると，
Error in mvrnorm(n = 100, mu = mux, Sigma = sigma, empirical = FALSE) :
'Sigma' is not positive definite
というエラーメッセージが返されて，データが生成されません。

この原因の確認と対処法をいかに述べようと思います。

エラーの原因の特定

メッセージを読めばわかるので言うまでもないのですけれど，相関行列の設定が不適切ということになります。
ただ，それで終わっては芸がないので，もう少し細かく見ていきます。

まず，mvrnorm関数の中身は次のようになっていることを確認します。
見やすさを重視して，スクリーンショットにしてみました。

この中で今回関心があるのは，
stop("'Sigma' is not positive definite")
という条件式までの部分になります。

かつ相関行列の指定までですので，次の部分が問題となっていると考えてよいでしょう。

eS <- eigen(Sigma, symmetric = TRUE)
ev <- eS$values
if (!all(ev >= -tol * abs(ev[1L])))
stop("'Sigma' is not positive definite")

これを，日本語に翻訳してコメントをつけてみます。

#与えた相関行列／共分散行列の固有値を求める
eS <- eigen(Sigma, symmetric = TRUE)

#固有値の計算結果（固有値，固有ベクトルが含まれる）の中から固有値だけを取り出す
ev <- eS$values

#もしすべての固有値が条件式を満たさなかった場合にはエラーメッセージを返す
if (!all(ev >= -tol * abs(ev[1L])))
stop("'Sigma' is not positive definite")

これらの内容から，今回のエラーは
if (!all(ev >= -tol * abs(ev[1L])))
という固有値に関する条件式を満たせなかったことが原因であることがわかります。

この条件式についてもうすこし細かく見てみます。
tolというのは上のスクリーンショットを見ていただければわかりますが，1e-06*1が割り当てられています。
abs関数は絶対値を返すものなので，abs(ev[1L])は1つ目の固有値の絶対値という意味になります。
また，!は否定演算子です。
したがって，この条件式は(-0.000001)*一つ目の固有値の絶対値を，得られた10個の固有値のすべてが越えなければFALSEを返すif条件式ということになります。

ためしに今回の状況における固有値の計算結果を見てみましょう。

> eS <- eigen(sigma, symmetric = TRUE)
> ev <- eS$values
> ev
[1] 2.9524175 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 -0.1524175

> tol=1e-06
> -tol * abs(ev[1L])
[1] -2.952417e-06

> ev[10]>=-tol * abs(ev[1L])
[1] FALSE

ひとつでもFALSEを返す場合にはエラーを返すようになっているので，かくして'Sigma' is not positive definiteというメッセージが表示されるわけです。

対策（力技）

さて，対策なのですが，要は与えた相関行列の組み合わせによっては固有値を計算する際に不適切なものが生じる恐れがある，ということなので，与える相関行列を最初から修正してしまえばいいじゃないかということになります。
そこで使うのが，{sfsmisc}パッケージのnearcor関数です。
nearcor {sfsmisc} | inside-R | A Community Site for R

これは不適切な相関行列を，cor関数でとりうる最も近い値にスムージングする関数とのことなので，素直に先ほどの相関行列をnearcor()の中に入れてからmvrnorm関数の相関行列として指定してあげれば，とりあえずエラーは消えるという力技的な解決法となっています。

お粗末様でした。

*1:=0.000001

加工肉は発がんリスクを上げるのか？

WHOが「加工肉には発がん性がある」といった旨の声明を出したことが世間で話題となっているようですね。www.bbc.com

日本人への影響は小さいだの，続報がどんどん出てきているようですが，個人的にこの声明自体どこまで本当なのかしらと疑問に思ったので元ネタを探してみようと思いました。
おそらくそうした声明を出すに至った研究結果があったのだろうと考えてちらっと探してみたところ，いくつかヒットしました。

その中でも，今回はCross, Leitzmann, Gail, Hollenbeck, Schatzkin, & Sinha (2007) という論文*1をつまみぐいしてみようかと思います。journals.plos.org

タイトルを眺めるに，赤身肉・加工肉の摂取と発がんリスクの関係に関する前向き調査といったところでしょうか。
簡単にですが概要を訳してみました。

論文の概要

Abstract

Background

赤身肉と加工肉はいくつかの解剖学的観点から発がんと関連付けられてきたが，しかしながら悪性の範囲と肉の摂取との関係についての前向き調査はなされてこなかった。
本論文では，赤身肉ないしは加工肉の摂取ががんのリスクを高めるかどうか，様々な観点から調査した。

Methods and Findings

• 1995-1996をベースラインとした，50～71歳の500,000人に対するコホート調査

• 食肉の摂取はベースラインで質問紙調査を行い推定

• Cox比例ハザード回帰で，赤身肉と加工肉を摂取している人の5分位におけるハザード比と95%信頼区間を推定

• 8.2年のフォローアップで，53,396件のがんが確認された

• 赤身肉摂取の最も高い20%と最も低い20%との間で，統計的に有意なリスク増加（20%～60%のレンジ）が食道・結直腸・肝臓・肺について見られた

• 加えて，加工肉を摂取する人の上位20%は，結直腸がんのリスクが20%，肺がんのリスクが16%高まることが示された

Conclusions

赤身肉と加工肉の両方の摂取が，結直腸と肺のがんとの間に正の関連が見られた。
加えて，赤身肉の摂取は食道と肝臓のがんのリスク増加に関連が見られた。

今回のお話し

事前のお勉強

まず，そもそも私の不勉強さの問題ではあるのですが，Cox比例ハザードモデルについてはあまり知りませんでした。
医療系で使うらしい手法だなーとか，そのくらいの印象です。
なので，これについて簡単に，手元にあった森田（2014）の『実証分析入門　データから「因果関係」を読み解く作法』の第15章で勉強してみようと思います。

実証分析入門 データから「因果関係」を読み解く作法

サバイバル分析

Cox比例ハザードモデルというのは，サバイバル分析の文脈で用いられるモデルだそうです。
サバイバル分析（Survival Analysis）*2について，これはある個体が死亡してしまった場合そこから先は観測されなくなってサンプルから落ちてしまう，という特徴を持ったデータを扱うときの手法の枠組みです。
サバイバル分析に関する詳細な説明はここでは割愛しますが，重要なのはこの枠組みではハザード率（Hazard Ratio; HR）というものを分析対象とすることだと言えます。
ハザード率は，森田（2014）では次のように定式化されています。

たとえば，a歳のある病気に罹患した患者を集めて30年間の継続調査を行ったとしましょう。
この病気が5年以内死亡率が非常に高いものだったとしたとき，分子や分母にあたるある時点での死亡者数や生存者数（確率）だけで考えると，25年目まで生きている人自体がそもそも非常に少ないためその時点での死亡率などを考えることに実質的な意味がなくなってしまいます。
そこで，ハザード比というものに変換する，要するに割合に直すことで刑事変化の問題を回避して本当の死亡率に焦点化するのがサバイバル分析で扱う対称の特徴と言えるかもしれません。

Cox比例ハザードモデル（Cox Proportion Hazard model）って？

さて，本題のCox比例ハザードモデルですが，このモデルは次のように定式化されます。

ここではベースラインハザードという，全個体で共通のものです。
また，右辺の指数の中には時間を表す添え字tが存在していないことが見てわかると思います。
つまり，ベースラインハザードという共通のものと，時点で不変である様々な変数=要因を含んだ部分の積としてハザード率が表現されているということになります。
比例ハザードという名前の意味ですが，今回のお勉強本の記述がわかりやすかったので（一部変更はありますが）そのまま引用します。

各個体のハザード率が共通要員であるベースラインの何倍か（倍率は説明変数で決まる），というベースラインハザードに比例する形で決まるので，比例ハザードと呼ばれる

このモデルは時間部分，すなわちベースラインハザードは推定せず，どの要因の影響が大きいかという観点に絞った分析を行うことが特徴的であると言えます。
ベースライン部分の情報として，存続期間を順序情報のみ用いたうえでパラメトリックな解析を行うことから，セミパラメトリックな手法であると考えることもできるようですね。
ハザード率の経時変化を考えずに要因の影響という関心のある部分のみをモデリングするという便利な手法ではありますが，欠点として次の3つがあげられていました。

1. 一部使わない情報がありもったいない
2. 比例ハザードという仮定自体が強いものである
3. 複数のリスク要因を独立とする仮定も強いものである

まとめると，ハザード率が全員共通のベースラインから，説明変数によって何倍になるかをセミパラメトリックな形でモデリングしたサバイバル分析の一手法ということですかね。

論文の図表を簡単に眺めてみる

Cox比例ハザードモデルをよく知らなかったのでいい勉強の機会だと思って記事を書き始めた，という事情がありここまででほとんど目的は達成してしまったのですが，せっかくなので少し今回取り上げた論文の結果を眺めてみようと思います。

たとえば，加工肉に関する分析を行ったTable 3ですが，基準となるQ1，すなわち加工肉を食べる頻度でソートした時の下位20%の人に対して，上位20%の人であるQ5の結直腸がんに関するハザード比は1.20でした。
表の脚注によると，説明変数としては年齢，性別，教育歴，結婚歴，がん経験，人種，BMI，喫煙，運動習慣，総摂取エネルギー量，飲酒，果物・野菜消費量を投入しているとのことです。

この論文ではがん発症をイベントとしているようですから，それまでにがんにかかっていなかった人がその時点でがんになる確率がハザード比としたとき，Q1の人に対してQ5の人は1.2倍（20%）多く発症リスクがあるということになるのでしょう。
ただ，信頼区間を見てみると1.07-1.30となっていますから，ほぼ変わらないかもしれないし，大きくて30%ほどなのかもしれません。
要旨で触れられていた肺がんについても，Q1に対してQ5のハザード比が1.16とのことですから16%とのことでしたが，信頼区間を見てみると1.06-1.26ということで個人的にはなんとも言い難いように思います。

検定結果は統計的に有意であったかもしれませんが，有意であることを以て真実であると主張できるわけではないですし，メタ分析などを通じて効果の大きさのより厳密な評価は必要かもしれませんね。
まだまだCox比例ハザードモデルへの理解が浅く，使い方も含めわかっていない部分が多いことから手続き自体に関する批判的検討ができないのが悔やまれますが，私個人の印象としてはそこまで気にしなくても大丈夫ではないかと感じました。
フォロー期間や説明変数を見ていて，赤身肉や加工肉の摂取が因果関係として発がんリスクを高めるという主張は，この結果からはまだ強すぎるようにも思われます。
もちろん論文では強い表現はつかっていないようですから，報道の仕方の問題もあるのかしら。

メタ分析をしている論文もあるようですから，後日気が向いたらabstractだけでも眺めてみようかと思います。journals.plos.org

*1:PLoS Med, 4

*2:イベントヒストリー分析（Event History Analysis），存続モデル（Duration model）とも

【R】任意の因子負荷から架空データを発生させる（1因子モデル編）

久しぶりの更新です。
更新したいネタはたくさんあるのですが，なかなか筆を執ることができません。

そんなどうでもいい個人的事情はおいておきまして，今回は備忘もかねて，統計ソフトRの{psych}パッケージを利用して任意の因子負荷行列から架空データを発生させる方法を書いておこうと思います。
{psych}パッケージのマニュアルはこちら。（※リンク先pdf注意）

大まかに次のような手順を踏みます。

1. 因子負荷行列を用意する
2. {psych}パッケージのsim関数を用いて相関行列に変換
3. {MASS}パッケージのmvrnorm関数で正規乱数を発生

とりあえずということで，まずは1因子モデルでの方法をこの記事にまとめます。
まとめるといっても手順は少ないので，ごく簡単に使う関数を載せて，簡単な具体例を付記するのみですが。
ちょっと手抜きをしており，各関数の引数は必要最小限のみ指定しています。

わかりやすさを重視して各セクションごとにコードを分割して記載，一番最後にコピーペーストや参照しやすさの観点からひとつにまとめたものを記載というスタイルといたします。

因子負荷の状況設定～相関行列への変換

まずはここからです。
せっかくなので，CiNiiでオープンアクセスできる心理尺度を参考に状況設定してみます。
今回は（自分の性格を鑑みて，ダメージ的な意味で）心に突き刺さるものがあったので教育心理学研究に投稿された次の論文を選んでみました。

小浜　駿（2010）．先延ばし意識特性尺度の作成と信頼性および妥当性の検証．教育心理学研究，58，325-337．ci.nii.ac.jp

この論文で作成された尺度の第1因子である「先延ばし中の否定的感情」を構成するとされる7項目について因子負荷を引用すると，.88, .82, .79, .77, .59, .57, .56という数字が得られます。
これを因子負荷ベクトルとして与えて相関行列を生成し，そこからデータを発生させるというわけです。*1

#psychパッケージの読み込み
library(psych)

#引っ張ってきた因子負荷をベクトルとして変数flの中に格納
fl<-c(0.88, 0.82, 0.79, 0.77, 0.59, 0.57, 0.56)

#psychパッケージのsim関数を使う
#今回は1因子モデルなのでsim.congeneric()を用いる
#生成した相関行列を変数cor.matの中に格納
cor.mat<-sim.congeneric(loads=fl)

任意の平均，相関行列から架空データ生成

ここまでで相関行列の生成ができましたので，このセクションではそこからデータを生成します。
{MASS}パッケージのmvrnorm関数で一発です，便利な時代になりました。

#MASSパッケージの読み込み
library(MASS)

#平均値のベクトルを用意する
#今回は7項目なので，各項目について計7つ平均があればよい
#今回は5件法とのことなので便宜的に3としました
mux<-c(3, 3, 3, 3, 3, 3, 3)

#架空データの発生
#作成したデータをデータフレームとするためにdata.frame()をかませています
#empiricalはTRUEならば指示通り（与えた通り）のデータを生成する模様
#nはサンプルサイズなので適当に与えてください
data<-data.frame(mvrnorm(n=100, mu=mux, Sigma=fl, empirical=TRUE))

これで，サンプルサイズが100の，7項目に関するデータが生成されているはずです。

コードのまとめ

2つに分割して記載しているので，ひとまとめにしたものも以下に用意しておきます。
説明のためのコメントが多いのは備忘録的性質の記事のためご愛敬ということで。

#パッケージの読み込み
library(MASS)
library(psych)

#引っ張ってきた因子負荷をベクトルとして変数flの中に格納
fl<-c(0.88, 0.82, 0.79, 0.77, 0.59, 0.57, 0.56)

#psychパッケージのsim関数を使う
#今回は1因子モデルなのでsim.congeneric()を用いる
#生成した相関行列を変数cor.matの中に格納
cor.mat<-sim.congeneric(loads=fl)

#平均値のベクトルを用意する
#今回は7項目なので，各項目について計7つ平均があればよい
#今回は5件法とのことなので便宜的に3としました
mux<-c(3, 3, 3, 3, 3, 3, 3)

#架空データの発生
#作成したデータをデータフレームとするためにdata.frame()をかませています
#empiricalはTRUEならば指示通り（与えた通り）のデータを生成する模様
#nはサンプルサイズなので適当に与えてください
data<-data.frame(mvrnorm(n=100, mu=mux, Sigma=fl, empirical=TRUE))

そのうち1～4因子まですべてまとめて因子負荷行列として与えたうえで架空データを生成できるようなものをまとめようという気持ちはあります。

論文で明記されているデータを再現するようなデータを発生させて，それを分析してみるなんていうのは勉強になりそうなのでいいですよね。

この記事では既存の関数だけを使ってやりくりしていますが，そのうち目的に応じて関数を組む必要も出てきそうです。

次の本，『みんなのR』と言いますが，バイブルクラスに出来がいいらしいので私も早く買います。
友人も即買って，べた褒めいていました。
なんだかんだ，Rもまだまだできない方なのが悩みの種です。

みんなのR -データ分析と統計解析の新しい教科書-

*1:本当は第4因子まであるので，4因子のモデルとして因子負荷行列を与えてデータ生成するほうが望ましいのだと思いますが，今回は単純化のために中から第1因子のみをピックアップして記事にしています。

【書評】靜哲人 (2007). 基礎から深く理解するラッシュモデリング 項目応答理論とは似て非なる測定のパラダイム.

IRT（項目反応理論）についての書籍は多数ありますが，Raschモデルに特化した書籍というのは見たことがなかったので，買って読んでみました（買ったのはだいぶ前なので，積本を消化したかたちになりますが……）。

基礎から深く理解するラッシュモデリング―項目応答理論とは似て非なる測定のパラダイム

本書の章立ては次のようになっています。

はじめに

第1部　数学的準備

第1章 　加算記号
第2章 　確率
第3章 　指数と対数
第4章 　微分

第2部　テスト理論の基礎

第5章 　平均・分散・共分散
第6章 　古典的信頼性

第3部　ラッシュ測定理論

第7章 　ラッシュモデルの導出
第8章 　ラッシュモデルの姿
第9章 　受験者能力と項目難度の分離
第10章　受験者能力と項目難度の推定
第11章　情報・誤差・信頼性
第12章　データのモデルへの適合度
第13章　多カテゴリーの項目
第14章　ラッシュモデルと項目応答理論

ラッシュモデルは多くのIRT関連本では1PLモデル（1母数ロジスティックモデル）として位置づけられ，項目反応理論における2PLなどの下位モデルであるとされています。
しかしながら，サーストン系IRTモデルとRaschモデルとでは数学的表現は類似しているものの成立した経緯が大きく異なり，したがってモデルとデータの関係について取る立場も大きく異なります。*1

そうした状況の中で，本書はRaschモデルが「項目反応理論とは似て非なる測定のパラダイム」であることを主張し，わかりやすく説明をしている一冊であるといえましょう。

Raschモデルは言語テスト系の方々からは今もアツい支持を得ているようで，本書の著者である靜先生も英語関係の学問領域の方であるようです。

Raschモデルの特徴は「客観測定（objective measurement）」という言葉で表現されることが多いです。
というのも，IRTではデータをモデルにフィットさせるように（とは本書の言葉）識別力や当て推量などのパラメータを追加し，データをよりよく説明できるモデルを利用するという意味でデータ本位ですが，一方でRaschモデルではモデルに適合しないデータが得られた場合は項目などをより鍛えるべきだ，のようにモデル本位の考え方をします。

その意味で，従来のモデリングが行う「データにモデルを合わせる」という考え方とは真逆の「モデルにデータを合わせる」という考え方は測定のパラダイムシフトである，というのが副題の示す事柄であるわけです。

書評と銘打っているので個人的な感想をいくつか。

まずはモデルの導出についてですが，非常にわかりやすく丁寧です。
具体例をふんだんに用いて説明がなされており，式展開も省略がほとんどないので（一部好みの問題で納得できない式変形はありますが）文系人間であっても問題なく追えるように思います。
不要な方は読み飛ばせばいいと思いますが，第一部として数章を割いて数学的準備をしてくれています必要な方にとってはありがたいですね。
統計的な予備知識のある方は第6章までは読み飛ばして，第7章から読まれても問題はないように思います。
ただ，古典的信頼性の話などはわかりやすく，自分としても再確認したものがありましたので既習の方にも一読の価値はあるかと思いました。

モデルの数式的な導出だけでなく，被験者能力や項目母数の推定についても数値例を豊富に載せて，わかりやすく書かれていました。
エクセルで計算例を多く提示している点も，イメージがしやすくなるためわかりやすさに寄与していると言えるでしょう。
エクセル部分は個人的には天丼な印象だったので，最初はきちんと追いかけたものの2回目からはざーっと流してしまいましたが。

Raschモデルの入門書として，またテスト理論の復習としてよい一冊だと思いますが，そんな本の中に多カテゴリーについてもきちんと触れてあるのは好感が持てました。

この本の一番のキモともいうべきは最終章である第14章の，「ラッシュモデルと項目応答理論」という部分なのだと思いますが，ここは私としてはあまりピンときませんでした。
著者のアツい思いは十分に伝わってくるのですが，「客観測定」という言葉やデータとモデルの関係に対する考え方の違いは説明されれば，「まあそうかな」という感じなのですが，いかんせん主張が強すぎるように感じられて一歩引いてしまいました。
私自身も本書を手に取ったときには，RaschとIRTは歴史的経緯が違うのだから別物なはずだと思っていたはずなのですが，読んでいたときの心情としてはそこまでアツく主張するほどのものなのだろうかという疑問や気持ち悪さのほうが大きかったように思います*2。

特に，個人的にわかりにくかったのは最後の章における項目識別力（本書では弁別力と呼んでいます）が全項目で等しいという制約について述べている部分です。
ここは後で考えてみると，得られるデータはガットマンパタンが望ましく，したがってレベルの高い／低い受験者が困難度の低い／高い項目について期待される反応でないものをとることは望ましくない。
言い換えれば，項目同士の関係が受験者の能力で変わるようなことはあってはならない*3ので，項目特性曲線が交わるのは都合が悪い，したがって平行＝全項目で識別力が等しいという制約が必要になる。
こんなロジックなんだと書評を書きながら気づきました。
この部分の説明については，積読をしていたからかはわかりませんが，最終章が言葉足らずであるように思われます。

また，索引がないことについては個人的に不満が残りました。

さて，本書の主張と概ね同意見であり，RaschモデルとIRTとは異なるものであると思っているはずなのにこのモヤモヤは何であろうか，と疑問に思ったので手元のIRT関連本をのぞいてみて，すこしだけすっきりしたのでそれについても簡単に述べておきます。

次に述べる事柄については，加藤・山田・川端 (2014) の『Rによる項目反応理論』を引用・参照しました。

Rによる項目反応理論

• Raschモデルの考え方は「モデルは常に正しい。その正しいモデルに合うデータを持ってくることが肝要だ」

つまり，測定における理想をモデル化したものであり，現実の測定も理想であるRaschモデルに近づけるべきであるという主張です。
そのために項目を選定・吟味し，モデル制約を満たす測定を行っていこうというわけです。
guessingやslipといったパラメータはデータをうまく説明するために恣意的に付け足されたものである，と靜（2007）*4は述べます。

• Raschモデルは「不変性（invariance）」を強く強調する立場

不変性とは古典的テスト理論と項目反応理論の差異ともいうべき特徴で，次の2点を満たすことを言います。

1. 項目パラメタが，テストを受けた受験者集団に依存せずに推定される
2. 能力パラメタが，テストに用いた項目に依存せずに推定される

すなわち，項目や能力が相互に依存せず，分離されていることを指します。

これを以てRaschモデルでは，モデルにフィットしていれば測定の客観性が担保される，という主張をしていると靜（2007）は述べます。
しかし，加藤・山田・川端2（2014）が引用している村木（2011）は，IRTにおける識別力パラメータを例に，それらに標本の影響が全くないわけではないことを述べて，パラメータの普遍性あるいは不変性を強調しすぎることには警鐘を鳴らしています。
Raschモデルの場合にはモデルへのデータのフィットを考えることでこの問題に一定の意思表示をしているものとは思われますが。

モヤモヤを解決するために脇道にそれましたが，閑話休題，書評をまとめたいと思います。
本書に関する個人的な感想を次のように箇条書きしてみました。

• モデルの導出が丁寧である

• 推定の説明もわかりやすく，エクセル計算の提示によりイメージがしやすい

• 索引，お願いします

• 最終章ではアツい主張がみられる（すこし言葉足らずな印象を受けました）

• 入門書としては，数学的準備も完備されており良い

以上，僭越ながら書評でした。

*1:どこかのタイミングで，両モデルの成立というか導出というかについては軽くまとめたいと考えています。

*2:もちろん個人の感想であり，特定の人物・領域などに向けた悪口ではありません。

*3:2PL以上のIRTでは識別力が各項目で異なるため，項目1と項目2の難しさが受験者のレベルによって異なるということが許容される

*4:書評を書いている「本書」にあたります。加藤・山田・川端（2014）の引用部での記述時に紛らわしさがないようこの表現を用いました。

【書評】涌井良幸 (2009). 道具としてのベイズ統計.

わかりやすさに定評のある涌井先生の，ベイズ入門本です。
途中まで読んで積んでいたものをようやく読み終えたので，今更ですが書評など。
章立ては次の通り。

序章：GoogleもMSもベイズ統計！
1章：ベイズ統計の準備をしよう
2章：ベイズの定理とその応用
3章：ベイズ統計学の基本
4章：ベイズ統計学の応用
5章：MCMC法で解くベイズ統計
6章：階層ベイズもExcelで
付録

この本の優れた点を次のように箇条書きしてみました。

• ベースとなるベイズの定理に関する例題が豊富
• ベイズ推定と最尤推定とを対比して学習できる
• 共役事前分布に関する説明が丁寧
• 昨今普及したMCMC法に関する説明が圧倒的にわかりやすい
• 各章にベイズ推定を行う例題がついており，イメージがわきやすい（特に経験ベイズ，階層ベイズはわかりやすいです）
• （Excelでの計算ワークシートを組めるよう丁寧に解説している）*1

本書は伝統的な頻度論的統計学の既習者で，かつベイズ統計に関するスキーマを作りたいという人にはうってつけの一冊だと思います。
この一冊をとっかかりにしてイメージを作った後に，久保先生の緑本や数ある入門書（渡辺先生，松原先生など）を読むと学習がスムーズになるかもしれません。

この本の要旨として，次のようなことをおさえておくとベイズ統計を理解するための枠組みを構成できるかもしれません。

• ベイズの定理とベイズ統計学の関係
• 最尤推定法とベイズ推定との関係
• ベイズ推定を行うための2つの方法
• MCMC法（マルコフ連鎖モンテカルロ法）とは何か？

次からひとつずつ，部分的に所見を交えながらですが述べていきたいと思います。

ベイズの定理とベイズ統計学の関係

（中学や）高校の教科書にも掲載されているベイズの定理とは次のように定式化される定理です。

• ：事象Aが起こったことを所与としたときに事象Bが起こる確率

• ：事象Aが起こる確率

• ：事象Bが起こる確率

この定理より，事象Bが起こったというデータが観測されたときに，その下で事象Aが原因である確率を求めることができると理解されます。

これをすこしだけ統計学風にアレンジして，事象Aを仮説H，事象BをデータDという形で読み替えて表現すると，次のようになります。

たとえば，Kさんが使っている携帯電話がiPhoneである，というデータが得られたとします。
国内のシェアをA社が30％，D社が45％，S社が25％としたときに，KさんのキャリアがA社である確率を求める，ということを考えます。*2
このとき，どの携帯会社もiPhoneを販売していますから，仮説としては次の3つが考えられます。

1. ：Kさんが契約しているのはA社である
2. ：Kさんが契約しているのはD社である
3. ：Kさんが契約しているのはS社である

このとき，数式上の確率については次のような意味となります。

• ：キャリアがの条件の下で機種がiPhoneである確率
• ：キャリアがである確率
• ：携帯電話がiPhoneである確率

このうちには「尤度（liklihood）」，には「事前確率（prior probability）」という名前がついています。

こうしたことがベイズの定理について言えるかと思いますが，これらはあくまでベイズの「定理」についての話であり，ベイズ「統計学」の話ではありません。

ベイズ統計学を考えるにあたって，定理をより使いやすく書き直すと次のように表現できます。

この式は，（事後分布）が（尤度関数）と（事前分布）の積に比例することを意味します。*3
この式を利用してパラメータの事後分布を求めることがベイズ推定であると言えるでしょう。

最尤推定法とベイズ推定との関係

ベイズに興味を持つ方であればご存知とは思いますが，頻度論者はパラメータを所与のもの，すなわち「神のみぞ知る定数」であると考え，一方でベイジアンはパラメータが分布をもつ「確率変数」であると考えます。

根本的な態度こそ違いますが，最尤推定とベイズ推定は関連付けておくと理解が深まる印象を個人的に持っています。
というのも，無情報事前分布である一様分布を事前分布としてベイズ推定を行ったときに，事後分布のモード（事後モード）は最尤推定値と一致することが知られているためです。

誤解を恐れずに言うと，ベイズ推定は最尤推定の拡張である，と考えることができると思います。
上述したベイズ統計学での基本式を再掲し参照することで，それを確認します。

この式のうちは頻度論的な統計学が扱う尤度関数と同一のものであることから，頻度論的な統計学での最尤推定は次のようにも表現できます。

この式はベイズでいうところの事後分布は尤度関数と同じである（正確には定数倍である）ことを示しています。
したがって，ベイズ統計学というのは，最尤推定に事前分布を仮定すること・導入することを以てパラメータの事後分布を推定するものだと考えることができます。

ベイズの定理の側からも確認をしてみます。

このうち，はデータとして観測されており所与ですから，事前分布にあたる分子のを消してしまえば尤度関数を最大化するを推定するという最尤推定の問題に帰着しますよね。

ベイズ推定を行う2つの方法

1. 共役事前分布の利用
2. MCMC法の利用

このうち，MCMC法の利用が昨今爆発的に増えている印象です。
それを行うための環境*4が整備されてきているのが理由でしょうか。

MCMC法（マルコフ連鎖モンテカルロ法）とは何か？

ベイズ統計では確率分布を扱うため，積分計算が必須になります。
モデルが複雑になればなるほど，積分計算も難しくなり，臨終します。
そうした煩雑さを避けるために，周辺分布を利用したリサンプリングを重ねることで数値計算で力押しして疑似的に分布を求めてみよう，というのがモンテカルロ法の発想であると言えるかと思います。
その際，マルコフ性という性質を持たせることでリサンプリング効率を上げよう，というのがマルコフ連鎖モンテカルロ法の考え方といっていいでしょう。

ベイズ統計自体は今後も長く使われていくことが予想されますから，頻度論者であっても知識として知っておいて損はしないのでしょうね。
私個人としては今まで頻度論の中で統計を学んできたので，宗旨替えをするか否かを考えるためにももうすこしベイズについての勉強をしていこうと思っています。
頻度論的統計学とベイズ統計学では母数に対する考え方が違いますから，現象への態度の違い，モノの見方の違いなど，そうした数理的な面や実用性以外の面も併せて最終的な判断をしたいですよね。
頻度論的統計学の態度自体は非常に科学的であり，捨て去るべきものではないと考えます。
一方でベイズ統計の考え方は状況によってはより良く現象を記述できる可能性を有しており，量子論的な世界観との親和性は高いかもしれません。
頭を悩ませる夜は当分終わらなさそうです。

*1:これは個人的には重要視しませんが，ビジネスの世界では統計ソフトの導入やパッケージのインストールが個人の意のままというわけではないですから，大事なのだと思います。

*2:話がややこしくなるのでSIMフリーのことは考えません。

*3:文系人間としては，「パラメータの」事前／事後分布と呼ぶ方が丁寧でわかりやすかったです。

*4:BUGS, JAGS, Stanなど