Rを使って統計処理 その1

　C#で統計処理(T検定とか重回帰分析とか）を行う必要が出てきたので統計処理を行うライブラリを探してみたのですが、フリーで使える統計処理ライブラリを見つけられませんでした。
　数学計算用ライブラリを使って重回帰分析を行ったりしている方は何人かいるようですが、統計にそんなに詳しくないのでややこしい数式を自分で計算できる自信がありません。（○○分布の計算とかね）
　そこで統計処理用言語(環境）として有名な「R」を利用することにしました。
　幸いC#からRを操作するためのライブラリ(R.NET)もあるようなので、このR.NETを使って統計処理を行ってみました。(R.NETはNuGetからインストールできます）
　以下のサンプルはR.NETを使って重回帰分析を行うコードです。
　R.NETを使って統計処理を行う為にはRをインストールしておく必要があります。
　Rの外部パッケージを使用したい場合はR側でパッケージをインストールしておく必要もあります。
　今回のサンプルコードではRのcarパッケージを使用したいので、R側でcarパッケージのインストールを行っています。

using RDotNet;
using System;
using System.IO;
using System.Linq;

namespace RNetDemo
{
    class Program
    {
        static void Main(string[] args)
        {
            // R.dllがあるフォルダを環境変数PATHに設定
            var envPath = Environment.GetEnvironmentVariable("PATH");
            envPath += Path.PathSeparator + @"C:\Program Files\R\R-3.2.3\bin\x64";
            Environment.SetEnvironmentVariable("PATH", envPath);

            using (var r = REngine.CreateInstance("R"))
            {
                r.Initialize(); // 初期化は必須

                // ライブラリを読み込む
                r.Evaluate("library(car)");     // VIF計算に必要

                // Rエンジンに数値ベクターを渡す
                var weight = new double[] { 50, 60, 65, 65, 70, 75, 80, 85, 90, 95 };
                var weightVector = r.CreateNumericVector(weight);
                r.SetSymbol("weight", weightVector);

                var height = new double[] { 165, 170, 172, 175, 170, 172, 183, 187, 180, 185 };
                var heightVector = r.CreateNumericVector(height);
                r.SetSymbol("height", heightVector);

                // 数値ベクターを作成するコマンドを実行させても良い
                r.Evaluate("waist <- c(65, 68, 70, 65, 80, 85, 78, 79, 95, 97)");

                // データフレームを作成する
                r.Evaluate("bodyData <- data.frame(weight, height, waist)");

                // 重回帰分析を行い、結果を取得する。
                var regResult = r.Evaluate("reg <- lm(weight ~ height + waist, bodyData)").AsList();
                var summaryResult = r.Evaluate("summary(reg)").AsList();

                // VIFを取得(carライブラリを使用)
                var vifResult = r.Evaluate("vif(reg)").AsList();

                var resultPath = Path.Combine(Environment.GetFolderPath(Environment.SpecialFolder.Desktop), "Result.txt");
                using (var writer = new StreamWriter(resultPath))
                {
                    // 自由度調整済み決定係数
                    var adjR2 = summaryResult["adj.r.squared"].AsNumeric().First();
                    writer.WriteLine("Adj-R2\t{0}", adjR2);
                    writer.WriteLine();

                    // 偏回帰係数
                    var coefficients = summaryResult["coefficients"].AsNumericMatrix();
                    var colNames = coefficients.ColumnNames;
                    writer.WriteLine("Name\t" + string.Join("\t", colNames));
                    foreach (var rowName in coefficients.RowNames)
                    {
                        var valueQuery = colNames
                            .Select(x => coefficients[rowName, x]);
                        writer.WriteLine("{0}\t{1}", rowName, string.Join("\t", valueQuery));
                    }
                    writer.WriteLine();

                    // VIF(carライブラリの関数結果)
                    writer.WriteLine("VIF");
                    writer.WriteLine(string.Join("\t", vifResult.Names));
                    var vifValues = vifResult.AsNumeric().ToArray();
                    writer.WriteLine(string.Join("\t", vifValues));
                }
            }
        }
    }
}

　以下は同じ計算をRで行った結果とこのサンプルコードの結果です。


　R.NETを使う為には環境変数PATHにR.dllのパスを通す必要があります。
　このサンプルではC#側で既存のPATH環境変数にR.dllへのパスを追加しています。

　Rへの操作はRエンジンを通して行います。
　REngine.CreateInstanceメソッドでRエンジンを作成した後、初期化を行います。
　Rエンジンへの基本操作はCreateNumericVectorメソッドなどを使用してデータをR用のベクターに変換し、SetSymbolメソッドでRエンジンにデータを設定します。
　Evaluateメソッドを使うことで指定したRコマンドを実行することができます。
　サンプルのようにEvaluateメソッドを使ってベクターを作成することもできます。
　Rエンジンではデータフレームは作れないようなので、データフレームに使う為のベクターを全て設定したのちEvaluateメソッドを使ってデータフレームを作っています。

　データが用意できたら重回帰分析を行います。
　Evaluateメソッドに重回帰分析を行うコマンドを指定します。その際に戻り値を受け取ることでEvaluateメソッドの計算結果を受け取ることができます。
　Evaluateメソッドの結果はSymbolicExpression型で帰ってきますが、このままだと何もできないのでAsXxxメソッドで適した形に変更する必要があります。
　t値やp値のように数値が返ってきているのであればAsNumericメソッドで数値用ベクターに変換します。
　重回帰分析の計算結果には様々なデータが入っています。
　このような場合はAsListメソッドでGenericVector型に変換することで各データにアクセスできるようになります。
　ここからデータを取り出したい場合はインデクサで取り出したいデータを指定し、AsXxxメソッドで変換を行います。
　例えば、summary結果からadj.r.squaredの値を取り出したい場合は以下の様にします。

var adjR2 = summaryResult["adj.r.squared"].AsNumeric().First();

　偏回帰係数のsummary結果のようにテーブル状になっている場合はAsXxxMatrixメソッドを使います。
　XxxMatrixインスタンスにはColumnNames、RowNamesというプロパティがあるので、この情報を元にデータにアクセスしてください。
　Rの標準機能以外の場合も同様の方法で結果を取得できます。
　このサンプルではcarライブラリのvif関数を使ってVIF値を取得していますが、やり方は全く同じです。

　大体こんな感じでR操作を行えます。
　ちょっとした統計計算を行いたい場合ならR.NETで処理をするのもありだとは思いますが、私は以下の理由からR.NETを使う方法を諦めました。
１）計算速度がかなり遅い。
　全く同じ計算でもRで直接計算したときよりも目に見えて遅くなります。
２）Rエンジンが１つしか作成できない。
　これは同時に１つだけしか作成できないのではなく、アプリケーション内で１つだけしか作れません。
　作成したインスタンスをDisposeして、新たにRエンジンを作った場合、初期化時にフリーズします。
　そのため、最初に作ったRエンジンをアプリケーションが終了するまで使いまわす必要があります。
３）並列処理できない。
　２）の問題があるため、並列処理ができません。
　（R側で並列処理できるのであれば何とかなるかもしれません）

　他にも細かい事はありますが、これが主な理由です。
　今回やりたかったことがブートストラップ法によるリサンプリング->解析なので試行回数が数千回とかになります。
　それを並列処理でやりたかったのですが、R.NETではこれができないので諦めました。
　次回は、R.NETを使わないでRで処理する方法（今回の方法の代替案）を紹介する予定です。