好きなRパッケージが`{psych}`と`{lme4}`という話

この記事はR言語 Advent Calendar15日目の記事です。

私は卒論、修論、炎上案件は全部Rで切り抜けてきてい程度にはRを触れていて、かれこれ10年はRと暮らしています。
そうしていると当然、様々なパッケージに出会うものでして、Rユーザは好きなパッケージを述べることがアイサツ代わりになるので*1、今年のアドカレはそういう話をしようかなと思います。

昨今はtidyverseã‚„tidymodelsといったパッケージ「群」が注目されています。
これらはRの良くも悪くもカオスなデータ分析のフォーマットに一貫した思想を持ってメンテナンスされています。 逆に言うと、これら以外のパッケージは、各領域の問題・関心に応じて異なります。 何よりR言語は「統計解析に特化した言語」であるはずなので、それはそれは多様な解析手法が実装されています。ただそこそこ残念なことに「統計解析」のパッケージに関してはそこまで注目されないのが、深層学習や大規模言語モデルが注目されるアルゴリズムになっている現代のデータ分析市場における、少しさみしいところです *2。

仕事が修羅場だし、忙しくてなに書こうかなんともわからないし、どうせn番煎じで、漸近的には内容のない記事になるのですが、せっかくなので世話になったり、シンプルに好きな統計解析パッケージを2つ紹介しようかなと思いました。 数理的な話は統計・機械学習の数理Advent Calendarで書いてほしいし、何かを書きます。

github.com

因子分析パッケージ{psych}

私も長いこと調査データの分析に染まっている人間なので、まずは調査データでよく用いるパッケージの話をします。 あらかじめいいますがこれは個人の見解であって、所属元を代表した意見ではありません。 また、これが正解だとか言うこともありません。好きなパッケージを言わせろ。{psych}。

{psych}は心理測定理論、実験心理学に用いる分析全般を幅広くカバーしたパッケージです。大学でも仕事でも、調査データを使った分析で世話にならないことがないくらい、よく使っています(いました)。

調査データの前提

調査観察データは、国勢調査を始めとした、総務省統計局が行う各種世帯調査はもちろん、 産業界においてもブランドのイメージや世代の価値観の変化、ユーザ観点からのプロダクト評価など、 データの量の膨大化・質の多様化が進んだ現代においても測定の難しい「人の意識や状態」 を測定するために重要な位置を占めています。「マーケティングリサーチ」という業界・業種があるように、 これを測定するための技術・ノウハウが産業的価値をもっていることは一定事実だと思います。

ここでは「消費者がものを買うときに考えていること」を知ることを例に取ります

調査データは、実施する側が知りたいこと、測りたいことを測ることができるという強みがある反面、 たとえば「あなたが商品を買うときに考えることはなんですか？」と聞いても、明確な答えを回答者から得ることは 難しいと思います。そのため「商品を買うときは価格を重視している」とか「有名なメーカーが出しているものを買う」 「環境に優しい商品を買う」など、いくつか選択肢を提示しながら、回答者が商品を買うときにどんなことを意識していそうかを 可能な限り引き出そうと試みます。

そんな中で起きる暗黙のやり取りに「回答者に可能な限りたくさんのことを聞きたい」という実施側のニーズと、 「長いアンケートはだるい」という回答者側の事情(回答負荷)で、これをなんとかうまく折衝すること、 つまり「可能な限り短いアンケートで、可能な限り多くの情報を引き出す」ことがアンケートデータでは求められます。

次元集約手法としての因子分析・主成分分析

こうした意識を測定するときによく用いられるアプローチは「実施する側で考えうる限りの 『ものを買うときに考えていること』の仮説を出し切って、何らかの形で集約する」ことです。 この「仮説を出し切る」というのはまあまあの難易度が高いのですが、今回は統計解析パッケージを紹介する記事なので、 この辺知りたい人はbob3さんとかに聞いて下さい調査会社とかに依頼してください。

さて、どうにか仮説を出し切って、回答してもらったことを仮定して話を進めると、 得られた回答をどのように集約するかを考える必要があります。 もちろん、事前に「価格系の質問」「環境系の質問」のようにきっちり設計して、それをまとめても良いのですが、 その時も「どの設問がどの程度の重みでまとめられうるのか」を解決することは難しいかもしれません。

こうした問題に対して、よく知られた解析手法は「因子分析」や「主成分分析」と言われたものが知られています。 因子分析は特に心理学の領域で見られますし、主成分分析も社会科学分野で用いられてきましたが、 今は機械学習での次元縮約アルゴリズムとしての位置づけが大きいかもしれません。

Rで因子分析をする場合、{psych}パッケージでの実装が知られていますし、多分これが一番シンプルに実行できると思います。 何なら主成分分析はデフォルトのパッケージ{stats}にprcomp()として実装されています。

今回は統計解析パッケージを紹介する記事なので、{psych}パッケージで因子分析をやってみたいと思います。

因子分析も主成分分析も共通して「因子(主成分)の数をどの程度にするか？」は人間が何らかの基準で決めることになります。因子分析の場合、良く用いられる手法はスクリープロットですが、{psych}であればこれも実施できます。 スクリープロットは第一因子の特異値を縦軸、因子数を横軸にとった折れ線グラフで、{psych}ではこれに加えて多変量正規分布に従う分散共分散行列の特異値分解をしたときのシミュレーションデータと比較する事もできます。

library(tidylog)
library(psych)

# psychにあるサンプルデータを使ってみる
data(bfi)

bfi_to_pca <- bfi |> 
  tidylog::select(-gender, -education, -age)

# スクリープロット
psych::fa.parallel(
  bfi_to_pca
)

R> Parallel analysis suggests that the number of factors =  6  and the number of components =  NA 

スクリープロットはこちら

スクリープロット

これを見て、適した因子数はシミュレーションで得られたデータ(赤線)を下回るあたりが適している6～7あたりがちょうどよいと判断します。なんならRも、因子数6を推奨しています*3。

因子数を6と設定したときに因子分析を実施してみましょう。

fa_result <- psych::fa(
  bfi_to_pca, 
  nfactors = 6,
  rotate = "promax",
  fm = "ml")

summary(fa_result)
R> Factor analysis with Call: psych::fa(r = bfi_to_pca, nfactors = 6, rotate = "promax", fm = "ml")
R>
R> Test of the hypothesis that 6 factors are sufficient.
R> The degrees of freedom for the model is 165  and the objective function was  0.36 
R> The number of observations was  2800  with Chi Square =  1013.79  with prob <  4.6e-122 
R> 
R> The root mean square of the residuals (RMSA) is  0.02 
R> The df corrected root mean square of the residuals is  0.03 
R> 
R> Tucker Lewis Index of factoring reliability =  0.922
R> RMSEA index =  0.043  and the 10 % confidence intervals are  0.04 0.045
R> BIC =  -295.88
R> With factor correlations of 
R> ML1   ML2   ML3   ML5   ML4  ML6
R> ML1  1.00 -0.36  0.41  0.33 -0.11 0.09
R> ML2 -0.36  1.00 -0.21 -0.18  0.11 0.26
R> ML3  0.41 -0.21  1.00  0.31 -0.15 0.20
R> ML5  0.33 -0.18  0.31  1.00  0.03 0.28
R> ML4 -0.11  0.11 -0.15  0.03  1.00 0.09
R> ML6  0.09  0.26  0.20  0.28  0.09 1.00

fa_result$weights  # 因子負荷量
R> 各入力変数が、各因子にかかる因子負荷量がでてくる

fa_result$scores   # 因子特点行列
R> サンプル別に、各因子に対する特典が得られる

こんな感じで、基本的な出力はある程度しっかり出してくれます。玉に瑕というものがあるならば、
因子のラベル順序が、必ずしも昇順ではない、というところでしょうか。バグなのか何なのか……。

ちなみに主成分分析と因子分析はどちらも特異値分解を元とした統計解析手法ですが、今回は統計解析パッケージを紹介する記事なので、これらの話は統計・機械学習の数理Advent Calendarに記載できたらいいですね。

回帰分析

機械学習と統計学の共通部分にある、統計解析っぽい問題に「回帰問題」があると思っています。 私は回帰分析が大好きで、回帰分析を食べて生きていると言っても過言ではないでしょう。

Rではデフォルトのパッケージ{stats}でいわゆる普通の回帰分析をlm()、ロジスティック回帰分析など いわゆる「一般化線形モデル」をglm()で実行できます。 何も追加パッケージを入れなくても回帰分析ができるプログラミング言語はそんなにないと思うので、これはすごいと思います。

私が好きなのは更に複雑なモデルで、一般化線形混合モデルと呼ばれるものがあります。 すでに過去に同様の記事を書いているんですが、何度書いても良いので書きます。

一般化線形混合モデルはマルチレベルモデルとかランダム効果モデルとか いろいろな領域で別の名前がついていて面倒ですが、この記事では 「特定のサブグループによってパラメータの効果が変わりうるときに採用するモデル」としてなんとか切り抜けます。

過去に書いた記事でも{lme4}を紹介して来ました。何なら修論の一部にも使ってきています。このパッケージはglmライクに一般化線形混合モデルを記述することができますが、固定効果の回帰係数の有意性検定を出しません。これは色々の議論があります。今の私の理解では「帰無仮説のときに従う分布」の仮定は、本当に実態として満たされるのか？満たされない例が作れていそう。というところから、検定に用いる量は出せるようにしながらも、具体的な検定はあえて実装しない、というのが{lme4}の開発側のスタンスであると考えています。
今回は統計解析パッケージを紹介する記事なので、具体的に回帰係数の統計的検定がどうなっているのかは端折ります。 いずれにしても、{lme4}は混合モデルを大きなモデルに対して高速な計算の元実現できることで、良くも悪くも尖っているところが好きです。なんか、かっこいい。 もしも{lme4}ベースで混合モデルの回帰係数の検定を実施したい場合、{lmerTest}というパッケージがあります*4。 実装については以下の過去記事を見たほうがはやいです。最近は過去の自分を越えられる事ができない。

socinuit.hatenablog.com

そこで出した結果。データが特定の群別に測定されたデータ(入れ子データ、ネステッドデータ)を、1つのモデルでうまいこと表現できるのが、
割とかっこいいな、と思いませんか？私はかっこいいと思います。

github.com

終わりに

私が好きなRパッケージは{psych}と{lme4}です。皆さんの好きなパッケージはなんですか？

参考文献

• 足立 浩平, 山本 倫生『主成分分析と因子分析　特異値分解を出発点として』

• Ben Bolker and others, "GLMM FAQ"

*1:そんなこともない

*2:これらの発展は人類に大きく寄与するものである確信はありますし、測定できるデータの種類と量が増えたことで、 数理統計学が発展した時代の制約の一部が取り払われたことは間違いはないと思います。 これらを認めたうえで、私はこの趨勢が、個人的に、気に食わないというだけです。

*3:実務上はより広めに、5～9くらいの因子数で一気に回して、適切な解釈が可能な要素を出すようなチェリーピッキングもまま行います。

*4:ただし、実施可能なのはlmerのみ。一般化線形混合モデルglmerにおける検定はこれでもやってくれないし、ランダム効果に対する検定は魔境である。