「なにかやりたい気持ち」改め、
「説明変数にノイズを含むモデルで遊んでみた」



この記事は Stan Advent Calendar 2016 の12月10日の記事です。
Stan Advent Calendar 2016 盛り上がってますね!!楽しいです!!

アヒル本 StanとRでベイズ統計モデリング の7.7「説明変数にノイズを含む」のモデルで遊んでみた、という内容です。モデルについての誤解がある気がしています。お気づきの点がございましたら、ご指摘いただければ幸いです。
 
 

説明変数にノイズを含むモデル


このモデルは説明変数をX、目的変数をYとしたときに、以下のようになります。
\begin{eqnarray}
X[n] & \sim & Normal(X_{true}[n], 1.0)\\
Y[n] & \sim & Normal(\alpha+\beta X_{true} [n], \sigma)
\end{eqnarray}


はじめの式は「説明変数が平均を真の値(Xtrue) 、標準偏差を1.0とした正規分布から発生している」、つまり「データとして得られた説明変数には標準偏差1.0のノイズが含まれている」ということを表現しています。ここではノイズについて1.0という固定値を与えていますが、もちろん変更できます。
二つ目の式では、説明変数の真の値で目的変数を回帰しています。アルファが切片、ベータが回帰係数、シグマが誤差標準偏差を示します。


構成概念をあつかう心理学では、説明変数に測定誤差などのノイズが含まれることが多いです。そこで複数項目で測定を行い、構造方程式モデリングで誤差と分離して目的変数への影響を検討しているわけなのですが、既存尺度を利用してデータを得た際に因子分析とかを各研究がそれぞれやり直しているような状況です。先行研究の因子分析の結果からは「どの項目にどの程度因子が寄与していて、ノイズがどのくらい含まれるのか」ということがわかるわけですから、そういった情報を利用して「少ない項目数、少ないサンプルサイズでも妥当な推定をする」ということができるといいなぁと思っていました。因子分析は潜在変数が説明変数、観測変数が目的変数となっている回帰分析と捉えることができます。はじめの式は「研究者が見たい真の値に誤差がのっかった状態で説明変数が生成されている」ということを意味していますから、アヒル本でこのモデルをみたときには興奮して、鼻血がでるんじゃないかと思いました。「そんなうまい話、、」というあれもあったので、鼻血はでなかったわけなのですが。


 

仮想データを発生させて遊んでみる


ということでデータの発生から。先に記述したモデルどおりです。

これをlmで分析すると下の結果。真の値 (t) を使って回帰分析をすると、(当たり前ですが) 切片も回帰係数も設定した値 (それぞれ3と5) に近い値になります。説明変数を使って分析をすると、切片は3.04、回帰係数は2.33となりました。

 
 
はたして、説明変数にノイズが含まれていることを仮定したモデルだと、真の値で回帰したときの値に近い推定ができるのでしょうか?発生させた仮想データに対して、説明変数にノイズが含まれていることを表現した以下のモデルで推定します。モデルのコードはそのまんまですね。ここでは説明変数のノイズの大きさに、データ発生の際に設定した値である1を固定値として与えていますが、実際に使うときには、先行研究の結果に基づいて分布を与えることになるのではないかと思います。

 

推定結果はこちら。見事にノイズを含む説明変数をつかってlmで分析したときと同じ値が得られており、回帰係数については設定した値である5とはほど遠い値が得られました。やはりそんなうまい話はないのか。。


 
 
ちょっと普通にX1で回帰分析するモデルで推定してみます。回帰分析における誤差の部分 (e) が、説明変数にノイズがあることを仮定するモデルよりも大きくなっています (ただし、lmで分析したときの残差の標準偏差は3.61なので、不当に大きいとか、そういうことではないです)。「説明変数にノイズが含まれている」という仮定は回帰係数の推定にではなく、誤差の推定に効いてくるようです。

 

よろしい、ならば3項目だ


変数が1つだけというのはさすがに無茶だったのかなと思い、変数を増やしてみました。データの発生は変数の数が増えただけで、先ほどと同様です。また、真の値 (t) を使って回帰分析をすると、やはり切片も回帰係数も設定した値 (それぞれ3と5) に近い値になります。3つの変数を平均した値を説明変数としてもちいると回帰係数は3.81となりました。


さて、このデータを以下のモデルで推定します。3項目の平均値をつかうモデルでも推定し、比較を行います。


結果はこちら。やはり、両モデルの違いは、切片や回帰係数の値ではなく、誤差の部分にあらわれています。


なお、このデータを普通にSEMで分析してみると、設定した値5に近い値が得られました。

 
 

まとめ


心理学のデータでは説明変数にノイズが含まれていることが多く、「ノイズを除去してより正確な変数間の推定ができれば最高だ!」と興奮したわけですが、回帰係数の推定に効くという感じではなさそうです (仮想データの発生の仕方が悪いのかなぁ)。このあたり、自分はまだまだモデルについての感覚が甘く、遊んでみないとわからないという部分でした。先行研究の因子分析の結果をつかって、「少ない項目数、少ないサンプルサイズでも妥当な推定をする」ということができるといいなぁと思ってます。

 
 

【追記】事前分布の設定をするとうまくいく!?


2016.12.12追記
関西学院大学の清水裕士先生に事前分布の設定をするとうまくいくというご助言をいただきました。
モデルにはパラメタの事前分布についての設定を追加しています。特に、パラメタtについては「平均mu、標準偏差sigmaに従う正規分布から発生している」という仮定をいれており、muとsigmaを新たにパラメタとして推定しています。

 
 

推定結果はこちら。データ発生の際に設定した値はa=3, b=5, e=0.5, mu=0, sigma=1なのですが、(特にtについての) 事前分布を設定しない場合に比べ、今回のモデルではそれぞれ近い値が得られています。ただし、今回は反復回数を10万、間引き数を5に設定していて、4万のmcmcサンプルがあるのですが、有効サンプル数 (n_eff) がとても少ないです。もともと説明変数にノイズが含まれていることを表現したモデルだと反復回数や間引き数を多くしないと収束していなかったのですが、事前分布の設定を追加すると、さらに収束のために反復回数や間引き数が必要な感じです。

春の方法論セミナーに参加してきました


日本社会心理学会の第3回春の方法論セミナーに参加してきました。
第3回のテーマは「統計的因果推論への招待」で、過去回と同じく、説明がわかりやすくてとても楽しかった。
以下は感想。これから勉強。

統計的因果推論では、「AがBの原因である」ということの条件に、AがあるときのBの確率(p(B|A))とAがないときのBの確率(p(B|¬A))で、「p(B|A) > p(B|¬A)」ということがある (ようだ)。この条件が必要条件と考えられているのか、十分条件と考えられているのかはわからない。

でも、「p(B|A) > p(B|¬A)」だけだと、「AがBの原因である」場合と「BがAの原因である」場合の区別がつけられない。なので、「p(B|A) > p(B|¬A)」は「AがBの原因である」ことの必要条件であって、十分条件ではない 。ということで、「AがBの原因である」というためには、「p(B|A) > p(B|¬A)」とは別の基準を満たす必要がある。たとえば、時間的先行性で、AがBよりも時間的に先行している必要がある、とか。あとは、AとBの共通原因の影響を排除できているかどうか、とか、だろうか。とりあえず、僕が持ってる素朴な因果概念はこの3つが満たされていたらOKだと思う。

個人的に目から鱗だったのが、星野先生のトークで出てきた「割り当てられなかった条件のデータは欠測データになる」ということ。

大塚先生のトークにあったように、「(Aさんは) 予防接種したからインフルにかからなかった」は「予防接種してインフルにかからなかった」と「予防接種しなくてインフルにかかった」を観察しないといけない。けど、現実ではどちらかしか観察できない。個別事象の因果についてはその検証をするために可能世界を想定しないといけない。心理学のように、個人レベルではなく個体群レベルで考えるなら可能世界を想定する必要はないのかもしれないけど、個体群レベルで考えてもやはり、割り当てられなかった条件のデータは (観察したいけど) 観察できず欠測データになる、というのは同じ (だろう)。

実験で無作為割り付けが重要なのは、実験操作のみが群間で異なるようにするため。だからこそ、従属変数の差を操作に帰属できる。もちろん、これは当たり前のことなんだけど、知識としてもっていただけというかなんというか、まぁ、「無作為割り付けしてるからいいじゃん」みたいな感覚だった。無作為割り付けは交絡要因を統制する消極的な方法、これも知識としてもってる。でも、「本当に実験操作のみが群間で異なるのですか?」といわれたら「せやかて!」と思うレベルには正当化している。

前置きが長くなったけど、要は、「割り当てられなかった条件のデータは欠測データになる」と聞いた時に、観測データを欠測セルに代入している操作が頭に浮かんできて (それ自体は間違ったイメージだろうけど)「うげっ」となった、という話です。傾向スコア万歳。積極的に共変量を除去していこう。

とはいっても、共変量を統計的に除去する、というのには昔から何か引っかかるものがあって、もうちょっと自分の考えを整理しないといけないなぁと思っているところ。でも、たぶん、実際の使用に引っかかってるだけのような気がする。共変量として何を入れるべきかということについて、従属変数と関連が強いが欠測と関連が低い変数を、といった基準もあるのだろうけど、概念的に共変量として考えられるか、ということも重要だろう。媒介変数として機能しているものを共変量として扱っちゃったらまずいんでしょう?「共変量」って言葉の (個人的) イメージが悪いのかもしれない。あと、共変量をめちゃくちゃ除去するという方法は、因果関係があるかないかの判断には使えると思うけど、そこで推定される値に心理学的な意味はあるのだろうか、という疑問。

あと思ったことは、
1) 測定の信頼性についてちゃんと勉強しよう、
2) 林先生のトークであった「理が知りたい」と「介入がしたい」の違いがわかってない、
3) 講演概要にあった「因果推論と統計推論とは、そもそも世界と認識についての考え方 (パラダイム) が根本的に異なっている」というのがわかってない、
の3本。

昨年のGLMMに引き続き、わかりやすく説明してもらったので、関連本を読める気がしてきた。
ありがとうございました。

広島ベイズ塾で発表しました


魁!!広島ベイズ塾でベータ分布について発表しました。今回の会合では「分布感をつかむ」といういまいちつかみきれないテーマで、いくつかの分布についてお勉強しました。正規分布、ポアソン分布、ベータ分布、ガンマ分布、コーシー分布、負の二項分布です。ガンマ分布、コーシー分布については、コワイ本読書会のなかでもよく事前分布としてつかわれていたんですが、どんな分布かまったく把握していなかったのでとても勉強になりました。

負の二項分布についてはあらためて復習しないとわからないですね。「2回くらい変身します」っていわれたときのフリーザ感ったらなかったですね。これが分布感か。

分布の勉強とは離れて、ベイズについてもいろいろ教えてもらいました。なにやら変分ベイズ(?)やノンパラベイズ(?)がrstanで近々できるようになるそうです(ノンパラベイズはもうちょっとかかるんだったかな)。変分ベイズができればmixtureがMCMCできるらしい。ノンパラベイズができればmixtureのときの妥当なクラスの数も推定可能になるらしい。原理も教えてもらったのですが、まったくわかりませんでした。でも、期待で胸が膨らみます。

一応、ベータ分布についての資料をあげておきます。


ハムストリングスと腰痛

ハムストリングスはふとももの裏の筋肉の総称です。大腿二頭筋、半腱様筋、半膜様筋で構成されています。股関節と膝関節をつなぐ二関節筋で、基本的には膝関節屈曲動作、股関節伸展動作ではたらきます。

この柔軟性が低下すると、わたしの場合(おそらく一般的にも)、腰痛がでやすいです。ハムストリングスは股関節についていますので、ここが硬くなると上半身を前傾させるときに、骨盤からうごいてくれなくなります。そうすると、腰の部分を大きく曲げる必要があるため、負担がかかって痛みにつながるということらしいです。

ということで、腰痛が気になる人は、一度ハムストリングスの柔軟性をチェックしてみてはいかがでしょうか。もちろん、腰痛の原因は他にもたくさんありますので、ここの柔軟性にだけ注意すればOKというわけではありませんが。

もも裏のストレッチといえば前屈が一般的ですが、「手をどこまで伸ばせるか」ということを意識しすぎて、意外ともも裏が伸びていないということがよくあります。しっかり骨盤から倒してやって、もも裏を伸ばすことを意識する必要があります。しかしながら、「ハム伸ばさないとなー」と思うころには、ある程度ハムストリングスの柔軟性は低下しており、前屈の際にしっかり骨盤から倒すということがすでに難しくなっている、ということが少なくありません。場合によっては、無理に前屈をすることで腰の負担を高めてしまうかもしれません。

おすすめは、しゃがんだ状態で足首をつかんで膝を伸ばしていく、という方法です。このとき、胸とふとももを密着させておくことに注意して、膝を伸ばしていきます。普段、柔軟体操をしていない方は、膝を少し伸ばそうとするだけで、強烈なストレッチを体感できることと思います。あくまで「ハムストリングスを伸ばす」ことが目的ですので、膝は伸びきらなくてOKです。胸とふとももが離れない程度に膝を伸ばしてください。文字で書くとわかりづらいですね。ジャックナイフストレッチという呼び方があるようなので、詳しくは検索していただければと思います。

立位での前屈との違いですが、「前屈は膝関節を伸ばした状態で股関節を屈曲させる。先ほどの方法は股関節を屈曲させた状態で膝関節を伸ばす」という違いがあります。膝関節を伸ばす際には、大腿四頭筋(ふとももの前側) の力を使えますので、ハムストリングスを伸ばしやすいです。