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

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

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

 ICC(2, 1) 二要因モデル  

Two-way random effects for Case2

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

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

平均平方和
 MSA : mean of square(評価者)
 MSB : mean of square (被験者)
 MSE : mean of square

 N : 全体のサンプル数
 n : 被験者数
 k : 繰り返し数(評価者A、B、C)

使用するパッケージ

library(irr)
library(tidyr)      

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

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

sampleは、下記のページより
統計学入門−第5章統計学入門−第5章
このページでもとめる値

ICC(2,1)0.94295%信頼区間:0.845  0.984
ICC(2,3)0.98095%信頼区間:0.942  0.995

データセットを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)

以下のようなグラフをイメージして進めます

plot(as.integer(dat2$ID), dat2$data, pch=19, col=c("red", "green", "blue"), xlab="被験者要因(i:1~10)", ylab="評価 (A,B,C) ")
abline(h=mean(dat2$data))
text(x=1.8,y=63,"全体平均=59.7")
legend("topleft", legend = c("A", "B", "C"), col = c("red", "green", "blue"), pch =19)

f:id:yoshida931:20210912221407p:plain:w500

要因:評価者被験者
評価者は変量モデル:多くの評価者からたまたま選ばれた3人
複数の評価者 (k = A, B, C ) が複数の被験者 ( n = 1, 2, … , 10 ) に評価したときの評価者間( A, B, C ) の信頼性

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

fit2 <- lm(data ~ group + factor(ID)  , data = dat2)
summary(fit2)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   16.467      4.625   3.560 0.002235 ** 
groupT2       -1.600      3.270  -0.489 0.630571    
groupT3       -1.800      3.270  -0.550 0.588805    
ID2           12.000      5.971   2.010 0.059684 .  
ID3           22.000      5.971   3.685 0.001696 ** 
ID4           28.000      5.971   4.690 0.000183 ***
ID5           35.000      5.971   5.862 1.50e-05 ***
ID6           54.667      5.971   9.156 3.41e-08 ***
ID7           63.333      5.971  10.607 3.58e-09 ***
ID8           72.000      5.971  12.059 4.66e-10 ***
ID9           79.000      5.971  13.231 1.03e-10 ***
ID10          77.333      5.971  12.952 1.46e-10 ***

Intercept は当てはめ値(予測値)より
fitted(fit2)[1]

分散分析

fit2 <- lm(data ~ ID + group, data = dat2)    
anova(fit2)

Response: data
          Df  Sum Sq Mean Sq F value   Pr(>F)    
group      2    19.5    9.73   0.182   0.8351    
ID         9 22162.7 2462.52  46.051 1.33e-10 ***
Residuals 18   962.5   53.47                               

group = 評価者 (k:A, B, C )
ID = 被験者 ( n:1, 2, … , 10 )

自由度 平方和 平均平方和 F値
評価者 k - 1 SSA MSA=SSA/(k - 1 ) MSA/ MSE
被験者 n - 1 SSB MSB=SSB/(n - 1) MSB/ MSE
残差 (k - 1)(n - 1) SSE MSE=SSE/(k - 1)(n - 1)
全変動 k*n - 1 SST  
icc(dat1[,-1], model ="twoway", type ="agreement", unit = "single")
# ICC(A,1) = 0.942
# 95%-Confidence Interval for ICC Population Values:
# 0.844 < ICC < 0.984

算出方法
分散分析モデル
 y _ {ij}=μ+a _ {i}+b _ {j}+(ab) _ {ij} + e _ {ij}  \   \   \   (i=1, 2, 3 \   \   \   j=1,2,\cdots,10)
 a _ {i} ;評価者の効果
 b _ {j} ;被験者の効果
 (ab) _ {ij} ;被験者 jと評価者 iの交互作用

 ICC(2.1) とは、全体の分散に対する b _ {j}の分散の割合
この値が大きくなれば、評価者間の分散が小さい(信頼性がある)と判断する

 α _ {i}の分散=  V(a _ {i})
 b _ {i}の分散= V(b _ {i})
 (ab) _ {ij} e _ {ij} の分散= V(e _ {ij})

 ICC(2.1) =\frac{V(b_ {j})}{V(a _ {i}) + V(b _ {j}) + V(e _ {ij})}

 MSA ≒ k×V(a _ {i}) + MSE
 MSB ≒ n×V(b _ {i}) + MSE

 V(a _ {i}) = \frac{MSA - MSE}{k}
 V(b _ {i}) = \frac{MSB - MSE}{n}

  MSE = V(e _ {ij})

分散分析表より
k = 3
n = 10
MSA = 9.73
MSB = 2462.52
MSE = 53.47

va <- (MSA - MSE) / 10 
vb <- (MSB - MSE) / 3
(ICC_2.1 <- vb / (va + vb +MSE))

95%信頼区間

Fj <- MSA / MSE
c <- (k - 1)*(n - 1)*(k*ICC_2.1*Fj + n*(1 + (k - 1)*ICC_2.1) - k*ICC_2.1)^2
d <- (n - 1)*k^2*ICC_2.1^2*Fj^2+(n*(1+ (k - 1)*ICC_2.1) - k*ICC_2.1)^2
(FL2 <- qf(0.975, n-1, round(c/d,0)))
(FU2 <- qf(0.975, round(c/d,0), n-1))
(ICC_2.1_L <- (n*(MSB - FL2*MSE)) / (FL2*(k*MSA+(n*k - n - k)*MSE)+n*MSB))#下限
(ICC_2.1_U <- n*(FU2*MSB - MSE) / ((k*MSA + (n*k - n- k)*MSE) + n*FU2*MSB))#上限

 ICC(2, k) 二要因モデル  

k人の評価者が複数の被験者に対して1回実施した時の平均値の信頼性に関する指標

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

icc(dat1[,-1], model ="twoway", type ="agreement", unit = "average")
# ICC(A,3) = 0.98
#  95%-Confidence Interval for ICC Population Values:
#  0.942 < ICC < 0.995

算出方法

k = 3 #評価者数K1  
n = 10 #被験者数K2   
MSA =  9.73
MSB = 2462.52
MSE = 53.47 

ICC(2,k)の分散

va_k <- (MSA - MSE/3) / (10*3)  #評価者の分散
vb_k <- (MSB - MSE/3) / 3  #被験者の分散
(ICC_2.k <- vb_k / (va_k + vb_k +MSE/3))

95%信頼区間

(ICC_2.k_L <- (k*ICC_2.1_L / (1+(k-1)*ICC_2.1_L)))    #下限
(ICC_2.k_U <- (k*ICC_2.1_U / (1+(k-1)*ICC_2.1_U)))   #上限