クラスターRCTとかマルチレベルモデル、他にdesign effectとか測定の一致がどうとかで、グループ内での変量効果的な物を考えることがある。

この、ある種の変量効果の比を、グループ内相関だとかクラスター内相関とか、級内相関だとかICCだとかの名前で呼び、相関係数の一種として考えているみたい。

ところで、いわゆる一般的なピアソンの積率相関係数は[-1,1]であり、（標本相関係数の定義とCauchy–Schwarzの不等式により）その推定値も[-1,1]からはみ出さない性質がある。

ということは、級内相関も当然[-1,1]くらいに理論値も推定値も収まってるはず、くらいの先入観を持ってもおかしくない。*1

本記事では、この先入観について、簡単な例を考えてみよう、というもの。別にその辺の解説記事を否定するものではないから、自分の宗教に合わなければ信じなければよい。

今回対象とするiccについて

一元配置変量モデルにおけるICC(1,k)を考える。最も簡単なタイプで、繰り返し測定のモデル的なやつ。

つまり、観測値 $\displaystyle Y_{i,j}$ が次のモデルに従って生成されると考える。

$\displaystyle Y_{i,j}= \mu + T_{i}+ W_{i,j}$

ここで、 $\displaystyle i =(1,2,\cdots , n),j = (1,2,\cdots ,k)$ 、 $\displaystyle T_{i}\sim N(0, \sigma^2_{T})$ 、 $\displaystyle W_{i,j}\sim N(0, \sigma^2_{W})$ とする。モデル通りの設定なので特に変なことはしていない。

で、ICC(1,k)の定義は以下

$\displaystyle ICC(1,k):= \frac{ \sigma^2_{T} }{ \sigma^2_{T} + \frac{ \sigma^2_{W}}{k} }$

ここから、ICC(1,k)の定義から、次の性質を持つことがわかる。

分散に関わる量の比なので、必ず0以上
最大は1になる

なので、-1から1の値をとる、みたいな解説はこの時点で間違いであると考えてもよい。

相関係数なのに[0,1]なのかよ、と思わないでもないけど、回帰分析の重相関係数は[0,1]なのでまだ許せる。

iccの推定方法

定義の性質は分かった、では推定値の性質はどうか。

例えばRでは、irrパッケージやpsychパッケージでiccを計算できる。

iccが考えているのは変量効果であるため、本来ならば分散成分の推定が必要となる。分散成分は次の方法で推定される。

昔ながらの分散分析モデルで推定する
混合モデルをいわゆるMLEとかで頑張って推定する

前者はいわゆるrepeated measurement anovaな感じの計算であり、irrパッケージはこれしか対応しない。後者は少し新しい方法で、計算コストは少し高いがそれなりに良い推定ができる。現在のpsychパッケージはデフォルトでこちらの方法で計算する。

古典的な方法で計算する推定値の挙動

ここからが本題。「昔ながらの分散分析モデルで推定したicc(1,k)はどのような分布となるか」を考えてみよう、コードは以下

library(psych)
nrow<-10
ncol<-4

T_sd<-1
W_sd<-2

run <- 100
estimated_icc <- rep(0,run)
for (i in 1:run){
  #T_i+W_ijの行列を作成
  idat<-diag(rnorm(nrow,sd=T_sd)) %*% matrix(1,nrow=nrow,ncol=ncol) + matrix(rnorm(nrow*ncol,sd=W_sd),nrow=nrow,ncol=ncol)
  #平均のICC(1,k)を計算
  estimated_icc[i]<-ICC(idat,lmer=FALSE)$results$ICC[4]
}
hist(estimated_icc)