Rでパッケージを作れることは分かったので、自分(の授業、のゼミ生)に必要な関数・計算モデルを書き出してみた。
ふむふむ、どうしてもポリコリック相関係数は外せない(笑)
ということで、まずは関数を作るところからはじめにゃなー、とライブラリPolycorをみてみたら、つい先日バージョンアップがなされていたようだ(10/28/2008)。
みると、hetcorという関数がパッケージに追加されている。
これはすごいよ。複数の変数を与えると相関行列を作ってくれるんだけど、二変数が連続値ならピアソンの、一方が連続値で他方が順序レベルであればポリシリアルの、両方が順序レベルであればポリコリックの相関行列を出してくれるのだ。
(訂正)次の一文は、マニュアルの誤読による私の勘違いでした。撤回させていただきます。
行列までいっぺんに出してくれる機能が欲しかったのに、「カテゴリ(bin)が5以上になれば連続値と判断して・・・」というおまけ付き。これがあったら完璧じゃよー!
出力結果を乗せておきます。分かる人にはこの良さが分かるはず。#印は私のコメントです。
<span class="synComment">#二段階推定法をする場合</span> Two-Step Estimates <span class="synComment">#行列の種類と数値(!)</span> Correlations/Type of Correlation<span class="synSpecial">:</span> x1 x2 y1 y2 x1 <span class="synConstant">1</span> Pearson Polyserial Polyserial x2 <span class="synConstant">0.5933</span> <span class="synConstant">1</span> Polyserial Polyserial y1 <span class="synConstant">0.5953</span> <span class="synConstant">0.7409</span> <span class="synConstant">1</span> Polychoric y2 <span class="synConstant">0.6241</span> <span class="synConstant">0.6317</span> <span class="synConstant">0.5711</span> <span class="synConstant">1</span> <span class="synComment">#標準誤差</span> Standard Errors<span class="synSpecial">:</span> x1 x2 y1 x1 x2 <span class="synConstant">0.02051</span> y1 <span class="synConstant">0.0309</span> <span class="synConstant">0.02295</span> y2 <span class="synConstant">0.02025</span> <span class="synConstant">0.01994</span> <span class="synConstant">0.03738</span> <span class="synComment">#サンプル数。例では1000サンプルの乱数を発生させて求めています。</span> n = <span class="synConstant">1000</span> <span class="synComment">#二変数<a class="keyword" href="http://d.hatena.ne.jp/keyword/%C0%B5%B5%AC%CA%AC%C9%DB">正規分布</a>の有意性検定</span> P-values <span class="synStatement">for</span> Tests of Bivariate Normality<span class="synSpecial">:</span> x1 x2 y1 x1 x2 <span class="synConstant">0.5029</span> y1 <span class="synConstant">0.3407</span> <span class="synConstant">0.878</span> y2 <span class="synConstant">0.1022</span> <span class="synConstant">0.5255</span> <span class="synConstant">0.05615</span> <span class="synComment">#<a class="keyword" href="http://d.hatena.ne.jp/keyword/%BA%C7%CC%E0%BF%E4%C4%EA">最尤推定</a>の場合</span> Maximum-Likelihood Estimates Correlations/Type of Correlation<span class="synSpecial">:</span> x1 x2 y1 y2 x1 <span class="synConstant">1</span> Pearson Polyserial Polyserial x2 <span class="synConstant">0.5933</span> <span class="synConstant">1</span> Polyserial Polyserial y1 <span class="synConstant">0.5949</span> <span class="synConstant">0.7406</span> <span class="synConstant">1</span> Polychoric y2 <span class="synConstant">0.6243</span> <span class="synConstant">0.6321</span> <span class="synConstant">0.5715</span> <span class="synConstant">1</span> Standard Errors<span class="synSpecial">:</span> x1 x2 y1 x1 x2 <span class="synConstant">0.02051</span> y1 <span class="synConstant">0.03105</span> <span class="synConstant">0.02301</span> y2 <span class="synConstant">0.02079</span> <span class="synConstant">0.02044</span> <span class="synConstant">0.03842</span> n = <span class="synConstant">1000</span> P-values <span class="synStatement">for</span> Tests of Bivariate Normality<span class="synSpecial">:</span> x1 x2 y1 x1 x2 <span class="synConstant">0.5029</span> y1 <span class="synConstant">0.347</span> <span class="synConstant">0.8607</span> y2 <span class="synConstant">0.1046</span> <span class="synConstant">0.5139</span> <span class="synConstant">0.05675</span>
さぁこれを使えるようにすれば・・・うーん、イエス!!(マトリックスのタンクの感じで)
しかし、本業としての論文書きや処理するデータがたまっていることも事実。Rで遊ぶ時間は一日一時間まで、とかってきめないとにゃー・・・。