## tl;dr > [[平均への回帰と回帰希釈を直感的に理解する]] の主張を、答えの分かっているデータを作って確かめます。確かめるのは次の 4 点です。 > > 1. 効果を全員一律に設定しても、「下位ほど伸びた」が出てくる > 2. そのとき観測される回帰係数は 1.00 ではなく ICC になる > 3. $\beta_{\text{観測}} = \beta_{\text{真}} \times \mathrm{ICC}$ が数値として成り立つ > 4. ICC で割り戻せば $\beta_{\text{真}}$ に戻る > > コードは R をメインにし、同じ内容の Stata 版を末尾に置きました。乱数の種を固定しているので、コピペすれば同じ数字が出ます。 理屈のほうは [[平均への回帰と回帰希釈を直感的に理解する]] に、その土台になる ICC は [[ICC を直感的に理解する]] に書きました。このノートはコードと数値の検証だけを扱います。 ## シミュレーションで何をしているか 現実のデータでは「下位ほど伸びた」が本物か見かけかを判断できません。答えを知らないからです。そこで、**答えのほうをこちらで決めてしまいます。** シナリオA では「講習の効果は全員一律 +10 点。成績によって効き方が変わることは一切ない」と決めました。これは期待でも仮定でもなく、コードにそう書いた以上、この条件では動かせない事実です。 ```r kai2_A <- jitsuryoku + 10 + rnorm(n, 0, 5) # ↑ # 実力が 40 点でも 80 点でも足されるのは同じ 10。 # jitsuryoku に応じて上乗せ量が変わる項はどこにもない ``` ここで「下位ほど伸びた」が出たら、それは 100% 見かけです。本物かどうかを議論する余地がありません。真実が分かるのは、私たちがそう作ったからです。**見かけがどれくらいの大きさで出るかを測っておいて、その物差しを持って現実のデータへ戻ります。** シナリオB は逆に、「本当に下位ほど効く」と決めた条件です。 ```r kai2_B <- 60 + 0.50*(jitsuryoku - 60) + 10 + rnorm(n, 0, 5) ``` 真ん中の項が、クラス平均 60 点からのズレを半分に縮めます。実力 40 点の生徒は 60 − 10 + 10 = 60 点、実力 80 点の生徒は 60 + 10 + 10 = 80 点に着地するので、伸びはそれぞれ +20 点と 0 点になります。実力が 20 点違う 2 人の差が 2 回目には 10 点差になるので、$\beta_{\text{真}} = 0.50$ です。 ## 第1部: 模試の点数 ### 設定 | 要素 | 設定 | コード上の記号 | |---|---|---| | 真の実力 | 平均 60 点、SD 10 点 | `jitsuryoku` | | 1 回のテストのたまたま | SD 5 点 | `rnorm(n, 0, 5)` | | 夏期講習の効果 | 全員一律に +10 点 | シナリオA | | 下位ほど効く効果 | ズレを半分に縮める | シナリオB、$\beta_{\text{真}} = 0.50$ | 2 つの SD から ICC が決まります。 $ \mathrm{ICC} = \frac{\sigma^2_{\text{実力}}}{\sigma^2_{\text{実力}} + \sigma^2_{\text{たまたま}}} = \frac{10^2}{10^2 + 5^2} = \frac{100}{125} = 0.80 $ 分母の 125 は答案に実際に現れる点数の分散です。観測される点数の SD は $\sqrt{125} = 11.18$ 点になるはずで、この予測が当たるかを最初に確かめます。 ### コード ```r # 400_part1.R 模試の点数で平均への回帰を体感する set.seed(20260814) n <- 5000 jitsuryoku <- rnorm(n, mean = 60, sd = 10) # 真の実力 kai1 <- jitsuryoku + rnorm(n, 0, 5) # 1 回目の答案 kai2_A <- jitsuryoku + 10 + rnorm(n, 0, 5) # シナリオA: 一律 +10 点 kai2_B <- 60 + 0.50*(jitsuryoku - 60) + 10 + rnorm(n, 0, 5) # シナリオB nobi_A <- kai2_A - kai1 nobi_B <- kai2_B - kai1 # 予測（実力 SD 10 点、答案 SD 11.18 点）が当たっているか c(sd_jitsuryoku = sd(jitsuryoku), sd_kai1 = sd(kai1)) # 1 回目の点数帯ごとの伸び band <- cut(kai1, breaks = c(-Inf, 45, 55, 65, 75, Inf), labels = c("<45", "45-55", "55-65", "65-75", ">=75")) data.frame( band = levels(band), n = as.vector(table(band)), kai1 = round(tapply(kai1, band, mean), 2), nobi_A = round(tapply(nobi_A, band, mean), 2), nobi_B = round(tapply(nobi_B, band, mean), 2), row.names = NULL ) summary(lm(kai2_A ~ kai1)) summary(lm(kai2_B ~ kai1)) ``` ### 結果1: ばらつきの予測が当たったか ``` sd_jitsuryoku sd_kai1 9.904 11.112 ``` 乱数で 5000 人を発生させると実力の SD は 10 点ちょうどにならず、今回は 9.904 点でした。実際に出たこの値から計算し直します。 $ \sqrt{9.904^2 + 5^2} = \sqrt{98.1 + 25} = \sqrt{123.1} = 11.095 \text{ 点} $ 観測された 11.112 点とほぼ一致します。実力のばらつきに運のばらつきが分散として上乗せされる関係が、数値として成り立っています。 ### 結果2: 点数帯ごとの伸び ``` band n kai1 nobi_A nobi_B 1 <45 428 40.15 14.35 21.87 2 45-55 1217 50.62 12.03 15.60 3 55-65 1692 60.09 10.14 10.11 4 65-75 1195 69.42 8.29 4.14 5 >=75 468 80.00 5.51 -2.10 ``` `nobi_A` は講習の効果を全員一律 +10 点に設定したほうです。45 点未満の生徒は平均 14.35 点上がり、75 点以上の生徒は 5.51 点しか上がっていません。**2 倍以上の差**がついています。この列だけを見せられたら「下位者に特に効く」と結論したくなりますが、設定は全員 +10 点なので、差は全部偽物です。 `nobi_B` は本物の収束がある条件で、下位 +21.87 点、上位 −2.10 点とさらに大きく開きます。**「下位ほど伸びた」という定性的な観察は A でも B でも同じように出るので、この表を眺めても区別できません。**区別できるのは回帰係数の数値だけです。 対照群のない前後比較も同じ理由で壊れます。 ```r cut20 <- quantile(kai1, 0.20) mean(nobi_A[kai1 <= cut20]) ``` ``` 1回目 50.7 点以下の生徒の平均上昇 = 13.31 点（真の効果は全員一律 10 点） ``` 「成績不振の生徒だけに講習をしたら 13.3 点上がりました」という報告になります。3.3 点は水増しで、効果がゼロでも同じ理屈で上がって見えます。 ### 結果3: 回帰係数は 1.00 ではなく ICC になる ``` === シナリオA === Estimate Std. Error t value Pr(>|t|) (Intercept) 22.9183659 0.527714388 43.42949 0 kai1 0.7866751 0.008624135 91.21786 0 R2 = 0.6247 === シナリオB === Estimate Std. Error t value Pr(>|t|) (Intercept) 46.3216850 0.419470477 110.42895 0 kai1 0.3944488 0.006855166 57.54036 0 R2 = 0.3985 ``` シナリオA は 0.787 です。$\beta_{\text{真}} = 1$ に設定したのに 1.00 が出ません。標本 ICC は 0.7944 で、これとほぼ一致します。シナリオB は 0.394 で、$0.50 \times 0.80 = 0.40$ にほぼ一致します。 | | $\beta_{\text{真}}$（設定値） | 理論値 $\beta_{\text{真}} \times \mathrm{ICC}$ | 観測された回帰係数 | |---|---|---|---| | A: 下位ほど効く効果ゼロ | 1.00 | 0.794 | 0.787 | | B: 本当に下位ほど効く | 0.50 | 0.397 | 0.394 | $\beta_{\text{観測}} = \beta_{\text{真}} \times \mathrm{ICC}$ が数値で確認できました。 ### 結果4: 係数 0.79 は「上位の人が下がる」ではない 係数が 1 を下回ったことを「下位は上がり、上位は下がる」と読むと後半が外れます。回帰直線に値を入れて確かめます。 ```r b <- coef(lm(kai2_A ~ kai1)) for (x in c(40, 60, 80, 100)) cat(x, "→", round(b[1] + b[2]*x, 2), "\n") b[1] / (1 - b[2]) # 予測が下向きに転じる境界 max(kai1) sum(nobi_A < 0 & kai1 >= 75); sum(kai1 >= 75) ``` ``` kai1 = 40 → 予測 kai2 = 54.39 伸び +14.39 kai1 = 60 → 予測 kai2 = 70.12 伸び +10.12 kai1 = 80 → 予測 kai2 = 85.85 伸び +5.85 kai1 = 100 → 予測 kai2 = 101.59 伸び +1.59 予測が下がる境界 x* = 107.4 実データの最大 kai1 = 101.74 上位帯(>=75) のうち実際に下がった人: 89 / 468 ``` **5000 人全員が上がります。**平均へ寄っているのはズレのほうです。 ``` 1回目 40 点(平均 -20.17) → 2回目は平均 -15.87 比 0.787 1回目 60 点(平均 -0.17) → 2回目は平均 -0.14 比 0.787 1回目 80 点(平均 +19.83) → 2回目は平均 +15.60 比 0.787 ``` 係数は差の情報を持ち、切片は水準の情報を持ちます。上位の人が水準でも下がるかは、係数の不足分が全体の底上げを上回るかで決まります。シナリオB（係数 0.394、切片 46.32）では境界が 76.5 点まで下がるので、上位帯の伸びが −2.10 と実際に減少へ転じます。 係数だけから必ず言えるのは、伸びとベースラインの関係だけです。 ```r coef(lm(nobi_A ~ kai1))[2] # -0.2133 coef(lm(kai2_A ~ kai1))[2] - 1 # -0.2133 ``` 恒等的に「係数 − 1」なので、係数が 1 を下回れば低い人ほど伸びが大きく出ます。 ### 結果5: 平均への回帰は「引き戻す力」ではない 全員が平均へ吸い寄せられるのなら、測るたびに散らばりが縮んでいずれ全員 60 点になるはずです。実力を固定したまま測定を 6 回繰り返します。 ```r x <- kai1 for (r in 1:6) { x <- jitsuryoku + 10*r + rnorm(n, 0, 5); cat(r, sd(x), "\n") } ``` ``` kai1 の SD = 11.112 1回目 11.141 2回目 11.101 3回目 11.154 4回目 11.124 5回目 11.135 6回目 11.030 ``` **縮みません。**毎回あたらしい運が乗るので、分布の幅は保たれます。縮むのは「1 回目に 80 点だった人たち」という条件つきの平均だけです。 さらに、時間を逆にしても同じ係数が出ます。 ```r coef(lm(kai2_A ~ kai1))[2] # 0.7867 coef(lm(kai1 ~ kai2_A))[2] # 0.7942 mean(kai1[kai2_A >= quantile(kai2_A, .75)]) # 71.48（全体平均 60.17） ``` 「2 回目に高かった人は、1 回目には平均寄りだった」も同じ強さで成り立ちます。過去に向かっても回帰するのですから、人を平均へ引き戻す力ではありません。**ノイズを含む変数で人を選んだこと自体の帰結**です。ここを取り違えると、「治療によって差が縮まった」という因果の話に滑ります。 ### 2 段階に分けて確かめる なぜ ICC 倍になるのかは、差を 2 回続けて絞ると考えると分かります。 ``` 1回目の点数の差 真の実力の差 2回目の点数の差 10 点 ─段階1→ 8 点 ─段階2→ 8 点 ICC 倍 β真 倍 ``` 段階1 は「1 回目の点数差のうち、本物の実力差はどれだけか」、段階2 は「実力差が 2 回目にどれだけ引き継がれるか」です。シミュレーションでは真の実力 `jitsuryoku` を持っているので、この 2 段階を別々に計算できます。 ```r coef(lm(jitsuryoku ~ kai1))[2] # 段階1 coef(lm(kai2_A ~ jitsuryoku))[2] # 段階2 coef(lm(kai2_A ~ kai1))[2] # 2 段階まとめて（実データで見えるのはこれだけ） ``` ``` 段階1 0.7943 ← 標本 ICC 0.7944 とほぼ同じ 段階2 0.9961 ← β真 = 1.00 とほぼ同じ 掛け算 0.7943 × 0.9961 = 0.7912 まとめて 0.7867 ← 掛け算の 0.7912 とほぼ同じ ``` 掛け算の 0.7912 と直接回帰の 0.7867 がわずかにずれるのは、有限標本で誤差項が実力とぴったり無相関にならないためです。理論上は完全に一致します。 実データでは真の実力を測れないので、見えるのは掛け算の結果だけです。だから $\beta_{\text{観測}}$ が小さく出たとき、それが段階1 のせい（測定誤差）なのか段階2 のせい（本物の収束）なのかを、ICC を使って切り分けます。 ### 共分散から導く 回帰係数の定義式からも同じ結果が出ます。記号を置きます。 | 記号 | 中身 | |---|---| | $T$ | 真の実力 | | $e_0, e_1$ | 1 回目・2 回目の測定誤差 | | $x = T + e_0$ | 観測された 1 回目の点数 | | $y = \beta_{\text{真}} T + \text{定数} + e_1$ | 観測された 2 回目の点数 | 分母は実力のばらつきに誤差のばらつきが上乗せされて $\mathrm{Var}(x) = \sigma^2_T + \sigma^2_e$ になります。分子のほうは、$e_0$ と $e_1$ が互いにも $T$ とも無関係なので共分散に何も残さず、$\mathrm{Cov}(x, y) = \beta_{\text{真}} \sigma^2_T$ だけが残ります。 $ \beta_{\text{観測}} = \frac{\beta_{\text{真}}\,\sigma^2_T}{\sigma^2_T + \sigma^2_e} = \beta_{\text{真}} \times \frac{\sigma^2_T}{\sigma^2_T + \sigma^2_e} = \beta_{\text{真}} \times \mathrm{ICC} $ 一文にすると、**測定誤差は分母だけを膨らませて分子には効かない**ので、割り算の答えが ICC 倍に縮みます。数字を取り出して確かめます。 ```r cov(kai1, kai2_A); var(jitsuryoku); var(kai1) ``` ``` Cov(kai1, kai2_A) = 97.137 Var(jitsuryoku) = 98.096 Var(kai1) = 123.478 回帰係数 = 97.137 / 123.478 = 0.7867 ICC = 98.096 / 123.478 = 0.7944 ``` シナリオA では $\beta_{\text{真}} = 1$ なので、分子の共分散 97.137 は $\sigma^2_T$ の 98.096 と一致するはずで、実際ほぼ同じです。分母の 123.478 が実力の 98.096 に誤差のぶん約 25 を足した値である点も見てください。この上乗せが回帰係数を縮めています。 ### ものさしを変えると基準が壊れる 「回帰係数 = ICC」という基準は、2 回の測定が同じものさしであることを前提にします。実力も測定誤差も講習の効果もそのままにして、2 回目の点差の開き方だけを変えます。 ```r kai2_C <- 60 + 1.2*(jitsuryoku - 60) + 10 + rnorm(n, 0, 6) # 点差が 1.2 倍開く kai2_D <- 40 + 1.5*(jitsuryoku - 60) + rnorm(n, 0, 7.5) # 1.5 倍開き、平均も 40 点へ ``` | 2 回目の試験 | 点差の開き方 | 観測される回帰係数 | |---|---|---| | 同じ形式の模試 | 1.0 倍 | 0.787 | | 期末試験 | 1.2 倍 | 0.960 | | 本番 | 1.5 倍 | 1.186 | 学力も誤差も効果も一切変えていません。動かしたのは 2 回目の問題の作りだけで、それだけで回帰係数は 0.79 から 1.19 まで動きます。一般に 2 回目の真の値のスケールが $k$ 倍になると $ \beta_{\text{観測}} = k \times \beta_{\text{真}} \times \mathrm{ICC} $ となり、$k$ が未知である限り ICC を基準に使えません。 ### 偏差値化すると収束が消える スケールの違いは標準化で消せますが、**探している収束も一緒に消えます。** ```r z <- function(x) (x - mean(x)) / sd(x) * 10 + 50 # 偏差値 coef(lm(z(kai2_B) ~ z(kai1)))[2] ``` ``` SD: kai1 11.11 kai2_A 11.06 kai2_B 6.94 素点 偏差値化後 シナリオA 0.7867 0.7904 シナリオB 0.3944 0.6312 基準線(ICC) 0.79 0.79 ``` 本物の収束とは「2 回目にばらつきが縮むこと」そのものです。実際 `kai2_B` の SD は 6.94 で、`kai1` の 11.11 から縮んでいます。偏差値化はそのばらつきを強制的に元の幅へ引き伸ばす操作なので、探しているものを自分で潰しにいくことになります。基準線からの距離は 0.39 から 0.16 へ、半分以下になりました。 GLS は同じ装置で同じ指標を 2 回測るので、ものさしの条件はもともと満たされています。実測値のまま扱ってかまいません。逆に、**各時点で GLS を標準化してはいけません。**上の表と同じことが起きます。 ## 第2部: GLS に置き換える TAVI 後 1 年の GLS をベースラインの GLS に回帰したら回帰係数が 0.42 だった、という状況を想定します。この 0.42 が本物の収束なのか測定誤差の産物なのかを判定します。 ### 設定 | 要素 | 設定 | |---|---| | 真のベースライン GLS（絶対値） | 平均 12%、SD 4% | | GLS の再現性 ICC | 0.85（文献値） | | TAVI の一律の改善 | +3.0% | | シナリオB の真の回帰係数 | 0.50 | 測定誤差の SD は ICC から逆算します。 $ \sigma_e = \sigma_{\text{真}} \sqrt{\frac{1 - \mathrm{ICC}}{\mathrm{ICC}}} = 4 \times \sqrt{\frac{0.15}{0.85}} = 1.680 $ シナリオB の理論値は $0.50 \times 0.85 = 0.425$ になります。判定したい実測値 0.42 とほぼ一致するので、この条件は「本物の収束があるとしたら、こう見えるはずだ」という参照になります。 ### コード ```r # 401_part2.R set.seed(20260814) n <- 200000 mu <- 12 # 真のベースライン GLS（絶対値, %）の平均 sdt <- 4 # 個人間 SD icc <- 0.85 # GLS 測定の再現性 delta <- 3.0 # TAVI による一律の改善量（%） btrue <- 0.50 # シナリオB の真の回帰係数 sde <- sdt * sqrt((1 - icc) / icc) # ICC = sdt^2/(sdt^2+sde^2) を sde について解く shin0 <- rnorm(n, mu, sdt) shin1_A <- shin0 + delta # 一律改善のみ shin1_B <- mu + btrue*(shin0 - mu) + delta # ズレを btrue 倍に縮める kansoku0 <- shin0 + rnorm(n, 0, sde) kansoku1_A <- shin1_A + rnorm(n, 0, sde) kansoku1_B <- shin1_B + rnorm(n, 0, sde) summary(lm(shin1_A ~ shin0)) # 真値どうし（誤差ゼロ） summary(lm(kansoku1_A ~ kansoku0)) # 観測値どうし・シナリオA summary(lm(kansoku1_B ~ kansoku0)) # 観測値どうし・シナリオB summary(lm(I(kansoku1_A - kansoku0) ~ kansoku0)) # 変化量をベースラインに回帰 summary(lm(I(kansoku1_B - kansoku0) ~ kansoku0)) ``` ### 結果 ``` === 真値どうし・シナリオA（誤差ゼロ）=== Estimate Std. Error t value Pr(>|t|) (Intercept) 3 3.495554e-16 8.582331e+15 0 shin0 1 2.764533e-17 3.617247e+16 0 === 観測値どうし・シナリオA === Estimate Std. Error t value Pr(>|t|) (Intercept) 4.7633676 0.015048863 316.5268 0 kansoku0 0.8526799 0.001179921 722.6584 0 R2 = 0.7231 === 観測値どうし・シナリオB === Estimate Std. Error t value Pr(>|t|) (Intercept) 9.9202533 0.0121868532 814.0127 0 kansoku0 0.4238325 0.0009555224 443.5610 0 R2 = 0.4959 ``` 誤差をゼロにすると回帰係数はきっちり 1.000 で、$\beta_{\text{真}}$ の設定値そのものが出ます。誤差を入れた途端に 0.853 まで落ちます。 | シナリオ | $\beta_{\text{真}}$ | 理論値 $\beta_{\text{真}} \times \mathrm{ICC}$ | 観測された回帰係数 | |---|---|---|---| | A: 一律改善のみ | 1.00 | 0.850 | 0.853 | | B: 真の収束あり | 0.50 | 0.425 | 0.424 | ### 変化量をベースラインに回帰してはいけない 「改善量 ~ ベースライン」という解析はよく見かけますが、前後比較よりさらに大きく偏ります。 ``` === 変化量 ~ ベースライン。シナリオA（効果は一律）=== Estimate Std. Error t value Pr(>|t|) (Intercept) 4.7633676 0.015048863 316.5268 0 kansoku0 -0.1473201 0.001179921 -124.8559 0 === 変化量 ~ ベースライン。シナリオB === Estimate Std. Error t value Pr(>|t|) (Intercept) 9.9202533 0.0121868532 814.0127 0 kansoku0 -0.5761675 0.0009555224 -602.9869 0 ``` シナリオA は「下位ほど効く効果ゼロ」に設定した条件なのに、係数は −0.147 で $p < 0.001$ です。理論値 $\mathrm{ICC} - 1 = -0.15$ とぴったり一致しています。シナリオB は −0.576 で、理論値 $0.50 \times 0.85 - 1 = -0.575$ に一致します。 $ \beta_{\text{変化量}} = \beta_{\text{真}} \times \mathrm{ICC} - 1 $ $\beta_{\text{真}} = 1$（ベースライン依存性ゼロ）でも係数は $-(1-\mathrm{ICC})$ という負の値になり、n が大きければ必ず有意になります。さらに両辺に同じ測定値 $x$ が入る数学的な連動（mathematical coupling）もあるため、この解析は二重に汚染されています。**改善量をベースラインに回帰した有意な負の係数は、それ自体では何の証拠にもなりません。** ### 希釈補正 ICC で割り戻せば真の回帰係数に戻ります。 ```r coef(lm(kansoku1_A ~ kansoku0))[2] / icc coef(lm(kansoku1_B ~ kansoku0))[2] / icc ``` ``` シナリオA: 0.8527 / 0.85 = 1.0032 （真値 1.00） シナリオB: 0.4238 / 0.85 = 0.4986 （真値 0.50） ``` 信頼区間まで欲しいときは Stata の `eivreg`（errors-in-variables regression）が使えます。末尾の Stata 版に出力を載せました。R で本格的にやるなら `simex` パッケージ（シミュレーション外挿法）がありますが、ICC が既知で単純な線形回帰なら上の割り算で十分です。 ### n = 150 でも区別できるか 実データは 20 万例ではありません。n = 150 の標本を 1000 回取り直します。 ```r set.seed(20260814) res <- replicate(1000, { i <- sample.int(n, 150) c(coef(lm(kansoku1_A[i] ~ kansoku0[i]))[2], coef(lm(kansoku1_B[i] ~ kansoku0[i]))[2]) }) ``` ``` シナリオA: 平均 0.852 2.5-97.5%tile 0.764 〜 0.933 シナリオB: 平均 0.424 2.5-97.5%tile 0.360 〜 0.492 ``` 2 つの分布はまったく重なりません。実測 0.42 はシナリオB の中心にあり、シナリオA の範囲（0.764 以上）からは遠く外れています。n = 150 でも判定できます。 ### 感度分析と損益分岐 ICC 判定は ICC の値に乗っているので、ICC を振って結論が変わらないかを確かめます。 ```r for (i in c(0.75, 0.80, 0.85, 0.90, 0.95)) { cat("ICC =", i, "→ 補正後の回帰係数 =", round(0.42 / i, 3), "\n") } ``` ``` ICC = 0.75 → 0.560 ICC = 0.80 → 0.525 ICC = 0.85 → 0.494 ICC = 0.90 → 0.467 ICC = 0.95 → 0.442 ``` どの ICC を採っても 0.44〜0.56 の範囲で、1.00 から大きく離れています。結論は ICC の選び方に依存しません。 査読者に「GLS の再現性がもっと悪ければ 0.42 は全部誤差で説明できる」と言われたときは逆算して返します。$\beta_{\text{観測}} = 1 \times \mathrm{ICC}$ が成り立つには $\mathrm{ICC} = 0.42$ でなければならず、報告されている 0.85 前後とはかけ離れています。その水準の再現性なら GLS は臨床指標として使い物になりません。 ### 誤差が相関している場合 ここまでベースラインと 1 年後の測定誤差を独立に発生させてきました。同じ読影者が両方を続けて測ると、誤差に正の相関 $\rho$ が生じます。このとき $ \beta_{\text{観測}} = \frac{\beta_{\text{真}} \sigma^2_{\text{真}} + \rho\,\sigma^2_e}{\sigma^2_{\text{真}} + \sigma^2_e} $ となります。シナリオA（$\beta_{\text{真}} = 1$）で $\rho$ を振ります。 ```r for (rho in c(0, 0.3, 0.6)) { sh <- rnorm(n, 0, sde) e0 <- sqrt(rho)*sh + sqrt(1-rho)*rnorm(n, 0, sde) e1 <- sqrt(rho)*sh + sqrt(1-rho)*rnorm(n, 0, sde) cat("rho =", rho, "→", round(coef(lm(I(shin1_A+e1) ~ I(shin0+e0)))[2], 4), "\n") } ``` ``` rho = 0.0 → 0.8497 rho = 0.3 → 0.8941 rho = 0.6 → 0.9406 ``` $\rho > 0$ なら回帰係数は**大きい**方へ動き、希釈は弱くなります。つまり実測 0.42 という小さい値が観測されたことは、この方向のバイアスを考えるとむしろ保守的です。論法にとって有利な向きなので、突かれたら説明できるようにしておくとよいです。 ## この論法の弱点 ### ICC は集団に依存する ICC は「個人間の分散 ÷ 全分散」なので、**対象集団のばらつきが小さいほど下がります**。文献の ICC が幅広い患者を含むコホートで算出されたもので、自分の研究が重症 AS 患者に限定されて GLS の分布が狭い場合、自分のデータでの実効 ICC は文献値より低くなります。文献値をそのまま持ち込むのは、厳密には楽観的です。 自施設で 2 回測定した部分集団があればそこから ICC を算出するのが最善です。なければ、損益分岐 ICC（0.42）が文献値からどれだけ離れているかを示して、多少の下振れでは結論が変わらないと述べます。 ### どの ICC を使うかで話が変わる 同じ「GLS の ICC」でも中身が複数あります。 | 種類 | 含まれるばらつき | 値の傾向 | |---|---|---| | 同一観察者・同一画像の再計測 | トレースの揺れのみ | 高い | | 観察者間・同一画像 | ＋読影者の違い | 中間 | | test-retest（別の日に撮り直し） | ＋撮像条件、生理的日内変動 | 低い | 回帰希釈で問題になるのは**ベースラインの測定値が真の値からどれだけずれているか**なので、本来は test-retest に近いものが正しい参照です。同一画像の再計測 ICC は高く出やすく、これを使うと希釈の影響を過小評価します。引用する文献がどの ICC を報告しているかは確認してください。 ### 線形性と正規性 $\beta_{\text{観測}} = \beta_{\text{真}} \times \mathrm{ICC}$ は、真の値と誤差が独立で関係が線形であることを前提にしています。GLS のように床効果・天井効果がある指標では、極端な値の付近で崩れます。散布図を見て明らかに曲がっていないかは確認してください。 ### 補正はあくまで一次近似 希釈補正は、指定した ICC の値が正しいという前提の上に立ちます。`eivreg` が返す信頼区間は「ICC = 0.85 が正しいとして」の区間で、ICC 自体の不確実性は入っていません。上の感度分析を添えるのが誠実です。 ## Stata 版 同じ内容を Stata で実行します。乱数生成器が R と違うため数字は完全一致しませんが、結論は同じです。 ### 表示桁の設定 Stata の既定では `.8009146` のように桁が多く出て読みにくくなります。do-file の冒頭で表示桁をそろえておきます。 ```stata set cformat %9.4f // 係数、標準誤差、信頼区間 set pformat %5.3f // p 値 set sformat %8.2f // t 値、z 値 ``` `regress` や `eivreg` の係数表に効きます。引数なしで `set cformat` と打つと既定に戻ります。恒久設定にもできますが、do-file の冒頭に書いておけば他の人が実行しても同じ出力になるので、こちらを勧めます。 幅の上限は 9 で、`%12.4f` のような指定は `width too large` というエラーになります。係数が 10000 を超えるデータでは `%9.4f` だと `1.13e+04` と指数表記に化けるので、`%9.3f` に落とします。 `summarize` の表はこの設定の対象外です。変数側の表示形式を決めてから `format` オプションを付けます。`tabstat` はどちらにも属さず、`format(%6.2f)` をコマンドに直接書きます。 ### 第1部: 模試の点数 ```stata * 402_part1.do 模試の点数で平均への回帰を体感する clear all set seed 20260814 set obs 5000 set cformat %9.4f set pformat %5.3f set sformat %8.2f gen double jitsuryoku = rnormal(60, 10) gen double kai1 = jitsuryoku + rnormal(0, 5) gen double kai2_A = jitsuryoku + 10 + rnormal(0, 5) gen double kai2_B = 60 + 0.50*(jitsuryoku - 60) + 10 + rnormal(0, 5) format jitsuryoku kai1 kai2_A kai2_B %9.2f summarize jitsuryoku kai1 kai2_A kai2_B, format gen double nobi_A = kai2_A - kai1 gen double nobi_B = kai2_B - kai1 gen byte band = . replace band = 1 if kai1 < 45 replace band = 2 if kai1 >= 45 & kai1 < 55 replace band = 3 if kai1 >= 55 & kai1 < 65 replace band = 4 if kai1 >= 65 & kai1 < 75 replace band = 5 if kai1 >= 75 label define bl 1 "<45" 2 "45-55" 3 "55-65" 4 "65-75" 5 ">=75" label values band bl tabstat kai1 nobi_A nobi_B, by(band) statistics(mean n) format(%6.2f) regress kai2_A kai1 regress kai2_B kai1 ``` ``` Variable | Obs Mean Std. dev. Min Max -------------+--------------------------------------------------------- jitsuryoku | 5,000 60.05 9.77 21.67 95.64 kai1 | 5,000 60.06 10.94 22.49 100.23 kai2_A | 5,000 70.13 11.02 27.57 110.52 kai2_B | 5,000 70.03 6.89 42.25 93.94 band | kai1 nobi_A nobi_B -------+------------------------------ <45 | 39.87 13.70 22.03 45-55 | 50.68 12.11 15.74 55-65 | 60.04 10.10 10.11 65-75 | 69.21 8.11 4.24 >=75 | 80.36 6.18 -2.35 -------------------------------------- . regress kai2_A kai1 ------------------------------------------------------------------------------ kai2_A | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- kai1 | 0.8009 0.0086 92.66 0.000 0.7840 0.8179 _cons | 22.0232 0.5277 41.74 0.000 20.9888 23.0577 ------------------------------------------------------------------------------ . regress kai2_B kai1 ------------------------------------------------------------------------------ kai2_B | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- kai1 | 0.3892 0.0070 55.53 0.000 0.3755 0.4030 _cons | 46.6584 0.4279 109.03 0.000 45.8194 47.4973 ------------------------------------------------------------------------------ ``` シナリオA が 0.8009、シナリオB が 0.3892 です。R の 0.7867 と 0.3944 に対応しており、理論値 0.80 と 0.40 のまわりに散らばっています。`kai2_B` の SD が 6.89 と `kai1` の 10.94 から縮んでいる点も、R と同じです。 ### 第2部: GLS ```stata * 402_part2.do clear all set seed 20260814 set obs 200000 set cformat %9.4f set pformat %5.3f set sformat %8.2f local mu = 12 // 真のベースライン GLS（絶対値, %）の平均 local sd = 4 // 個人間 SD local icc = 0.85 // GLS 測定の再現性 local delta = 3.0 // TAVI による一律の改善量（%） local btrue = 0.50 // シナリオB の真の回帰係数 local sde = `sd' * sqrt((1 - `icc') / `icc') gen double shin0 = rnormal(`mu', `sd') gen double shin1_A = shin0 + `delta' gen double shin1_B = `mu' + `btrue'*(shin0 - `mu') + `delta' gen double kansoku0 = shin0 + rnormal(0, `sde') gen double kansoku1_A = shin1_A + rnormal(0, `sde') gen double kansoku1_B = shin1_B + rnormal(0, `sde') regress kansoku1_A kansoku0 regress kansoku1_B kansoku0 gen double henka_B = kansoku1_B - kansoku0 regress henka_B kansoku0 * 信頼性を指定して希釈補正する eivreg kansoku1_B kansoku0, r(kansoku0 0.85) ``` ``` ------------------------------------------------------------------------------ kansoku1_A | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- kansoku0 | 0.8484 0.0012 720.93 0.000 0.8461 0.8507 _cons | 4.8165 0.0150 320.74 0.000 4.7871 4.8459 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ kansoku1_B | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- kansoku0 | 0.4231 0.0010 442.96 0.000 0.4213 0.4250 _cons | 9.9183 0.0122 813.66 0.000 9.8944 9.9422 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ henka_B | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- kansoku0 | -0.5769 0.0010 -603.89 0.000 -0.5787 -0.5750 _cons | 9.9183 0.0122 813.66 0.000 9.8944 9.9422 ------------------------------------------------------------------------------ ``` 観測値は 0.8484 と 0.4231 で、理論値 0.850 と 0.425 に一致します。変化量の回帰は −0.5769 で、理論値 $0.50 \times 0.85 - 1 = -0.575$ に一致します。 `eivreg` は信頼性を指定するだけで希釈補正をしてくれます。信頼区間まで出るので、手で割るよりこちらが実用的です。 ``` Errors-in-variables regression Assumed Variable reliability ---------------------------- Number of obs = 200,000 kansoku0 0.8500 F( 1, 199998) = 2.0e+05 * 1.0000 Prob > F = 0.0000 R-squared = 0.5826 Root MSE = 1.68653 ------------------------------------------------------------------------------ kansoku1_B | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- kansoku0 | 0.4978 0.0011 441.85 0.000 0.4956 0.5000 _cons | 9.0223 0.0142 637.31 0.000 8.9945 9.0500 ------------------------------------------------------------------------------ ``` 補正後の 0.4978 は真値 0.50 に一致しています。同じデータの通常 OLS が 0.4231 だったのと比べてください。 ## Reference - [[ICC を直感的に理解する]] — ICC そのものの説明 - [[平均への回帰と回帰希釈を直感的に理解する]] — このノートが検証している理屈 - [[R - ICC を計算してみよう（評価者間一致研究）]] - [[R - ICC を計算してみよう（試薬評価研究）]] - [[生物統計における仮定とは]] - [[R - 回帰係数を理解しよう]] - [[Stata - 簡単な統計量を計算しよう]] - Bland JM, Altman DG. Regression towards the mean. BMJ. 1994;308:1499. - Bland JM, Altman DG. Some examples of regression towards the mean. BMJ. 1994;309:780. - Hutcheon JA, Chiolero A, Hanley JA. Random measurement error and regression dilution bias. BMJ. 2010;340:c2289. - Tu YK, Gilthorpe MS. Revisiting the relation between change and initial value: a review and evaluation. Stat Med. 2007;26:443-57.