Kosugitti's BLOG

アンドロイドは正規分布の夢を見るか

この画面は、簡易表示です

R

分位点回帰Quantile Regressionは確かに面白い

社会心理学会大会2013@沖縄国際大学,はそこそこ楽しめたのだけど,中でも一番面白かったのが分位点回帰についてのWS。

前々から企画者のIさんに面白さは伺っていたのだけど(論文も書かれてましたね),WSで確認し,先ほど実際自分で触ってみて面白さを味わった。当然,Rでできるのである。

少し宣伝?解説?しておくと,従来の回帰分析が平均回帰を狙っているのに対し,分位点回帰は任意の分位点(30パーセンタイル点とか,第3四分位とか)での回帰係数を求めるというもの。しかもひとつではなくて,複数求めることができるから,第1四分位ではこういう傾きだったのに,第3四分位では傾きが変わりましたね,なんて事も分かる。

心理学変数のほとんどは正規分布するという便利な建前があるが,実態的データ,例えばネットワークサイズとか年収のようなもの,は当然偏った分布をするのであって,平均点目指して回帰する時点ですでに歪んでしまっているわけである。それを補正する可能性があるのがこの手法。

何より感心したのは,WSで「我々の仮説が平均回帰に縛られすぎてないか?もっと自由であるべきではないか?」という問いかけがあったこと。確かにそうだよなあ,実際のデータとか経験則は,必ずしも平均のことばかりではないからなあ。

 

さてRではquantregパッケージというのがあって,そのなかのrq関数が分位点回帰をしてくれる。QuantileRegressionだからqr関数じゃないの,と思うかもしれないが,qr関数はすでに固有値分解のQR法に割り当てられた名称なのでダメなのです(これは譲れない)。

quantregに入っているエンゲル係数データを使って,試しにちょっとやってみた。
[crayon-5bc47a3109e3a894607371/]
コードはfoodexp(食費)をincome(収入)で予測するだけの式だけど,tau=seq(0,1,0.25)のオプションで0,0.25,0.50.0.75,1.0の分位点回帰をせよ,という意味になっている。0とか1は無意味だけど,0.25,0.75は四分位だし,0.50は中央値回帰になっている。結果は明らかで,収入が下位25%のところだと0.47しかない係数が,中央値で0.56,第3四分位で0.64とあがっていき,上位階層の人間ほど食費にお金が回せるのだな,ということが明らかだ。

お金がらみでついでにもうひとつ。野球の年俸と成績の関係を分位点回帰にかけてみた。データはここから持ってきたものを使っている。
[crayon-5bc47a3109e4a251450134/]
これをみると,年俸が低いうちはホームランや打率がものを言うが,年棒が高い層は打率はむしろ負のパスになって,ホームランや三振のパス係数がおおきくなっていく。つまり,記録より記憶で儲けるようになると,一流ってことですかな。

人文社会科学系のデータは歪んでいるのも一杯あるだろうから,こういうモデルで自由な仮説検証がもっと進むといいですね。

>> See also RPubs http://rpubs.com/kosugitti/10473

 

参考文献はこちら。安くていいと思う。

[amazonjs asin=”1412926289″ locale=”JP” title=”Quantile Regression (Quantitative Applications in the Social Sciences)”]



SEMから教える。SEMだけで教える。

Yamadai.RというRの勉強会を定期的に開いているのだが,参加者が基本的にゼミ生なので,ゼミ生に対する統計の補講のようになっている。まあそれでもいいんですけど。そういう位置づけの会になっているので,授業では触れないやり方とか,新しい教え方を試す場にさせてもらっています。

 

普段の統計の講義(データ処理法)では,検定,回帰分析,因子分析,とオーソドックスな順で教えていくのだけど,いつも最後の方は時間が足りなくて,SEMまで言及できなかったりする。ということで,Yamadai.Rでは逆に,最初からSEMを教え始めた。回帰分析は説明したんだけど,lavaanパッケージを使ってやるのである。

重回帰分析を教えるときに,「独立変数間は無相関でないとねー」とか注意しないといけないが,lavaanを使っていると「このモデルは独立変数間に相関が考えられるよね,パスを足してみようか」という対応ができるわけで,わざわざ遠回りした上で「昔は苦労した」みたいな話をしなくていい。

さらに別の教育効果があることを先日発見した。

データとして,プロ野球選手の年棒とシーズン成績をここからとってきて使ってみた。打者のデータに限定して回帰分析につかうと,例えば「ホームラン一本打ったら年棒がいくらあがるか」という話ができて,学生もイメージしやすい。いい教材を見つけたと思う。

で,前回は潜在変数を仮定したモデル(因子分析)を教えていた。「安打数」「HR数」「四球数」「体重」を観測変数として,【腕力因子】を仮定し,「盗塁数」「犠打数」「体重」を観測変数とした【脚力因子】を仮定したモデルを描いてみたのである。このモデルが妥当かどうかは検証していくとして,まずこういう目に見えないものを背後においてみる,というお話なので,やってみて「おおできた」,「腕力と脚力が負の相関,なるほどな」みたいな実感を得た。いいかんじ。

そのあと,モデル修正指標を見てみると,腕力因子から盗塁数にパスを出すと良い,という事であった。ここでちょっと「?」となったのだが,学生が「腕を振るから走るんじゃないか」「スライディングのときに腕をのばすのではないか」「まさか逆立ちして走るのか」とか言い出したので,シメタ!と思ったのである。

【我々は,背後に潜在変数があると仮定してモデルを描いた。その潜在変数が「腕力」であるかどうかは,俺が名前をそうつけただけで,本当に腕の力なのかどうかわからない。なのに我々は,いったん腕力と命名してしまうと,まるでそれが実際に存在するかのように考えて,「腕力だとするとこのモデルはどういう意味があるか」と思考の向きが逆転させてしまう。それは(よくあることだが),よくないことだ。よくよく注意しなければならない。】

というようなことを教える好例となったのである。

これ,探索的因子分析から教えていると,ちょっと気付きにくいことなんじゃないだろうか。あるいは,教える側からは強調するのだけど,教わる側からするとイメージしにくい問題ではないだろうか。そもそも探索的因子分析から教える場合,共通性の推定方法とか,因子数の決め方とか,次々に新しい言葉,気をつける点が出てきてしまって,「因子の解釈と命名」まで注意が維持しにくいような気がするんだよね。これをいい例で教えられたと思ったわけです(間違えてくれた学生ありがとう)。

心理学の研究も気をつけないとそういうことになりがちで,自分で勝手に仮定した構成概念に関係する文言からなる質問紙をして,因子分析の結果,「ほら見ろ因子がでたぞー!」っていうのって,自作自演というか,マッチポンプというか,「てっちゃんの手品」になってるんだよね。

注)「てっちゃんの手品」というのは,愚息(てっちゃん)が「お父さん見てみて!いい?ぬいぐるみをこの箱の中に隠して・・・ジャーン!ぬいぐるみがでてきましたー!」という手品を披露してくれたことに起因する小杉ゼミ専門用語。

 

なので,SEMから教える,SEMだけで教えることを考え始めた。

回帰分析,重回帰分析,パス解析をへて潜在変数を仮定したパス解析,という順に教えていくと,流れが非常にスムーズだし,評価も適合度だけでいいので伝わりやすい。従来の方法でも,どうせ回帰分析,因子分析を教えたあと,構造方程式モデリングに言及するわけだし,しかも最後の一コマだけつかって,「今まではこういう方法があるって教えてきましたけど,これからはこちらに統合されていきます」と言ってもあまり感動してもらえない。今後はもうSEMだけでいいだろうし,いつまでも古典的な方法を従来の順番通り教えていく必要もないかなあ。学部教育なので「因子分析とか回帰分析とかは聞いたことがないです。構造方程式モデリングしかできません」という学生を作るのもちょっとな,と思っていたんだけど,これからはもうそれでいいかもしれない。

 

来年度のデータ処理法,シラバス大幅変更の予感,です。

 



RとMatlabとOctaveと精度

とある論文をもとにして,Rで計算できるようにコードを書いているんだけど,ちゃんと書けたかどうか自信がないw

ので,恥ずかしながら著者の先生に連絡して,検算してもらえないかお願いしてみたところ,快く引き受けてくださったが,何より元のプログラムがFortranで書かれているようで,現在の環境下では動かないとのこと。で,ありがたいことに,現在お使いのMatlabでコードを書いてみましょう,とお願いしてお返事をいただいたのが八月。ありがたいことに,こちらがお渡ししたサンプルコードにくわえてMatlabのコードも送ってくださったのです。

 

ところが,自分の計算結果と検算する間もなく学会シーズンが始まってしまった。半月ほどあいたが,今やっと時間が取れたので検算してみた。

 

案の定,しょうもないミスがあり(行と列を逆にしてた),それを修正するとだいたい合うんだけど,ちょっとやっぱり違うところが出てくる。

それがまあ奇妙なもので,行列は同じ,固有値は同じなんだけど,固有ベクトルが違うというもの。これは正規化するかしないかとか,計算アルゴリズムの問題だよなあ,どうやって検算しようかと頭をひねった。

ということで,まずMatlabのクローンといわれているOctaveをインストール。これをインストールするために,MacPortsをインストールして,失敗して,探しまわったらdmgファイルが落ちていて何じゃそらってなって,という時間を過ごした後,計算してみた。確かにクローンというだけあって,Matlabのコードはほぼコピペで動く。さて固有値・固有ベクトルというところで,同様の問題が出た。固有ベクトルがちょっと合わない。何だこれ,と思ってみたら,固有値分解するプログラムがMatlab(Octave)は二種類あるんですね。

eig関数とeigs関数で,使うアルゴリズムが違うみたい。後者はスパース行列等にも対応しており,計算パッケージとしてはARPACKをベースにしているとのこと。eigs関数の方を使うと,OctaveでもMatlabと近い(それでも完全一致はしない,小数第一桁目までは合う)数値がでる。

Rは基本,eigen関数だけで,LAPACKがベースになっているようだ。そしてRでARPACKを使うにはigraphパッケージを使う必要があり,そこのARPACK関数を使うとある。ところがこのARPACK関数,正方行列を渡したら固有値分解してくれるんじゃなくて,関数の形でベクトルを預けると,一般化固有値問題を解いて指定した数の固有値を出す・・・そんな使い方をするようだ。たしかに社会ネットワークなど,粗なネットワークならそういう使い方もいいだろうけど,こちらはそういう使い方は嬉しくないわけで。

ということで,ここで手をとめることにした。数値計算科学まで首を突っ込んでいられるほど,時間がないんですよw

まあプログラムの概略はわかったし,だいたい自分が間違えてないことが理解できたからいいかな。あとは俺がMatlabを買うか,Rがもっと発展するかですね!

 



階層線形モデル(Hieralical Liner Model)の実例

今度は午後の部,階層線形モデルの話。
ネストされたデータは全部HLMの土俵。反復測定も、個人と集団も。後者の方が得意らしいけど。 ネストされたデータで,HLMを使わないと、サンプルの独立性の検定に違反する。また、平均値が集団の性質を反映していない。後者はカップルデータのようなときに顕著。グループ内の類似性を評価し、それに合わせたモデリングをする。それがHLM,というお話。

ところで,階層回帰分析と階層線形モデリングは,同じ階層という言葉を使っているけど,意味が全然違う。

前者は手続きが順番に行われる,という意味で階層的であり,まず要因A,要因Bを入れ,次に交互作用ABを入れる,という順番でステップを踏む。あくまでも順番であって,あんまり層を積み重ねるという意味じゃないと思う。実際英語としては,Multiple Linear Regression,つまりただ「重回帰分析」と表されるね。重回帰分析は二つ以上の要因がある,という広い意味で使われる。交互作用項を特に【調整変数】と呼んで,センタリング等の適切な処置を経て投入する,というところが特徴。
個人的には「順番にやる回帰分析」「一歩ずつ重回帰」,といった名称にすればよかったんじゃないかと思う。
(実は内心,両者の区別が分からないやつが間違って使い始めたのが定着したのではないかと疑っている)

さて,HLMは本気で階層的。データがレベルを持っている。例えば個人のデータは集団の性質、個人の性質、誤差からなる。場合によっては集団と個人の交互作用も考えることがある。グループ内の類似性が高いことは、集団レベルの情報を多くもっているということ

ここでMLRと混乱させられるもう一つの秘密が。それは,HLMでもステップ1,ステップ2という用語があり得るのです。ほんとはレベル1,レベル2というのが正しい。レベル1は個人のデータ。個人レベルの情報。Within(群内)とも呼ぶ。レベル2は集団のデータ。集団レベルの情報で,Betweenとも。変数がどっちレベルで集まっているか,がしっかり把握できていないと混乱するぜ。HLMはレベル2の変数がレベル1の係数を予測する(回帰のパスが下位レベルの係数にささる)こともでき,それを「交互作用」と呼んじゃうからさらにMLRとの混同がおきやすいよね。要注意。

 

さて,実際にデータを触りながらやった方が分かりやすいかと思うので,Rソースを示し,一歩ずつ解説をしていきます。
サンプルデータは上の記事と同じで,使うデータは清水先生のサイトにあるHLM用サンプルデータをご利用ください。分析資料も同じサイトにあるよ。

RでHLMをするには,lme4というパッケージが必要です。

ソースはこんな感じ。

[crayon-5bc47a310b484938186367/]

さて解説。最初の数行は割愛させてください。センタリングの話から。
HLMでもセンタリングをやります。個人レベルの変数に行うセンタリングは,グループごとの平均値で中心化するセンタリングで「グループセンタリング」といいます。集団レベルの変数に行うセンタリングは,全体の平均値を使うやつで「グランドセンタリング」といいます。

それがここ。

[crayon-5bc47a310b496652546166/]

mean関数は平均値を出す関数。na.rmは欠損値を外せというオプション。もう一つ,グループセンタリングはave関数というのがあって,変数Group毎に平均値を出すという便利な関数。

次の行は,集団平均値をデータとして使う場合があるので,そのための細工。

[crayon-5bc47a310b4ad853579727/]

グループごとに平均を出して,さらにそれをセンタリングしたものを入れています。ave関数がちょっとおかしな形になっているけど,これはデフォルトだとna.rmが入らないから。
ave関数は本来,ave(x,id,FUN)と書いて,データに,群ごとに,ある処理FUNをする,というもの。FUNはデフォルトでmeanなんだけど,このmeanのデフォルトがna.rm=Fなので,それを教えてあげないといけない。つまり,FUNはfunction(x)を使うよ,そしてそれはmean(x,na.rm=T)だよ,という二度手間構造(もっといいやり方があったら教えてエロい人)。

ともかく,これで準備オーケー。
まずは群ごとに,idtを従属変数,talk,perを独立変数とした回帰分析をしてみる(清水先生の分析モデル1)。

[crayon-5bc47a310b4b8056827705/]

関数は,パッケージmlmRevに入っているlmerで,書き方はlmと同じ,チルダの左に従属変数,右に独立変数。独立変数はセンタリングしたtalk_cとper_gを入れる。階層性を表す変数Groupは係数1で関わってきますよ,というのが(1|Group)の意味。

次に,集団平均値を入れたモデル(清水先生資料の分析モデル2)
これはさっきの変数を使うだけだから簡単。

[crayon-5bc47a310b4c1644613613/]

Model.1に比べて,talk_g_mという独立変数が増えただけです。

最後に,ランダム係数モデル(清水先生の分析モデル3)。ここでのランダムは乱数という意味ではなく,無作為でもなく,確率変数(random variables)という意味でのランダムね。

[crayon-5bc47a310b4cb115627549/]

ここでは,Groupごとにtalk_cの係数が変わってくるよ,ということと,交互作用項talk_c*per_gが増えています。

いずれもHADというエクセルマクロで下準備し,SPSSで処理をするという(資料上の)一連の流れがこれだけで再現できるので,どこにどのような数値が表れているかと確認しながらやってみてください。

級内相関の出し方なんかについては,パッケージを使えば出来るんだけど,それまたちょっと調べてから後日ブログにアップする予定です。

ひとまず。

 

追記(2012/02/14)
清水先生のスライドだと,まず級内相関をだして階層性を入れる意義を検証しようね,という話がありました。
この級内相関,パッケージで簡単にでるかなとおもったけど,色々調べても手計算している例が多い。
清水氏のHADみたいに,まず各変数において級内相関を(アルファ係数までも!)出すようなかんすうがあればいいんですけどね。
ひとまず,手計算のやり方を書いておきます。

[crayon-5bc47a310b4d6167901905/]

とまぁ,このように,まずグループレベルの変数しか入ってない回帰モデルを作る

[crayon-5bc47a310b4e1256475265/]

と,その結果として集団レベル変数の分散と残差の分散がでるので,集団レベルの分散/(集団レベルの分散+残差の分散)
を手計算してやる

[crayon-5bc47a310b4eb094792883/]

いいってことになります。

以上。




top