Rによる級内相関係数

はじめに

本記事ではRでの級内相関係数の算出方法を紹介する。ICC(1,1)、ICC(1,k)、ICC(2,1)、ICC(2,k)、ICC(3,1)、ICC(3,k)のプログラムについては" Rによる算出"で記載した。算出するだけであればRを利用すれば容易である。しばしば目にする級内相関係数の算出方法は平均平方による式で書かれており、個人的にイメージが湧きづらく平均平方を用いることで級内相関係数が算出できる部分については数式を整理した。級内相関係数の大きさの解釈は応用領域やその状況に強く依存するため本記事内では触れない。

【参考文献】

• 実験計画法 (新統計学シリーズ 2)
• Shrout PE, Fleiss JL. Intraclass correlations: uses in assessing rater reliability. Psychol Bull 1979 Mar;86(2):420–428. PMID:18839484
• Koo TK, Li MY. A Guideline of Selecting and Reporting Intraclass Correlation Coefficients for Reliability Research. J Chiropr Med 2016 Jun;15(2):155–163. PMID:27330520
• Liljequist D, Elfving B, Skavberg Roaldsen K. Intraclass correlation - A discussion and demonstration of basic features. PLoS One Public Library of Science (PLoS); 2019 Jul 22;14(7):e0219854. PMID:31329615

記事内容は個人の意見であり、所属組織を代表しません。また記事内容の正確性については一切保証致しません。本記事は全て各自の責任のもとご利用下さい。

【変更履歴】
2023/1/5 「級内相関係数の概要>平均平方による算出方法」のICC(2,1)、ICC(2,k)の式を修正

目次

• 級内相関係数の概要
  • 6種類の級内相関係数
  • 平均平方による算出方法
• Rによる算出
  • Package
  • データ生成
  • psychによる算出
  • 平均平方による算出
    • Case 1
    • Case 2, Case 3
• 数理的な背景
  • ICC(,k)の式の導出
    • Case 1
    • Case 2, Case 3
  • 平均平方の期待値
    • Case 1
    • Case 2, Case 3
  • 平均平方による算出方法の導出
    • Case 1
    • Case 2, Case 3
  • シミュレーションによる平均平方の期待値の確認
• まとめ
  • （メモ）今後、追加する可能性ある解析

級内相関係数の概要

本記事では、Shrout and Fleiss (1979) で定義されている6種類の級内相関係数（ICC）を扱う。しばしばICC(1,1)、ICC(1,k)、ICC(2,1)、ICC(2,k)、ICC(3,1)、ICC(3,k)と表現される。以後、ICC(1,)をCase 1、ICC(2,)をCase 2、ICC(3,)をCase 3と呼ぶ。

Case 1のときに検者内相関係数、Case 2,3のときには検者間相関係数などと呼ばれたりもする。ICCの解釈等についてはKoo TK and Li MY A (2016)が参考になるかもしれない。

級内相関係数を考える場合に、検査値（測定値）$y_{ij}$を患者ごとのばらつき$p_{i}$、測定ごとのばらつき$r_{i}$、これらの交互作用$pr_{ij}$、誤差$e_{ij}$に分けて考えていくことになる。しかしICC(1,1)、ICC(1,k)、ICC(2,1)、ICC(2,k)、ICC(3,1)、ICC(3,k)は、一元配置分散分析もしくは繰り返しのない二元配置分散分析の状況が多く、これらでは交互作用と誤差を区別することができない。

そのため検査値のモデルを考えるにあたりShrout and Fleiss (1979)のような交互作用を考慮したモデルではなく、Liljequist D, Elfving B and Skavberg Roaldsen K (2019)で用いられているより簡便なモデルを採用した。

6種類の級内相関係数

検査値のモデルと患者によるばらつき、測定によるばらつき、誤差によるばらつきを用いて6種類の級内相関係数を表す。ICC(,1)については比較的理解しやすいが、ICC(,k)のときにどの分散にkの逆数がかかるのかについては後ほど説明する。

記法

• $i$ : 患者（i=1,2,…, n）
• $j$ : 測定回数（Case 1）、測定者（Case 2,3）(j=1,2,…,k)
• $y_{ij}$ : 検査値
• $\mu$ : 平均値
• $p_{i}$：患者による平均からのばらつき
• $r_{i}$：測定による平均からのばらつき
• $e_{ij}$：誤差
• $σ_{p}^2$：$p_{i}$の分散
• $σ_{r}^2$：$r_{i}$の分散
• $σ_{e}^2$：$e_{ij}$の分散
• $\rho$：級内相関係数

Case 1
〇検査値のモデル

$$y_{ij} = \mu + p_{i} + e_{ij}$$

〇ICC(1,1)の式

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+σ_{e}^2}$$

〇ICC(1,k)の式

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+\frac{1}{k}σ_{e}^2}$$

Case 2,3
〇検査値のモデル

$$y_{ij} = \mu + p_{i} + r_{j} + e_{ij}$$

〇ICC(2,1)の式

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+σ_{r}^2+σ_{e}^2}$$

〇ICC(2,k)の式

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+\frac{1}{k}(σ_{r}^2+σ_{e}^2)}$$

〇ICC(3,1)の式

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+σ_{e}^2}$$

〇ICC(3,k)の式

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+\frac{1}{k}σ_{e}^2}$$

平均平方による算出方法

級内相関係数の式を述べたが、得られたデータから分散を推定する必要がある。平均平方から各分散を推定することになり、級内相関係数は以下の式で求めることができる。各分散を推定する式は平均平方の期待値を考えればよいので、それについては後ほど述べる。

記法

• $MSp$ : 患者間の平均平方
• $MSr$ : 測定間の平均平方
• $MSe$ : 誤差の平均平方

〇ICC(1,1)

$$\rho=\frac{MSp-MSe}{MSp}$$

〇ICC(1,k)

$$\rho=\frac{MSp-MSe}{MSp+(k-1)MSe}$$

〇ICC(2,1)

$$\rho=\frac{MSp-MSe}{MSp+(k-1)MSe+\frac{k}{n}(MSr-MSe)}$$

〇ICC(2,k)

$$\rho=\frac{MSp-MSe}{MSp+\frac{1}{n}(MSr-MSe)}$$

〇ICC(3,1)

$$\rho=\frac{MSp-MSe}{MSp+(k-1)MSe}$$

〇ICC(3,k)

$$\rho=\frac{MSp-MSe}{MSp}$$

Rによる算出

数式での説明の前に、Rを用いて級内相関係数を求める方法をまとめる。パッケージを利用した方法（今回はpsychを利用）と平均平方を自ら組み合わせる方法を述べる。

Package

利用するPackageは以下の通り。

• psych

# package
library(psych)

データ生成

まず4人の評価者が10人の患者に対して1回ずつ測定する状況を考えてデータを作成した。これはCase 2もしくはCase 3にあたる。Case 1 とみなす場合には、rater を同一の評価者が同一の患者に対して測定した回数とみれば1人の評価者が10人の患者に対して4回ずつ測定した状況となる。

# データ生成
set.seed(1)

k <- 4
n <- 10

r <- rnorm(k, mean=0, sd=sqrt(2))
p <- rnorm(n, mean=0, sd=sqrt(3))
e <- rnorm(k*n, mean=0, sd=1)

data <- c()

for (j in 1:k){
  for (i in 1:n){
      y <- r[j] + p[i] + e[n*(j-1)+i]
      data <- rbind(data,c(j,i,y))
  }
}

data <- as.data.frame(data)
colnames(data) <- c("rater","patient","y")
data$rater <- as.factor(data$rater)
data$patient <- as.factor(data$patient)

data

##    rater patient           y
## 1      1       1  0.80971565
## 2      1       2 -2.35196601
## 3      1       3 -0.05787785
## 4      1       4  1.33671264
## 5      1       5  0.93256427
## 6      1       6 -0.82098636
## 7      1       7  2.65151969
## 8      1       8  0.57142512
## 9      1       9 -1.88739474
## 10     1      10 -6.71126390
## 11     2       1  1.45026083
## 12     2       2 -1.21751079
## 13     2       3  0.94816726
## 14     2       4  0.06777440
## 15     2       5  0.77884338
## 16     2       6  0.14870424
## 17     2       7  4.23687222
## 18     2       8  0.83215144
## 19     2       9 -0.42863776
## 20     2      10 -3.63006689
## 21     3       1 -1.98809267
## 22     3       2 -3.01784481
## 23     3       3 -0.73179539
## 24     3       4  0.03774519
## 25     3       5  0.91555061
## 26     3       6 -0.94752977
## 27     3       7  1.27220088
## 28     3       8 -0.75989070
## 29     3       9 -1.56081419
## 30     3      10 -4.46106685
## 31     4       1  2.13803625
## 32     4       2  0.12747966
## 33     4       3  3.46490159
## 34     4       4  4.30341657
## 35     4       5  3.14100409
## 36     4       6  2.60822727
## 37     4       7  5.27265542
## 38     4       8  2.31926965
## 39     4       9  1.52116719
## 40     4      10 -2.70926808

psychによる算出

psychパッケージのICCを用いれば、すぐに級内相関係数を算出することができる。ICC()へはマトリックスの形でデータをいれる必要があるため、今回はデータ整形が必要になる。

# データ整形
mtrx <- matrix(data$y, ncol = k, nrow = n)
df <- as.data.frame(mtrx)

print(df)

##             V1         V2          V3         V4
## 1   0.80971565  1.4502608 -1.98809267  2.1380363
## 2  -2.35196601 -1.2175108 -3.01784481  0.1274797
## 3  -0.05787785  0.9481673 -0.73179539  3.4649016
## 4   1.33671264  0.0677744  0.03774519  4.3034166
## 5   0.93256427  0.7788434  0.91555061  3.1410041
## 6  -0.82098636  0.1487042 -0.94752977  2.6082273
## 7   2.65151969  4.2368722  1.27220088  5.2726554
## 8   0.57142512  0.8321514 -0.75989070  2.3192696
## 9  -1.88739474 -0.4286378 -1.56081419  1.5211672
## 10 -6.71126390 -3.6300669 -4.46106685 -2.7092681

# psychによる算出
ICC(df)

6種類の級内相関係数について、信頼区間とともに算出されていることがわかる。

## Call: ICC(x = df)
## 
## Intraclass correlation coefficients 
##                          type  ICC    F df1 df2       p lower bound upper bound
## Single_raters_absolute   ICC1 0.58  6.5   9  30 4.3e-05        0.28        0.85
## Single_random_raters     ICC2 0.61 29.8   9  27 9.3e-12        0.16        0.88
## Single_fixed_raters      ICC3 0.88 29.8   9  27 9.3e-12        0.72        0.96
## Average_raters_absolute ICC1k 0.85  6.5   9  30 4.3e-05        0.60        0.96
## Average_random_raters   ICC2k 0.86 29.8   9  27 9.3e-12        0.42        0.97
## Average_fixed_raters    ICC3k 0.97 29.8   9  27 9.3e-12        0.91        0.99
## 
##  Number of subjects = 10     Number of Judges =  4
## See the help file for a discussion of the other 4 McGraw and Wong estimates,

平均平方による算出

パッケージを利用せずとも平均平方を組み合わせることで級内相関係数は計算することができる。
パッケージを利用した場合と同様の結果が得られていることが確認できる。

Case 1

# 平均平方による算出

anova <- summary(aov(y ~ patient ,data=data))
MS <- anova[[1]][["Mean Sq"]]

MSp <- MS[1]
MSe <- MS[2]

ICC1_1 <- (MSp-MSe)/(MSp + (k-1)*MSe)
ICC1_k <- (MSp-MSe)/MSp

print(c(ICC1_1,ICC1_k))

## [1] 0.5789260 0.8461425

Case 2, Case 3

anova <- summary(aov(y ~ rater + patient ,data=data))
MS <- anova[[1]][["Mean Sq"]]

MSr <- MS[1]
MSp <- MS[2]
MSe <- MS[3]

ICC2_1 <- (MSp-MSe)/(MSp + (k-1)*MSe + (k*(MSr-MSe)/n))
ICC2_k <- (MSp-MSe)/(MSp + (MSr-MSe)/n)

ICC3_1 <- (MSp-MSe)/(MSp + (k-1)*MSe)
ICC3_k <- (MSp-MSe)/MSp

print(c(ICC2_1,ICC2_k,ICC3_1,ICC3_k))

## [1] 0.6109442 0.8626619 0.8779913 0.9664256

数理的な背景

本章ではICC(,1)からICC(,k)の式を求める、また平均平方による計算で級内相関係数を求める式の導出を行う。

ICC(,k)の式の導出

ICC(,k)は測定に関する部分の平均をとっている場合の級内相関係数になる。Case 1、Case 2、Case 3でそれぞれ測定に関する平均をとると以下のような式になる。

Case 1

$$\frac{1}{k}\displaystyle\sum_{j=1}^ky_{ij} = \mu + p_{i} +\frac{1}{k}\displaystyle\sum_{j=1}^ke_{ij}$$

〇ICC(1,k)の式
つまり、$\frac{1}{k}\displaystyle\sum_{j=1}^ke_{ij}$の分散は$\frac{1}{k}σ_{e}^2$になるため、その部分だけ変えてあげればよい。

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+\frac{1}{k}σ_{e}^2}$$

Case 2,3

$$\frac{1}{k}\displaystyle\sum_{j=1}^ky_{ij} = \mu + p_{i} + \frac{1}{k}\displaystyle\sum_{j=1}^kr_{j}+\frac{1}{k}\displaystyle\sum_{j=1}^ke_{ij}$$

〇ICC(2,1)の式

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+σ_{r}^2+σ_{e}^2}$$

〇ICC(2,k)の式
同様に、$\frac{1}{k}\displaystyle\sum_{j=1}^kr_{j}$の分散は$\frac{1}{k}σ_{r}^2$、$\frac{1}{k}\displaystyle\sum_{j=1}^ke_{ij}$の分散は$\frac{1}{k}σ_{e}^2$ になるため、その部分を変える。

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+\frac{1}{k}(σ_{r}^2+σ_{e}^2)}$$

〇ICC(3,k)の式
Case 3では、もともと$σ_{r}^2$は考慮していないため$σ_{e}^2$だけ変更する。

$$\rho=\frac{σ_{p}^2}{σ_{p}^2+\frac{1}{k}σ_{e}^2}$$

今回は測定に関する平均、つまりjについて平均をとっているがその他の平均を考えるような場合でも同様のかたちで整理すればよいと思われる。

平均平方の期待値

ここでは次章に向けて、平均平方の期待値を整理する。

記法

• $ST$ ：全体の平方和
• $Sp$ : 患者間の平方和
• $Sr$ : 測定間の平方和
• $Se$ : 誤差の平方和

また基本的な仮定として、

• $\sum p_{i} = 0$
• $\sum r_{j} = 0$
• $\sum e_{ij} = 0$
• $p_{i}$, $r_{j}$, $e_{ij}$ はそれぞれ独立

をおいている。

Case 1

$$y_{ij} = \mu + p_{i} + e_{ij}$$

$$\begin{aligned} E[ST] &= E[\sum_i\sum_j (Y_{ij} - Y_{..})^2] \cr &= E[\sum_i\sum_j (p_{i}-p_.)^2 + (e_{ij}-e_{..})^2] \cr &= k*E[\sum_i (p_{i}-p_.)^2]+E[\sum_i\sum_j(e_{ij}-e_{..})^2] \cr &= k*(n-1)*σ_p^2+(kn-1)*σ_e^2 \cr \end{aligned}$$

$$\begin{aligned} E[Sp] &= E[\sum_i\sum_j (Y_{i.} - Y_{..})^2] \cr &= k*E[\sum_i((p_{i}-p_.)^2 + (e_{i.}-e_{..})^2)] \cr &= k*(n-1)*σ_p^2+(n-1)*σ_e^2 \cr \end{aligned}$$

$$\begin{aligned} E[Se] &= E[ST-Sp] \cr &= n*(k-1)*σ_e^2 \cr \end{aligned}$$

Case 2,3

$$y_{ij} = \mu + p_{i} + r_{j} + e_{ij}$$

$$\begin{aligned} E[ST] &= E[\sum_i\sum_j (Y_{ij} - Y_{..})^2] \cr &= E[\sum_i\sum_j (p_{i}-p_.)^2+(r_{j}-r_.)^2 + (e_{ij}-e_{..})^2] \cr &= k*E[\sum_i (p_{i}-p_.)^2]+n*E[\sum_j (r_{j}-r_.)^2]+E[\sum_i\sum_j(e_{ij}-e_{..})^2] \cr &= k*(n-1)*σ_p^2+n*(k-1)*σ_r^2+(nk-1)*σ_e^2 \cr \end{aligned}$$

$$\begin{aligned} E[Sp] &= E[\sum_i\sum_j (Y_{i.} - Y_{..})^2] \cr &= k*E[\sum_i((p_{i}-p_.)^2 + (e_{i.}-e_{..})^2)] \cr &= k*(n-1)*σ_p^2+(n-1)*σ_e^2 \cr \end{aligned}$$

$$\begin{aligned} E[Sr] &= E[\sum_i\sum_j (Y_{.j} - Y_{..})^2] \cr &= n*E[\sum_j((r_{j}-r_.)^2 + (e_{.j}-e_{..})^2)] \cr &= n*(k-1)*σ_r^2+(k-1)*σ_e^2 \cr \end{aligned}$$

$$\begin{aligned} E[Se] &= E[ST-Sp-Sr] \cr &= (nk-n-k+1)*σ_e^2 \cr \end{aligned}$$

各分散の推定方法

先ほど求めた平均平方の期待値を用いて分散分析表を作成して、各分散を求める式を整理した。

Case 1

分散分析表

変動因	自由度f	平均平方の期待値
全体	nk-1
患者 (MSp)	n-1	$kσ_{p}^2+σ_e^2$
誤差 (MSe)	n(k-1)	$σ_e^2$

よって、

$$σ_e^2 = MSe$$

$$σ_p^2 = \frac{(MSp-MSe)}{k}$$

として代入すればよい。

Case 2, 3

分散分析表

変動因	自由度f	平均平方の期待値
全体	nk-1
患者 (MSp)	n-1	$kσ_{p}^2+σ_e^2$
測定者 (MSr)	k-1	$nσ_{r}^2+σ_e^2$
誤差 (MSe)	nk-n-k+1	$σ_e^2$

よって、

$$σ_e^2 = MSe$$

$$σ_p^2 = \frac{(MSp-MSe)}{k}$$

$$σ_r^2 = \frac{(MSr-MSe)}{n}$$

として代入すればよい。

シミュレーションによる平均平方の期待値の確認

平均平方の期待値の式があっているかの確認として簡易的なシミュレーションを実施した。$σ_r^2 = 2$, $σ_p^2 = 3$, $σ_e^2 = 1$として、$k=4$, $n=10$で、Case2, 3に該当する平均平方の期待値を10,000回の試行回数の平均値として推定した。おおよそ期待値どおりの結果が得られたことが確認できた。

変動因	自由度f	平均平方の期待値	今回の値で予想される値
全体	nk-1
患者 (MSp)	n-1	$kσ_{p}^2+σ_e^2$	$4*3+1=13$
測定者 (MSr)	k-1	$nσ_{r}^2+σ_e^2$	$10*2+1=21$
誤差 (MSe)	nk-n-k+1	$σ_e^2$	$1$

# 平均平方の期待値のシミュレーション
set.seed(1234)
k <- 4
n <- 10


MS_re <- c()

for (try in 1:10000){
  r <- rnorm(k, mean=0, sd=sqrt(2))
  p <- rnorm(n, mean=0, sd=sqrt(3))
  e <- rnorm(k*n, mean=0, sd=1)

  data <- c()

  for (j in 1:k){
    for (i in 1:n){
        y <- r[j] + p[i] + e[n*(j-1)+i]
        data <- rbind(data,c(j,i,y))
    }
  }
  
  data <- as.data.frame(data)
  colnames(data) <- c("rater","patient","y")
  data$rater <- as.factor(data$rater)
  data$patient <- as.factor(data$patient)
  
  anova <- summary(aov(y ~ rater + patient ,data=data))
  MS <- anova[[1]][["Mean Sq"]]
  
  MS_re <- rbind(MS_re,c(MS[1],MS[2],MS[3]))
}

rbind(c(n*2+1,k*3+1,1),apply(MS_re,2,mean))

##          [,1]     [,2]     [,3]
## [1,] 21.00000 13.00000 1.000000
## [2,] 20.89972 12.89534 1.007152

まとめ

本記事ではRでの級内相関係数の算出方法の紹介および数式による整理を行った。級内相関係数は多くの定義があり、各状況にあわせて使い分ける必要がある指標になる。本記事が級内相関係数がどのような指標なのかを理解する助けになれば幸いである。もし誤りを見つけられた方はご連絡頂きたい。

（メモ）今後、追加する可能性ある解析

• 繰り返しのある2元配置の級内相関係数
• 階層性のある場合の級内相関係数

---

サポートして頂けるとモチベーションに繋がりますのでぜひ宜しくお願いします。データ解析や臨床研究でのご相談があれば、お気軽にTwiiterもしくはメールにてご連絡下さい。

作成者：Masahiro Kondo
作成日：2022/7/30
連絡先：m.kondo1042(at)gmail.com