## tl;dr > [[ICC を直感的に理解する]] の主張を、答えの分かっているデータを作って R で確かめます。 > > 1. ICC を 0.7 に決めて作ったデータから、相関・回帰係数・分散分析の 3 とおりで 0.7 を取り出す > 2. ICC を 0 から 1 まで振ると、回帰係数がその値に追随する > > 同じ内容の Stata 版は [[ICC を直感的に理解する（シミュレーション編、Stata版）]] にあります。乱数の種を固定しているので、コピペすれば同じ数字が出ます。 ## ICC 0.7 のデータを作り、3 とおりに取り出す ICC を 0.7 に決めて 100 人分を作り、相関・回帰係数・分散分析の 3 とおりで取り出します。同じ値が出ます。 ```r # 700_icc_kakunin.R ICC 0.7 のデータを作り、3 とおりに取り出す #----------------------------------------- ## 準備 #----------------------------------------- set.seed(398) n <- 100 #----------------------------------------- ## 全体の設定 # 平均、全体の SD、ICC を決める #----------------------------------------- heikin <- 60 sd_zentai <- 10 icc <- 0.7 icc_igai <- 1 - icc #----------------------------------------- ## データを作る # ICC は分散の比として定義されているので、 # SD に配分して使うときにはルートをとる必要がある #----------------------------------------- # 実力は人ごとに一度だけ作る。2 回の測定で共有する jitsuryoku <- heikin + rnorm(n, 0, sd_zentai*sqrt(icc)) # 運 (un) はテストのたびに作り直す un_1 <- rnorm(n, 0, sd_zentai*sqrt(icc_igai)) un_2 <- rnorm(n, 0, sd_zentai*sqrt(icc_igai)) test_1 <- jitsuryoku + un_1 test_2 <- jitsuryoku + un_2 #----------------------------------------- ## 解析 #----------------------------------------- round(sapply(data.frame(jitsuryoku, un_1, un_2, test_1, test_2), function(x) c(mean = mean(x), sd = sd(x), min = min(x), max = max(x))), 1) # ICC である 0.7 が相関として示される # ICC は同じ人を 2 回測った値の相関である（正確には、2 回を入れ替え可能として扱った相関） cor(test_1, test_2) # 回帰でも同じ値 (0.7) が出る # 係数が 1 を下回るので、平均あたりを境に、上の人は下がり下の人は上がる # これで「もとが高い人ほど下がる！」は間違い。効果も変化も入れていない fit <- lm(test_2 ~ test_1) summary(fit)$coefficients # 縮んだと言うなら、比べる相手は 0 ではなく 1 b <- coef(fit)[2] se <- summary(fit)$coefficients[2, 2] tval <- (b - 1) / se c(t = tval, p = 2*pt(-abs(tval), df.residual(fit))) # ICC を一元配置の平均平方から出す。ICC(1,1) にあたる id <- factor(rep(1:n, 2)) test <- c(test_1, test_2) ms <- summary(aov(test ~ id))[[1]][["Mean Sq"]] (ms[1] - ms[2]) / (ms[1] + ms[2]) # 測定は 2 回なので k = 2 #----------------------------------------- ## 変化量をベースラインに回帰しても、同じことをしている #----------------------------------------- # 係数は 0.70 - 1 になる。t 値は上の (b - 1)/se と完全に一致する henka <- test_2 - test_1 summary(lm(henka ~ test_1))$coefficients ``` ``` jitsuryoku un_1 un_2 test_1 test_2 mean 59.5 0.7 0.3 60.2 59.8 sd 8.6 5.1 5.0 9.5 9.6 min 42.5 -11.5 -14.9 41.6 41.2 max 80.0 12.9 12.3 81.7 80.8 > cor(test_1, test_2) [1] 0.6972493 > summary(fit)$coefficients Estimate Std. Error t value Pr(>|t|) (Intercept) 17.5925074 4.43677951 3.965153 1.395321e-04 test_1 0.7014747 0.07284946 9.629099 7.727968e-16 > c(t = tval, p = 2*pt(-abs(tval), df.residual(fit))) t.test_1 p.test_1 -4.097838e+00 8.600319e-05 > (ms[1] - ms[2]) / (ms[1] + ms[2]) [1] 0.6991731 > summary(lm(henka ~ test_1))$coefficients Estimate Std. Error t value Pr(>|t|) (Intercept) 17.5925074 4.43677951 3.965153 1.395321e-04 test_1 -0.2985253 0.07284946 -4.097838 8.600319e-05 ``` 3 とおりの求め方が 0.70 の近くに揃います。 | 求め方 | 値 | |---|---| | `cor(test_1, test_2)` | 0.6972 | | `lm(test_2 ~ test_1)` の係数 | 0.7015 | | 平均平方から計算した ICC | 0.6992 | 設定した 0.70 に対して三桁目が揺れるのは、n = 100 だからです。種を引き直せば ±0.1 ほど動きます。回帰係数が相関とほぼ同じ値になるのは、`test_1` と `test_2` の SD が 9.5 と 9.6 でほぼ揃っているためです。 R には Stata の `icc` にあたる標準の関数がないので、一元配置の平均平方から直接計算しました。この値は Stata の `icc`（one-way）と一致します。同じデータを両方に渡して 0.6991731 と 0.6991732 になることを確認しています。`psych` パッケージを入れているなら `psych::ICC()` の `ICC1` が同じものです。 ### 落とし穴 係数が 1 を下回るので、上の人は下がり下の人は上がります。境目は 58.93 点です。 | test_1 | 予測 test_2 | 差 | |---|---|---| | 40 点 | 45.65 点 | +5.65 | | 58.93 点 | 58.93 点 | 0.00 | | 80 点 | 73.71 点 | −6.29 | **このデータには治療も変化も入っていません。**同じ人を 2 回測っただけで、こうなります。境目が平均の 60.2 点より少し低いのは標本のぶれで、`test_2` の平均が `test_1` より 0.4 点低く出たためです。母集団では境目が平均に一致します。全体を底上げする効果がないからです。 ただし、回帰表の $p < 0.001$ を根拠にはできません。あれは「係数がゼロでない」の検定で、1 回目が 2 回目をいくらかでも予測すれば通ります。縮んだと言いたいなら比べる相手は 1 で、$(b - 1)/\mathrm{se}$ を計算します。 そしてその検定は、変化量をベースラインに回帰したときの検定と完全に一致します。どちらも $t = -4.0978$、$p = 8.6 \times 10^{-5}$ です。変化量の回帰は別の解析に見えて、$\beta - 1$ を計算しているだけだからです。測定誤差があれば $\beta$ は必ず 1 を下回るので、この検定は何も起きていなくても通ります。詳しくは [[平均への回帰と回帰希釈を直感的に理解する]] を参照してください。 ## ICC を 6 条件で振ってみる 点数全体の SD を 10 に固定したまま、実力と運への配分だけを変えます。回帰係数が設定した ICC に追随します。 ```r library(ggplot2); library(dplyr); library(tidyr) set.seed(251) n <- 100 sd_total <- 10 # 点数全体の SD。どの条件でもこれに固定する mean_score <- 50 # 点数の平均 ICCs <- c(1.0, 0.9, 0.8, 0.6, 0.3, 0.0) # 点数の SD を 10 に固定し、実力と運に配分する tsukuru <- function(icc) { sd_jitsuryoku <- sd_total*sqrt(icc) # 実力の個人差 sd_un <- sd_total*sqrt(1 - icc) # 測るたびに変わる分（運） jitsuryoku <- rnorm(n, mean_score, sd_jitsuryoku) data.frame(id = 1:n, icc = icc, kai1 = jitsuryoku + rnorm(n, 0, sd_un), # 1 回目 kai2 = jitsuryoku + rnorm(n, 0, sd_un)) # 2 回目 } d <- bind_rows(lapply(ICCs, tsukuru)) # 回帰係数と相関を出す d %>% group_by(icc) %>% summarise(回帰係数 = coef(lm(kai2 ~ kai1))[2], 相関 = cor(kai1, kai2)) ``` ``` icc 回帰係数 相関 1 0.0 -0.026 -0.026 2 0.3 0.286 0.297 3 0.6 0.595 0.644 4 0.8 0.802 0.702 5 0.9 0.890 0.942 6 1.0 1.000 1.000 ``` 回帰係数と相関は、原理的には同じ値になります。一致しないのは、1 回目と 2 回目の SD が標本では微妙に違ってくるためです。人数を増やせば一致します。 ## 次に読む - [[ICC を直感的に理解する]] — この記事が確かめている理屈 - [[ICC を直感的に理解する（シミュレーション編、Stata版）]] — 同じ内容の Stata 版 - [[平均への回帰と回帰希釈を直感的に理解する]] — 「ベースラインが悪い人ほど改善が大きい」が測定誤差だけで出てしまう話 - [[R - ICC を計算してみよう（評価者間一致研究）]] — 実データから ICC を計算する手順