Rで共分散構造分析 {lavaan}パッケージ編 #rstatsj
Rで共分散構造分析を行うにはこれまでは John Fox さんの{sem}パッケージを使うのが普通でした。
私も以前、{sem}パッケージの使い方についてエントリーを書いています。
id:bob3:20091226
id:bob3:20091227
id:bob3:20100530
しかし最近は{OpenMX}や{lavaan}など共分散構造分析のためのパッケージが他にも出てきました。
特に{lavaan}はまだβ版ながら評判が良いようです。
{lavaan}の公式サイトはこちら→ lavaan - latent variable analysis
というわけで、ブルーバックス『原因をさぐる統計学』に載っている事例をいくつが{lavaan}でやってみましょう。
- 作者: 豊田秀樹,前田忠彦,柳井晴夫
- 出版社/メーカー: 講談社
- 発売日: 1992/07/15
- メディア: 新書
- 購入: 6人 クリック: 30回
- この商品を含むブログ (17件) を見る
まずは{lavaan}と関連パッケージのインストール
install.packages(c("lavaan", "psych", "qgraph"))
パッケージの呼び出し。
library(lavaan) library(psych) library(qgraph)
まずは『原因を探る統計学』のp129、p237に載っている「ヘッド・スタート計画」の事例から。
まず、分析する相関係数行列を入力します。
# ヘッドスタート計画の相関係数行列 HS <- matrix(c(1.000, 0.484, 0.224, 0.268, 0.230, 0.265, 0.484, 1.000, 0.342, 0.215, 0.215, 0.295, 0.224, 0.342, 1.000, 0.387, 0.196, 0.234, 0.268, 0.215, 0.387, 1.000, 0.115, 0.162, 0.230, 0.215, 0.196, 0.115, 1.000, 0.635, 0.265, 0.295, 0.234, 0.162, 0.635, 1.000), 6) colnames(HS) <- rownames(HS) <- paste("x", 1:6, sep="")
それぞれの変数の意味は、
x1:母親の学歴
x2:父親の学歴
x3:父親の職業
x4:家庭の収入
x5:ＭＲＴ（知能検査）
x6:ＩＴＰＡ（知能検査）
次にモデルの記述です。
「多重指標モデル」と「外生潜在変数が２つのモデル」の２パターンでやってみます。
# 多重指標モデル HS.mdl.1 <- ' # 測定方程式の記述 f1 =~ x1 + x2 + x3 + x4 f2 =~ x5 + x6 # 構造方程式の記述 f2 ~ f1 ' # 外生潜在変数が２つのモデル # 因子間相関は記述しなくてもOK。 HS.mdl.2 <- ' # 測定方程式 f1 =~ x1 + x2 f2 =~ x3 + x4 f3 =~ x5 + x6 # 構造方程式 f3 ~ f1 + f2 '
{lavaan}パッケージではモデルの記述に以下のような記号の使い方をしています。
=~ 測定方程式
~ 構造方程式（回帰）
~~ 残差の共分散（相関）
それでは分析の実行です。
HS.fit.1 <- sem(HS.mdl.1, sample.cov=HS, sample.nobs=1800) summary(HS.fit.1, fit.measures=TRUE, standardized=TRUE) HS.fit.2 <- sem(HS.mdl.2, sample.cov=HS, sample.nobs=1800) summary(HS.fit.2, fit.measures=TRUE, standardized=TRUE)
出力結果は割愛しますが、CFI、TLI、AIC、BIC、RMSEAといった指標値に続き、推定された係数が出てきます。（なぜかメジャーなGFI、AGFIは出力されない…）
Std.allというのが標準化係数のようです。
また、sem()を実行したときに「相関係数行列からの分析はまだサポートしてないよ！」という警告メッセージが出ます。
このあたりはβ版ゆえ仕方がないでしょう。
それでもテキストと同じ数値が出てきます。
一応、パス図の描画もしてくれます。
ただ、本当に“一応”という感じで、モデルの確認程度にしか使えないと思います。
diagram(HS.fit.1, "多重指標モデル", errors=TRUE, cut=0, lr=FALSE)
qgraph(HS.fit.1, titles=FALSE)
※はてなフォトライフがメンテ中で画像を載せられませんでした…
diagram(HS.fit.2, "潜在変数が２つのモデル", errors=TRUE, cut=0, lr=FALSE)
※はてなフォトライフがメンテ中で画像を載せられませんでした…
qgraph(HS.fit.2, titles=FALSE)
※はてなフォトライフがメンテ中で画像を載せられませんでした…
といった具合で、まだまだ不具合や不満点はあるものの、多母集団同時分析や平均構造モデルにも対応するということで、{lavaan}は今後が楽しみなパッケージです。
時間ができたら他の事例でも試してみます。
あ、そうそう、転職しました。
第3回さくさくテキストマイニング勉強会 #sakutextmining
2011年6月4日(土)に開催された第3回 さくさくテキストマイニング勉強会に参加してきました。
過去２回は青山にあるオラクル社の超豪華な会場をお借りして開催されたのですが、今回からは数理システム社のご厚意により新宿にあるセミナールームをお借りしての開催となりました。
数理システム様に感謝。
内容についてはすでに参加報告されている方々の記事をご参照ください。
id:showyou:20110604 [twitter:@showyou]さん
第3回さくさくテキストマイニング勉強会に参加してきた [twitter:@holidayworking]さん
あの日見た膨大なデータを僕達はまだ知らない [twitter:@marusankak]さん
ここでは発表資料とust録画、写真を記録しておきます。
開場が13:30からだったので、会場近くのトルコ料理店<ボスボラス・ハサン>でランチ。
おいしかった。
会場へ到着〜。
id:AntiBayesianさんによるテキストマイニング入門セッションからスタート。
第三回さくさくテキストマイニング勉強会 入門セッション
View more presentations from AntiBayesian
入門セッションが終わるころには参加者も揃ってほぼ満席。
続いて、[twitter:@toilet_lunch]さんによる「単語重要度入門 〜テキストをダイエットさせよう〜」。
単語重要度入門 〜テキストをダイエットさせよう〜
View more presentations from toilet_lunch
会場の雰囲気はこんな感じでした。
次は[twitter:@Taka_Kuni]さんによる「特徴抽出からクラスタリング」。
３９３９
View more presentations from Takahito Kuniyoshi
次が[twitter:@gepuro]さんの「テキストマイニングの前のコーパス収集」。
Collecting corpus
View more presentations from gepuro
さらに[twitter:@s_wool]さんの「とりあえずTwitterで日本語を集めてみよう」。
そして２回目の登板になるid:AntiBayesianの「KH Coderで3分間テキストクッキング♪」。
今回のust配信も担当していただいた[twitter:@tks]さんの「テキスト/データマイニングと業務」。
最後が「テキストマイニングのマーケティングへの活用について」。
…なのですが、残念ながら大人の事情で非公開です。
ですが、数理システムさん（ sales@msi.co.jp ）宛に「ユーザーコンファレンス2010論文集希望」とメールを送ると近い資料が無料で入手できるかもしれません（笑）。
最後に数理システム社さんからのお知らせ。
紹介された書籍はこちら。
[rakuten:book:14112031:detail]
- 作者: 渕上美喜,末吉正成,高山泰博,今村誠,小木しのぶ,村田真樹,上田太一郎
- 出版社/メーカー: 共立出版
- 発売日: 2008/01/09
- メディア: 単行本
- 購入: 1人 クリック: 96回
- この商品を含むブログ (9件) を見る
[rakuten:neowing-r:10376872:detail]
懇親会も賑やかに盛り上がりました。
Ｒでクロス集計表の残差分析
ぼやーっとツイッターを眺めていたら、[twitter:@R_Linux]先生のツイートが目に留まった。
んん、ということはクロス集計表の残差分析がさくっとできそうだな。
クロス集計表の残差分析はあまりメジャーではないかもしれませんが、χ二乗検定の下位検定として位置付けられています。
端的に言うと「ある一つのセルの観測値が期待値に比べて有意に大きいか小さいか」を検定する手法です。
詳細な解説は下記を参照してください。
群馬大青木先生による解説
雪本さんによる解説
Ｒでの残差分析には id:langstat さんも挑戦しています(id:langstat:20110319)。
青木先生謹製の関数もあります。 カイ二乗分布を用いる独立性の検定
ただ、どちらも調整済み残差を出力するところまでに留まっているので、ここではｐ値の算出と多重比較によるｐ値の調整までを行ってみたいと思います。
サンプルデータは クロス集計表の有意差検定 から引用させていただきます。
# クロス集計表の入力 X <- matrix(c(58, 11, 10, 35, 25, 23), nrow=3, dimnames=list(c("賛成", "中立", "反対"), c("男性", "女性"))) X # 男性 女性 # 賛成 58 35 # 中立 11 25 # 反対 10 23 # 比率（縦％）を確認 round(t(t(X) / colSums(X)) * 100, 1) 男性 女性 賛成 73.4 42.2 中立 13.9 30.1 反対 12.7 27.7 # カイ二乗検定 RES <- chisq.test(X) # 調整残差 RES$stdres # 男性 女性 # 賛成 4.020507 -4.020507 # 中立 -2.478523 2.478523 # 反対 -2.377772 2.377772 # 残差分析のｐ値を算出 p.value.matrix <- pnorm(abs(RES$stdres), lower.tail=FALSE)*2 round(p.value.matrix, 4) # 男性 女性 # 賛成 0.0001 0.0001 # 中立 0.0132 0.0132 # 反対 0.0174 0.0174 # ホルムの方法による有意水準の調整（多重比較） Dim <- dim(X) p.value.matrix.holm <- matrix(p.adjust(p.value.matrix), Dim[1], Dim[2]) dimnames(p.value.matrix.holm) <- dimnames(X) round(p.value.matrix.holm, 4) # 男性 女性 # 賛成 0.0003 0.0003 # 中立 0.0528 0.0528 # 反対 0.0528 0.0528
探索的な分析に結構使い勝手が良いのですよねー。
針穴写真に挑戦中
K10Dで針穴写真に挑戦中。
そのままだといまいちな感じだったのでトイカメラ風にエフェクトをかけています。
レンズ？は韓国製のPinhole Lens 025。
画像の加工にはTOYCAMERA ANALOGCOLORを使いました。
センサーに付いたゴミが悩みの種…
ジャンクレンズ遊び
近所のハードオフで見つけたジャンクレンズ。
Ｋマウントで「SIGMA 1:3.5 f=200mm MULTI-COATED」と表記がある。
レンズにカビのような斑点があってボンヤリ曇っているのがジャンクの理由。
200mmの単焦点はすでにPentaxのK200mm/F4を持っているけど、若干コンパクトなのと半段明るい写りはどんなもんなのかと思って購入。
前玉２枚と後玉２枚は精密ドライバで回したらあっさり外れたので、きれいに中性洗剤で洗ってスッキリ汚れを落とした。
絞りの粘りも何とかしたかったけど、これ以上レンズを外すのは難しそうなので断念。
破損していたフレアカッターを外して、厚紙で作った自作フレアカッターをとりあえずハメ込む。
フードはTakumar200mm/F4のメタルフードを流用。
で、渋谷と池袋で試し撮りしてきた。
開放だと全体にもやーっと白い霧がかかったような写りで、ファインダーを覗いていてもピントが合ってるのか合ってないのか分からない。
順光ならまだましだが、逆光だとLens in a Capの方が良いんじゃないかというぐらい。
f5.6まで絞るとだいぶキリッとする。
そして特徴的なのがぐるぐるボケ。
これはオールドレンズとしてぐるぐるボケを活かした使い方をするのが良さそうですね。
