Rでちゃんとしたクロス集計表を出力する関数

Rで困ることはもうほとんどないのですが,唯一致命的なのが,まともなクロス集計表を出力する関数がないことです。論文で標準的な形式の,周辺度数込みで度数と行パーセントを並べているだけのものでいいんですが,tableやらprop.tableやらaddmarginsやらを駆使してもなかなか思い通りのものが作れない。gmodelsとかのCrossTable関数は悪くはないんですが,返り値がtableじゃない,オプションをいっぱい指定しないと好みのものができない,3重クロス表に非対応,など不満点はけっこう多いです。

パッケージでいい関数が作られてないか,数年間探し続けてたんですが,ラチがあかないので作りました。

使い方は簡単です。

 

> # 擬似データ作成
> 性別 <- factor(rbinom(1:1000, 1, 0.5), labels=c("男性", "女性"))
> 年齢 <- factor(rbinom(1:1000, 2, 0.4), labels=c("若年層", "中年層", "高年層"))
> 身長 <- factor(rbinom(1:1000, 2, 0.6), labels=c("低", "中", "高"))
> 
> # CrossTableオブジェクト作成用関数。使い方もオプションもtable関数とまったく同じ。
> t1 <- crossTable(性別, 年齢, 身長)
> 
> # クロス集計表出力
> summary(t1)
===================================================
                               身長                
                       ---------------------       
                          低     中     高   Total 
---------------------------------------------------
性別 男性 年齢  若年層      41     99     68    208
                         19.7%  47.6%  32.7%   100%
                中年層      29    116     85    230
                         12.6%  50.4%  37.0%   100%
                高年層      13     47     32     92
                         14.1%  51.1%  34.8%   100%
          -----------------------------------------
          Total             83    262    185    530
                         15.7%  49.4%  34.9%   100%
     ----------------------------------------------
     女性 年齢  若年層      22     75     62    159
                         13.8%  47.2%  39.0%   100%
                中年層      33    111     89    233
                         14.2%  47.6%  38.2%   100%
                高年層       8     35     35     78
                         10.3%  44.9%  44.9%   100%
          -----------------------------------------
          Total             63    221    186    470
                         13.4%  47.0%  39.6%   100%
===================================================

Chi-Square Test for Independence

性別 : 男性

Number of cases in table: 530 
Number of factors: 2 
Test for independence of all factors:
    Chisq = 4.494, df = 4, p-value = 0.3432

性別 : 女性

Number of cases in table: 470 
Number of factors: 2 
Test for independence of all factors:
	Chisq = 1.4735, df = 4, p-value = 0.8313

crossTable関数は使い方もオプションもtable関数とまったく同じです。crossTable関数はCrossTableオブジェクトを返しますが,CrossTableクラスはtableクラスを継承しているので,tableとまったく同じように使えます。summary関数に渡すと,クロス集計表と,独立性の検定の結果を出力します。クロス集計表の形式は完全に僕の好みに合わせているので,異論があれば言ってください。なにかオプションを付けるかもしれません。

これで,Rでもカテゴリ変数を心置きなく使えるようになりました。

semパッケージが激しく進化!version3になり超実用的に

semパッケージはオワコン。通はOpenMxかlavaanを使う。そう思っていた時期が私にもありました。

ところが,semパッケージがこの半年の間に2回もメジャーバージョンアップを行っており,凄まじい進化を遂げていました。こういう定番パッケージは大幅な更新はないものと思い込んでいましたが,恥ずかしい見当違いでした。

ざっとヘルプに目を通したところ,以下のようなすばらしい更新点があります。

  • 待ち望まれていた多母集団同時解析に対応!
  • CFA用のモデル記法の追加
  • LISREL風のモデル記法の追加
  • モデルの結合や更新の機能
  • 誤差項を自動的に追加するオプション
  • その他多数

以下,Osaka.R #4の発表資料に対応させる形で解説していきます。

LISREL風記法について

まず,LISREL風記法については,回帰式をそのまま書くような記法が使えるようになりました。

ものすごく簡潔でわかりやすいですね!説明の必要もあまりないくらいですが,左辺に従属変数,右辺にパラメータ名*独立変数の形で回帰式を書いていきます。分散(誤差項)はV(変数)=パラメータ名で表し,共分散はC(変数1, 変数2)=パラメータ名で表します。ちなみにこのときのVとCは大文字でも小文字でもかまいません。

モデルの更新

次に,モデルを更新してみましょう。前述の重回帰分析にパスを足して,パス解析にしてみます。

addはパスを追加する,という命令です。他にもパスを削除するdelete,変数やパラメータ名を置換するreplace,パスを固定するfix,パスを自由にするfreeがあります。ここで注意点として,追加したパラメータ名が元のモデルのパラメータ名と(意図せず)重複しないようにしなければなりません。パラメータ名を重複させた場合,それらのパスには等値制約を課すという意味になります。

モデルの結合

ちなみにupdateでは従来の記法しか使うことができませんが,ただパスを足すだけならば,モデルを結合するcombineModelsを使うことで簡潔な表現が可能になります。

さて,お気づきのとおり,上記のモデルでは教育年数と職業威信の誤差項を書いていません。endog.variancesにTRUEを指定すると,内生変数に誤差項が付いていない場合に自動的に誤差項を付けてくれるようになりました(初期値がTRUEなので実際には上記のようにわざわざ指定する必要はありません)。Mplusでは標準でそうなっているのですが,この機能があるのとないのとではモデル構築の効率が大きく変わってきます。

combineModelsは名前の通りの関数ですが,これもupdateと同じくパラメータ名の重複に注意が必要です。また,この関数はなぜか返り値がmatrixです。ただ,matrixのままでもsem関数に使うことはできます。とはいえ,結合したモデルにupdateを使う予定がある場合などは,ちょっと不便です。そこで,removeRedundantPaths関数を使うと,標準的なsemmodクラスに変換することができます。removeRedundantPaths関数の本来の機能は,名前の通り重複しているパスの削除なのですが,こういう使い方もできるようです。

ここまでLISREL風記法を中心に見てきましたが,注意しなければならないのは,この記法では,構造方程式(回帰の部分)は書けますが,測定方程式(因子分析の部分)は書けません。厳密には書けなくはないのですが,左辺は必ず従属変数かつ1変数でなければならないので,因子分析を表現するためには測定変数の数だけ式を書かなければなりません。

しかし,がっかりする必要はありません。強化されたsemパッケージには,CFA(検証的因子分析)用の記法も用意されています。

CFA用の記法

ものすごく簡潔ですね。reference.indicatorsはTRUEにすると,モデルの識別のために1つめの測定変数のパラメータを1に固定します。初期値はFALSEで,この場合は因子の分散を1に固定します。これだけだと因子間相関や誤差相関はどうすればいいんだ,という話になりますが,covsオプションに"変数1, 変数2"という形で共分散を指定することができます。c("変数1", "変数2")ではないことに注意してください。1つの文字列で1つの共分散を指定するので,複数の共分散を指定するときにc関数などを使ってベクトルで指定する,ということになります。ちなみにcovsオプションはcfa関数限定ではなく,specifyModelやspecifyEquationsでも使うことができます。

応用

ここまでのモデルを結合し,少し大きめのモデルを作ってみます。

多母集団の同時解析

最後に多母集団同時解析の解説です。

まずグループのカテゴリ数だけモデルを用意します。基本は同じで,同じパラメータ名の場合に等値制約が課せられます。いままでは1つのモデル内だけの話でしたが,多母集団同時解析なので,複数のモデル間でも等値制約が課せられることになります。

次に,multigroupModel関数で複数のモデルをひとまとめにします。groupsオプションは単にグループに名前を付けるだけで,サンプルを分割する変数のカテゴリ名に合わせる必要はありません。もっとも合わせない理由はほとんどないと思いますが,カテゴリの順番に気をつけてください。ちゃんと確かめていませんが,分析の出力はlevels関数の出力の順序で行われていると思います。

あとはsem関数を走らせるだけですが,データの渡し方に注意が必要です。共分散行列で渡す場合には,matrixのlistで渡します。生データを使う場合には,dataオプションにdata.framを指定し,formulaオプションに分析で使用する変数を指定します。groupオプションにはグループ分けに使う変数を文字列で指定します。

まとめ

さて,ここまで解説しつつ試した感じでは,いままで私が使ってきたソフトウェアやパッケージの中で,もっともモデルが書きやすかったです。その理由として考えられるのは,「構造方程式,測定方程式,共分散がそれぞれ別の書き方で書けること」なのではないかと思います。複雑なモデルを統一的な数式で表せるのがSEMのすばらしさなのですが,実際に分析するときにはそれが徒となって,モデルが複雑になるにつれどの式が何を表していたのかさっぱりわからなくなってきます。semパッケージの新機能では,因子は因子,回帰式は回帰式ではっきり書き分けることができ,また研究上大した意味を持たないことが多い(そのわりに分量が多い)共分散や誤差項も書き分けで隠蔽したり隔離したりできます。モデルの結合や更新も,SEMの鉄則である,小さなモデルを組み合わせて徐々に大きくしていく,という過程をそのままスクリプトで表現できるので,後から見ても非常に読みやすいように自然となります。

semパッケージは初学者の勉強用にはいいけど研究に使うには物足りないな,とずっと思っていたのですが,これならば万人にお勧めできます。足りないのは,混合効果とカテゴリカルな分析ぐらいでしょうか。そのあたりを必要とするのはかなりマニアックな人たちになってくるので,まあおとなしくMplusを使っておくべきでしょう。

これらの新機能の情報はウェブ等にはほとんどないのですが,ヘルプにはexampleが豊富に付いているので,あまり困ることはありませんでした。help(sem)とhelp(specifyModel)を一通り読んでexampleを走らせれば,だいたい理解できると思います。あとはhelp(package="sem")をざっと見て興味をひかれる関数のヘルプを読めば完璧です。進化したsemパッケージをぜひお試しください!

Rzパッケージ0.4-0リリース

バージョンアップしました。かなり多くの機能を追加したので,足りていない,と思っていた機能が増えているかもしれません。中心的な追加機能を以下に紹介します。

  • プロットビューを大幅更新(実行・クリア・履歴ボタン,変数選択ドロップダウンボックス,ggplot2 0.9.0対応,その他多数)
  • テーマの添付(kde42-oxygen)
  • 変数ビューのレイアウトを変更
  • 尺度の表示および変更方法を変更
  • 変数の削除・複製機能

他にもたくさんありますが(特にプロットビュー周り),自分でもよくわからなくなってしまったので…。パッケージ添付のNEWSにもう少し詳しい更新情報があります。

プロットビューがどんな感じのものなのか,いまいち文章でうまく書けないので,動画を作りました。自分ではそこそこ高機能かつけっこう使いやすいと思っているのですが,いかがでしょうか?

インストール時の注意点として,ggplot2が0.9.0へのアップデートでR (≥ 2.14)を要求するようになりました。それにともない,Rzも ggplot2 (≥ 0.9)かつR (≥ 2.14)が必須となります。そのため,環境によっては,Rzインストールする前にRおよびggplot2のバージョンアップが必要な場合があります。

また,前バージョンから更新した場合には,今バージョンから同梱しているテーマは自動的には適用されません。設定から手動で変更してください。

さらに,更新とは直接関係はありませんがチュートリアルを移転しました。また,Gumroadでの寄付の受付を開始しました。 このパッケージがあなたの研究,仕事,趣味に役立ったのなら,それだけで十分うれしいのですが,寄付していただけるともっとうれしいです。 金額は自由に設定できます。 クレジットカードで簡単に寄付できます。 いただいた寄付は学費に充てさせていただきます。そういえば今月末に結婚式を挙げる予定なのですが,ご祝儀代わりに(ry

ぜひ感想をお聞かせください。

Rzパッケージ0.3-8リリース

Rzパッケージ更新しました。ラベル編集機能,キーボードショートカット,そしてクイックエディタが目玉です。

f:id:phosphor_m:20120207134539p:plain

f:id:phosphor_m:20120207135046p:plain

f:id:phosphor_m:20120207140814p:plain

クイックエディタは複数変数の定型加工やでかいデータを探索的に分析するときには非常に強力です。詳しくはチュートリアル( http://m884.hateblo.jp/entry/2012/02/03/081945 )で。

ちなみに,バージョンアップ頻度が高すぎる,とCRANのメンテナに怒られたので,しばらくバージョンアップはありませんw どっちにしろ,自分が欲しかった機能はほぼ実装できたので,一段落です。

Rzパッケージ最新版 on CRAN Now!

大幅更新を行ったRzパッケージの最新版(0.3-6)が,CRANで公開されました。今のところ,まだ統数研CRANミラー(tokyo)にしかミラーリングされていませんが,今日明日中には筑波と兵庫にも行き渡ると思います(注:0.3-3というのが筑波と兵庫にミラーリングされていますが,これはバグがあって動かないバージョンです)筑波と兵庫にも行き渡りました。

かなり気合いの入った更新でソースコード的にはあまり原型を留めていないので,どこが更新されたのか列挙するのは難しいのですが,チュートリアル(リンクがうまく貼れていません,チュートリアルは1つ前の記事です)を作成したのでこちらを見てもらえればだいたい機能は網羅しています。プロット周りはまだぜんぜん書けていませんが…。

一言(?)で言えばこんなこと↓ができるようになりました。

f:id:phosphor_m:20120203083440p:plain 

あと地味に大きいのが「data.setとして同期」という機能で,これによってコンソールとRzを自由に行き来できるようになりました。