Spearmanの順序相関係数の計算とその検定

あるn組のデータセット $x_i$ ， $y_i$ ， $(i=1,\dots,n)$ について，その値に相関があるかどうかを調べる場合，Pearsonの積率相関係数を利用するのが最も標準的な手続きです．とはいえ， $x_i$ ， $y_i$ が非線型な関数で記述される場合，または外れ点が存在する場合などはSpearmanの順序相関係数を用いた方がより好ましい場合があります．そこで，今日はこのSpearmanの順序相関係数を計算する関数を実装してみます．

Spearmanの順序相関係数はデータ $x_i$ ， $y_i$ をそれぞれ $r_{x,j}$ ， $r_{y,j}$ と昇順，または降順に並べかえたデータについてPearsonの積率相関係数に相当する評価量 $r_s$ を計算します．
$r_s=\frac{\sum_{j=1}^n\left(r_{x,j}-\overline{r_x}\right)\left(r_{y,j}-\overline{r_y}\right)}{\sqrt{\sum_{j=1}^n\left(r_{x,j}-\overline{r_x}\right)^2}\sqrt{\sum_{j=1}^n\left(r_{y,j}-\overline{r_y}\right)^2}}$
ただし， $\overline{r_x}$ ， $\overline{r_y}$ はそれぞれ $r_{x,j}$ ， $r_{y,j}$ の平均とします．

相関係数の有意性の検定はデータに相関がないという帰無仮説 $H_0$ の元で行うことが考えられます．n > 30の場合，この仮説のもとでは以下の評価量が自由度n-2のt分布に従います．
$t = r_s\sqrt{\frac{n-2}{1-r_s^2}}$
したがって，自由度n-2のt分布の両側確率を計算し，この値が十分に小さければ $H_0$ を棄却でき，データセットに相関がある可能性が示されます．なお，n <= 30の場合は有意水準ごとに用意された表を用いて検定を行います．

サンプルコードは以下です．

#include <iostream>
#include <algorithm>
#include <iterator>
#include <vector>
#include <numeric>
#include <functional>

// DCDFLIB
#include "dcdflib"

// Spearman の順序相関係数の計算
// STLの使い方に無駄がある気がする
int Spearman(int n, const double *x, const double *y,
	      double &rs, double &pvalue){

  using namespace std;
  int info;
  vector< pair<double,int> > vx, vy;

  if (!(n > 30)){
    info = -1;
    return info;
  }

  // pairを作ってコンテナに格納
  for(int i = 0; i < n; i++){
    vx.push_back(make_pair(x[i], i));
    vy.push_back(make_pair(y[i], i));
  }

  // ソートする
  sort(vx.begin(), vx.end());
  sort(vy.begin(), vy.end());

  // 順位を格納
  vector<double> rx(n);
  vector<double> ry(n);
  for(int i = 0; i < n; i++){
    rx[vx[i].second] = (double)i;
    ry[vy[i].second] = (double)i;
  }

  // 平均を計算
  double mx = accumulate(rx.begin(), rx.end(), 0.0) / rx.size();
  double my = accumulate(ry.begin(), ry.end(), 0.0) / ry.size();

  // 平均減算する
  vector<double> rx2(rx.size());
  vector<double> ry2(ry.size());
  transform(rx.begin(), rx.end(), rx2.begin(), bind2nd(minus<double>(), mx));
  transform(ry.begin(), ry.end(), ry2.begin(), bind2nd(minus<double>(), my));

  // 順序相関係数を計算
  double cxy = inner_product(rx2.begin(), rx2.end(), ry2.begin(), 0.0);
  double cxx = inner_product(rx2.begin(), rx2.end(), rx2.begin(), 0.0);
  double cyy = inner_product(ry2.begin(), ry2.end(), ry2.begin(), 0.0);
  rs = cxy/(sqrt(cxx)*sqrt(cyy));

  // t検定 (両側検定)
  int which, status;
  double p, q, t, df, bound;

  df = n - 2;
  which = 1;
  t = (rs)*sqrt((n-2)/(1-(rs)*(rs)));
  cdft(&which, &p, &q, &t , &df, &status, &bound);
  // void cdft ( int *which, double *p, double *q, double *t, double *df,
  //               int *status, double *bound );
  // 目的         : (Studentの)t分布の累積分布関数の値を計算する
  // 引数         : *which  処理の内容を設定する
  //                        1: pとqの値をt, dfから計算する
  //                        2: tの値をp, q, dfから計算する
  //                        3: dfの値をp, q, tから計算する
  //              : *p      負の無限大からxまでの積分値
  //              : *q      1 - p
  //              : *t      積分値を計算する場合の上限
  //              : *df     t分布の自由度
  //              : *status 関数が終了した時点でのステータス
  //                        0 : 正常終了の場合は0が格納される
  //                        -i: i番目のパラメターがレンジ外になると負数iが格納される
  //                        1 : 解が探索領域の最小値を下回った
  //                        2 : 解が探索領域の最大値を上回った
  //                        3 : p + q が1とならなかった
  //              : *bound  正常終了時 (status = 0)は不定値
  //                        異常終了 (status < 0) ならレンジ外になったi番目のパラメターの限界値
  //                        異常終了 (status = 1) なら探索領域の下界
  //                        異常終了 (status = 2) なら探索領域の上界
  if (!status) {
    // p-valueを計算する
    if(t > 0.0){ // q を使う
      pvalue = 2*q;
    }
    else { // p を使う
      pvalue = p*2;
    }
    info = 0;
  }
  else{
    info = 1;
  }

  return info;
}

#define NS (50)
#define SL (0.05)

int main(void)
{
  using namespace std;
  double x[NS], y[NS], r, p;

  // 乱数源を初期化
  srand((unsigned int)time(NULL));

  // 一様分布の乱数を作る
  for (int i = 0; i < NS; i++){
    x[i] = (rand()/((double)RAND_MAX));
    y[i] = (rand()/((double)RAND_MAX));
  }

  // 数列の相関をSpearmanの順序相関係数で評価し，
  // その有意性をt検定で調べる
  // 帰無仮説H0：データは無相関
  if (!Spearman(NS, x, y, r, p)){
    cerr << "Spearman's rank correlation coefficient is " << r << "." << endl;
    cerr << "hypothesis is ";
    if(p < SL) {
      cerr << "rejected (" << p << " < " << SL << ")" << endl;
    }
    else{
      cerr << "not rejected (" << p << " >= " << SL << ")" << endl;
    }
  }

  return 0;
}

t分布の累積分布関数の計算はDCDLIBのcdftを用いました．STLの使い方に無駄がある気がしますが，これは著者がまだSTLの利用法に習熟していないためです．実行結果は以下です．

Spearman's rank correlation coefficient is 0.124418.
hypothesis is not rejected (0.389312 >= 0.05)

一様乱数は(極端にいい加減に作らない限り)無相関ですので，仮説を棄却しないという結果は概ね妥当でしょう．時折棄却されてしまうこともあるようですが…