Fortran - スピアマン順位相関係数の計算！

[ プログラミング, 数学 ] [ Fortran ]

Fortran 95 でスピアマンの順位相関係数(Spearman’s Rank Correlation Coefficient)の計算をしてみました。

0. 前提条件

  • Debian GNU/Linux 10.3 (64bit) での作業を想定。
  • GCC 9.2.0 (GFortran 9.2.0) でのコンパイルを想定。

1. スピアマンの順位相関係数について

各変量を順位に変換してピアソンの積率相関係数（いわゆる相関係数）を求めたものを スピアマンの順位相関係数(Spearman’s Rank Correlation Coefficient) と呼ぶ。

実際にはまず、 対の変数 のそれぞれに順位をつける。但し、同順位（タイ）がある場合は中央（平均）順位(mid-rank) で順位をつける。
（e.g. 2位が3個ある場合、 。3位が2個ある場合、 ）

そして、次の式によりスピアマンの順位相関係数 （または ）を求める。

(1) 同順位（タイ）が存在しない場合、

(2) 同順位（タイ）が存在する場合、

（注） 同順位が存在しない場合は (2) の が となり、結局 (1) になる。よって、同順位（タイ）が存在しない場合と存在する場合で場合分けをせず、全て (2) で計算しても、結果は同じになる。

また、スピアマンの順位相関係数は、値を順位に置き換えたもの（同順位（タイ）は中央順位法）の相関係数（ピアソンの積率相関係数）であるので、当然、以下の計算式でも計算できる。

さらに、同順位（タイ）が存在する場合の計算式を以下のように説明している文献（特に海外の文献）も多い。（但し、前述の計算式で計算した結果と若干の差異がある）

2. Fortran ソースコードの作成

  • 計算の本質は comp モジュール。
  • 使用する2つの配列 X, Y を適宜変更する。（calc_rcc_spearman モジュール内）
calc_rcc_spearman.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212

!****************************************************
! スピアマンの順位相関係数の計算
!
! DATE          AUTHOR          VERSION
! 2019.12.12    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2019 mk-mode.com All Rights Reserved.
!****************************************************
!
module cst
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module cst

module comp
  use cst
  implicit none
  private
  public :: rcc_spearman

contains
  ! スピアマンの順位相関係数の計算
  !
  ! :param(in) real(8)  xs
  ! :param(in) real(8)  ys
  ! :return    real(8) rcc
  function rcc_spearman(xs, ys) result(rcc)
    use cst
    implicit none
    real(DP), intent(in) :: xs(:), ys(:)
    real(DP)    :: rcc
    integer(SP) :: i, n, n3
    real(DP)    :: s_x, s_y, t_x, t_y, s
    real(DP), allocatable :: rank_x(:), rank_y(:)
    real(DP), allocatable :: tai_x(:, :), tai_y(:, :)

    rcc = -1.0_DP
    n = size(xs)
    if (size(ys) /= n) return
    allocate(rank_x(n))
    allocate(rank_y(n))
    allocate(tai_x(n, 2))
    allocate(tai_y(n, 2))
    ! ランク付け
    call rank(xs, rank_x)
    call rank(ys, rank_y)
    ! 同順位の数（インデックスが順位）
    call tai(rank_x, tai_x)
    call tai(rank_y, tai_y)
    ! Tx, Ty の sum 部分
    call sum_t(tai_x, s_x)
    call sum_t(tai_y, s_y)
    ! 相関係数
    n3 = n * n * n - n
    t_x = (n3 - s_x) / 12.0_DP
    t_y = (n3 - s_y) / 12.0_DP
    s = 0.0_DP
    do i = 1, n
      s = s + (xs(i) - ys(i)) *  (xs(i) - ys(i))
    end do
    rcc = (t_x + t_y - s) / (2.0_DP * sqrt(t_x * t_y))

    deallocate(rank_x)
    deallocate(rank_y)
    deallocate(tai_x)
    deallocate(tai_y)
  end function rcc_spearman

  ! ランク付け
  !
  ! :param(in)  real(8) src(:)
  ! :param(out) real(8) res(:)
  subroutine rank(src, res)
    implicit none
    real(DP), intent(in)  :: src(:)
    real(DP), intent(out) :: res(:)
    integer(SP) :: c, i, j, s
    real(DP), allocatable :: wk(:)

    s = size(src)
    allocate(wk(s))
    wk = 0.0_DP
    do i = 1, s
      c = 0
      do j = 1, s
        if (src(i) < src(j)) c = c + 1
      end do
      wk(i) = real(c + 1, DP)
    end do
    ! 同順位を中央（平均）順位(mid-rank)に
    do i = 1, s
      c = 0
      do j = 1, s
        if (wk(i) == wk(j)) c = c + 1
      end do
      res(i) = sum((/(j, j=int(wk(i)),int(wk(i))+c-1)/)) / real(c, DP)
    end do
    deallocate(wk)
  end subroutine rank

  ! 同順位の数
  !
  ! :param(in)  real(8)    src(:)
  ! :param(out) real(8) res(:, 2)
  subroutine tai(src, res)
    implicit none
    real(DP), intent(in)  :: src(:)
    real(DP), intent(out) :: res(:, :)
    integer(SP) :: c, i, j, s
    real(DP), allocatable :: wk(:)

    s = size(src)
    allocate(wk(s))
    wk = src
    res = 0.0_DP
    do i = 1, s
      c = 0
      do j = i, s
        if (wk(j) == 0.0_DP) cycle
        if (wk(i) == wk(j)) then
          c = c + 1
          if (i /= j) wk(j) = 0.0_DP
        end if
      end do
      res(i, :) = (/wk(i), real(c, DP)/)
    end do
    deallocate(wk)
  end subroutine tai

  ! Tx, Ty の sum 部分
  !
  ! :param(in)  real(8) src(:, 2)
  ! :param(out) real(8)       s_t
  subroutine sum_t(src, s_t)
    implicit none
    real(DP), intent(in)  :: src(:, :)
    real(DP), intent(out) :: s_t
    integer(SP) :: i, s
    integer(DP) :: t

    s = size(src) / 2
    s_t = 0.0_DP
    do i = 1, s
      t = int(src(i, 2), DP)
      if (t < 2) cycle
      s_t = s_t + t * t * t - t
    end do
  endsubroutine sum_t
end module comp

program calc_rcc_spearman
  use cst
  use comp
  implicit none
  real(DP)    :: rcc
  ! === タイ（同順位）が存在しない例
  !real(DP), parameter :: X(10) = (/ &
  !  & 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP,  5.0_DP, &
  !  & 6.0_DP, 7.0_DP, 8.0_DP, 9.0_DP, 10.0_DP  &
  !& /)
  !real(DP), parameter :: Y(10) = (/ &
  !  & 1.0_DP, 3.0_DP, 5.0_DP, 7.0_DP,  9.0_DP, &
  !  & 2.0_DP, 4.0_DP, 6.0_DP, 8.0_DP, 10.0_DP  &
  !& /)

  ! === タイ（同順位）が存在する例
  real(DP), parameter :: X(10) = (/ &
    & 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP,  5.0_DP, &
    & 5.0_DP, 7.0_DP, 8.0_DP, 9.0_DP, 10.0_DP  &
  & /)
  real(DP), parameter :: Y(10) = (/ &
    & 1.0_DP, 3.0_DP, 5.0_DP, 6.0_DP,  9.0_DP, &
    & 2.0_DP, 4.0_DP, 6.0_DP, 8.0_DP, 10.0_DP  &
  & /)

  ! === サイズが異なる例 ( => 実行時にエラー )
  !real(DP), parameter :: X(10) = (/ &
  !  & 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP,  5.0_DP, &
  !  & 6.0_DP, 7.0_DP, 8.0_DP, 9.0_DP, 10.0_DP  &
  !& /)
  !real(DP), parameter :: Y(9) = (/ &
  !  & 1.0_DP, 3.0_DP, 5.0_DP, 7.0_DP,  9.0_DP, &
  !  & 2.0_DP, 4.0_DP, 6.0_DP, 8.0_DP           &
  !& /)

  ! === X のサイズがゼロの例（ => 当然、コンパイルエラー ）
  !real(DP), parameter :: X(0) = (//)
  !real(DP), parameter :: Y(9) = (/ &
  !  & 1.0_DP, 3.0_DP, 5.0_DP, 7.0_DP,  9.0_DP, &
  !  & 2.0_DP, 4.0_DP, 6.0_DP, 8.0_DP, 10.0_DP  &
  !& /)

  ! === 数値以外のものが存在する例（ => 当然、コンパイルエラー ）
  !real(DP), parameter :: X(10) = (/ &
  !  & 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP,  5.0_DP, &
  !  & 6.0_DP, 7.0_DP, 8.0_DP, 9.0_DP, 10.0_DP  &
  !& /)
  !real(DP), parameter :: Y(10) = (/ &
  !  & 1.0_DP, 3.0_DP, 5.0_DP, 7.0_DP,  9.0_DP, &
  !  & "ABC", 4.0_DP, 6.0_DP, 8.0_DP, 10.0_DP  &
  !& /)

  rcc = rcc_spearman(X, Y)
  print '("X = ", 10F5.1)', X
  print '("Y = ", 10F5.1)', Y
  if (rcc == -1.0_DP) then
    print *, "[ERROR]"
  else
    print '("Spearman''s RCC = ", F11.8)', rcc
  end if
end program calc_rcc_spearman

3. コンパイル＆リンク

1

$ gfortran -Wall -O2 -o calc_rcc_spearman calc_rcc_spearman.f95

4. プログラムの実行

1
2
3
4

$ ./calc_rcc_spearman
X =   1.0  2.0  3.0  4.0  5.0  5.0  7.0  8.0  9.0 10.0
Y =   1.0  3.0  5.0  6.0  9.0  2.0  4.0  6.0  8.0 10.0
Spearman's RCC =  0.70731707

以上。