理学療法士がまとめたノート

統計学備忘録（R言語のメモ）

since2016　ときどきTEXのメモ

級内相関係数：ICC(1, 1)、ICC(1, k)

分散分析 信頼性

級内相関係数 (ICC：Intraclass Correlation Coefficient)

下記の文献のICC(1, 1)、ICC(1, k)をRでやってみます
Shrout PE, Fleiss JL. Intraclass correlations: uses in assessing rater reliability. Psychol Bull. 1979;86(2):420-8.

SST : total sum of squares, sum of squares total
SSA : sum of squares among groups
SSE : sum of squares error
MSA : mean of square
MSE :mean of square

N : 全体のサンプル数
n : 被験者数
k : 繰り返し数（1回～3回）

sampleは、下記のページより
統計学入門−第5章統計学入門−第5章

ICC(1,1)≒0.942、95％信頼区間：0.848  0.984
ICC(1,3)≒0.980、95％信頼区間：0.943  0.995

使用するパッケージ

library(irr)
library(tidyr)

下記のような引数が用意されています

icc(ratings, model = c("oneway", "twoway"), 
    type = c("consistency", "agreement"), 
    unit = c("single", "average"), r0 = 0, conf.level = 0.95)

データセットをRにコピーします

ID <- 1:10
T1 <- c(15,30,34,52,58,69,76,88,91,95)
T2 <- c(10,14,42,38,51,78,88,90,94,87)
T3 <- c(21,38,36,40,42,63,72,84,98,96)
dat1 <- data.frame(ID, T1, T2, T3)

分散分析のために縦のデータに変換

dat2 <- dat1 %>% gather(group, data, -ID)

変数を名義変数へ変換

dat2$group <- as.factor(dat2$group)
dat2$ID <- as.factor(dat2$ID)

ICC(1, 1) 一要因モデル

stripchart(data ~ group, data = dat2, method="jitter",  jitter = 0.05, vertical=TRUE, pch=20, cex=1.5, xlim=c(0.5,3.5), ylim=c(0, 100), xlab="評価3回実施", ylab="被験者10名")

f:id:yoshida931:20210913230228p:plain:w500

一般的な分散分析の散布図、この図では評価（T1, T2, T3）が要因となります。
1人の評価者が被験者 ( n = 10 ) に対して複数回 ( k = 3 ) 評価を実施した時の信頼性を求めます。

icc(dat1[,-1], model ="oneway", type ="consistency", unit = "single")  
# ICC(1) = 0.942
# 95%-Confidence Interval for ICC Population Values:
# 0.848 < ICC < 0.984

分散分析の前に線形回帰について
（要因を被験者としていることに注意）

fit1 <- lm(data ~ factor(ID), data=dat2)
summary(fit1)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    15.333      4.046   3.790  0.00115 ** 
factor(ID)2    12.000      5.721   2.097  0.04887 *  
factor(ID)3    22.000      5.721   3.845  0.00101 ** 
factor(ID)4    28.000      5.721   4.894 8.77e-05 ***
factor(ID)5    35.000      5.721   6.117 5.60e-06 ***
factor(ID)6    54.667      5.721   9.555 6.76e-09 ***
factor(ID)7    63.333      5.721  11.070 5.58e-10 ***
factor(ID)8    72.000      5.721  12.585 5.84e-11 ***
factor(ID)9    79.000      5.721  13.808 1.10e-11 ***
factor(ID)10   77.333      5.721  13.517 1.61e-11 ***

Intercept は当てはめ値（予測値）より
fitted(fit1)[1]
ID=1 の平均値と同値

この回帰分析のモデルは以下のようになる
\boldsymbol{B ^ T}= (b1, b2, b3, \cdots , b10)
\boldsymbol{X ^ T}=(1, x2, x3, \cdots ,x10) ダミー変数
Y = \boldsymbol{B} ^ T*\boldsymbol{X} + ei
ID1: x2～x10=0
ID2: x2=1, x3=\dots=x10=0
ID3: x2=0, x3=1, x3=\dots=x10=0
以下続く・・・

ここでは分散分析表を使います

anova(fit1)
#            Df Sum Sq Mean Sq F value    Pr(>F)    
# factor(ID)  9  22163  2462.5  50.153 9.014e-12 ***
# Residuals  20    982    49.1     

要因が被験者となるので繰り返し数は評価の数
N = 30  #(総サンプル数)    
k = 3   #繰り返し3回（T1, T2, T3）
n = 10  
SSAの自由度
n - 1 = 9    
SSEの自由度
N - n = 20

★ここでは被験者が要因になっている

自由度 平方和 平均平方和 F値
群間 n - 1 SSA MSA=SSA/(n - 1) MSA/ MSE
郡内 N - n SSE MSE=SSE/(N - n)
自由度 平方和 平均平方和 F値
群間 9 22163 2462.5 50.153
郡内 20 982 49.1

注意）一般的な分散分析と異なり被験者が要因となっている。以下の図をイメージしておく。

plot(as.integer(dat2$ID), dat2$data, pch = 20,  xlab = "被験者要因(n : 1～10)", ylab = "k = 3 (T1, T2 ,T3) ")
#　注意）IDが名義のままだとボックスプロットになります

f:id:yoshida931:20210924214059p:plain:w500

ICC1.1の算出方法

MSA = 22163/9    
MSE = 982/20

分散分析モデル
y _ {ij}=μ+a _ {i}+ e _ {ij} \ \ \ ( i=1,2,\cdots,10\ \ \ j=1,2,3 )
a _ {i}=i番目の被験者の効果

全体の分散 V( y _ {ij})
a _ {i}の分散 V(a _ {i})
e _ {ij} の分散を V(e _ {ij})

V( y _ {ij}) = V(a _ {i}) + V(e _ {ij})

ICC(1.1) とは、全体の分散 V( y _ {ij}) に対する分散 V(a _ {i}) の割合
ICC(1.1) ＝\frac{ V(a _ {i}) }{ V(a _ {i}) + V(e _ {ij}) }

MSA ≒ k × V(a _ {i}) + V(e _ {ij})
V(e _ {ij}) = MSE
V(a _ {i}) = \frac{MSA - MSE}{k}
ICC(1.1) = \frac{V(a _ {i})}{V(a _ {i}) +V(e _ {ij})} = \frac{MSA - MSE}{MSA + (k - 1)MSE}

( ICC1.1 <- ( MSA - MSE ) / (MSA + (k - 1) * MSE) )

ICC(1,1)の95%信頼区間の求め方 (分散比の信頼区間より)

F1 <- MSA / MSE
FL1 <- F1 / qf(0.975, n - 1, N - n)
FU1 <- F1 / qf(0.025, n - 1, N - n)
#繰り返し数 : 3 ( k = 3 )
( ICC_1.1_L <- (FL1 - 1) / (FL1 + (k - 1)) )#下限
( ICC_1.1_U <- (FU1 - 1) / (FU1 + (k - 1)) )#上限

ICC(1, k) 一要因モデル 　　

One-way random effects for Case1
1人の評価者が複数の被験者に対して複数回（k = 3）の評価を実施した時の平均値の信頼性に関する指標

k回の平均値が対象となるため誤差の分散が1/k倍となる

ICC(1, 1) より
V(a _ {i}) = \frac{MSA - MSE}{k}
誤差の分散
V(e _ {ij}) = \frac{MSE}{k}

ICC(1.1) = \frac{V(a _ {i})}{V(a _ {i}) +V(e _ {ij})} = \frac{ \frac{MSA - MSE}{k}} {\frac{MSA - MSE}{k} + \frac{MSE}{k} }
= \frac{MSA - MSE}{MSA} 

# unit = averageを使用する
icc(dat1[,-1], model ="oneway", type ="consistency", unit = "average")
# ICC(3) = 0.98     
# 95%-Confidence Interval for ICC Population Values:       
# 0.943 < ICC < 0.995

算出方法

(ICC_1.k <- ( MSA - MSE ) / MSA)
# ICC(1,k)の95%信頼区間の求め方 (分散比の信頼区間より)
(ICC_1.k_L <- (FL1 - 1) / FL1)#下限
(ICC_1.k_U <- (FU1 - 1) / FU1)#上限