LLMとGLSと共分散とBLUP

1 やること

Mrodeの本(Linear models for the prediction of animal breeding values)の輪読をしていて、

3章のBLUPにでてくる V

\[ V = ZGZ^\top + R \]

がBLUPの予測に実際どう効くのかという疑問があったのでシミュレーションで確認してみた。

2 分散共分散の置き方（先に定義）

まず、誤差と（遺伝）効果の分散共分散を次で定義する：

\[ \mathrm{Var}(e) = I\sigma_e^2 = R \]

\[ \mathrm{Var}(a) = A\sigma_a^2 = G \]

ここで \(A\) は（例：血縁行列 / GRM などの）既知の関係行列を表す。

---

3 Linear Mixed Model (LMM) とは

LMM (Linear Mixed Model) は、
線形回帰に「相関構造（ランダム効果）」を明示的に組み込んだモデルである。
モデルは次で定義される：

\[ y = Xb + Za + e \]

\[ a \sim N(0, G), \qquad e \sim N(0, R), \qquad a \perp e \]

• \(b\)：固定効果
• \(a\)：ランダム効果（ここでは遺伝効果など）
• \(G\)：ランダム効果の分散共分散行列（例：\(G=A\sigma_a^2\)）
• \(R\)：残差の分散共分散行列（例：\(R=I\sigma_e^2\)）

---

4 周辺モデルとして見た LMM

ランダム効果 \(a\) を積分消去すると，

\[ y \sim N(Xb,\ V) \]

\[ V = ZGZ^\top + R \]

となる。

4.1 重要な帰結

• LMM は GLS (Generalized Least Squares) と同じになる（\(y\) の共分散を \(V\) として扱う線形モデル）：

\[ y = Xb + \tilde e \]

where

\[ \tilde e \sim N(0,\ V) \]

---

5 固定効果 \(b\) の推定（BLUE）

\[ \hat b = (X^\top V^{-1} X)^{-1} X^\top V^{-1} y \]

---

6 ランダム効果 \(a\) の推定（BLUP）

\[ \hat a = G Z^\top V^{-1}(y - X\hat b) \]

LMM では：

• \(b\)：BLUE (Best Linear Unbiased Estimator)
• \(a\)：BLUP (Best Linear Unbiased Predictor)

が同時に得られる。

---

7 \(V^{-1}\) は BLUP にどう効いているか

\[ \hat a = \underbrace{G Z^\top}_{\text{構造}} \;\underbrace{V^{-1}}_{\text{重み調整}} \;\underbrace{(y - X\hat b)}_{\text{固定効果で説明できない残差}} \]

\(V^{-1}\) の役割は：

• 分散の大きい（信頼性の低い）方向を 弱める
• 冗長な情報を 割り引く

---

8 Simulation：\(V^{-1}\) による BLUP の変化（\(h^2=0.7\)）

このセクションでは，LMM における BLUP

\[ \hat a = GZ^\top V^{-1} r \qquad (r := y - X\hat\beta) \]

が，観測間共分散 \(V\) の仮定によってどのように変化するかを，最小例（\(n=4\)）で確認する。

---

8.1 設定（観測は常に 1）

ここでは 残差（固定効果で説明できない成分）を一定にするため，

\[ r = y - X\hat\beta = (1,1,1,1)^\top \]

と固定する。
この \(r\) は「全個体が同じ方向（平均方向）にずれている」状況を表し，
\(V\) に強い相関がある場合に \(V^{-1}\) がどの程度その方向を抑制するかが分かりやすい。

また，ランダム効果の割当は

\[ Z = I \]

とし，遺伝構造（関係行列）も

\[ A=I \]

とする（すなわち個体間の遺伝相関は無い，独立な遺伝効果）。

---

8.2 遺伝率（heritability）の設定：\(h^2=0.7\)

本例では分散成分を

\[ \mathrm{Var}(a) = G = \sigma_a^2 I,\qquad \mathrm{Var}(e) = R = \sigma_e^2 I \]

とし，遺伝率を

\[ h^2 = \frac{\sigma_a^2}{\sigma_a^2+\sigma_e^2} = 0.7 \]

に固定する。 ベースラインの \(V\) はスカラー倍の単位行列に一致する：

\[ V = ZGZ^\top + R = 0.7I + 0.3I = I \]

---

8.3 比較する 3 つのケース

本コードでは，次の 3 パターンの \(V\) を用意し，それぞれの

\[ \hat a = GZ^\top V^{-1}r \]

を比較する。

8.4 Case A（独立・同分散）

\[ V_1 = ZGZ^\top + R = I \] この場合 \(V_1^{-1} = I\) であり，

\[ \hat a = G r = 0.7\,r \]

となる。

---

8.5 Case B（観測 1 と 4 のみ相関を付加）

Case A の \(V\) をベースに，(1,4) と (4,1) のみ共分散を付加する：

\[ (V_2)_{14} = (V_2)_{41} = \rho\quad (\rho=0.7) \]

ここでは，観測 1 と 4 が「冗長な情報」を共有していると仮定したことになる。
その結果，\(V_2^{-1}\) は 1 と 4 の方向の情報を独立とみなして過大評価せず，
\(\hat a\) の推定値も Case A から変化する。

---

8.6 Case C（ほぼ完全相関：対角 1，非対角 0.95）

\[ (V_3)_{ii}=1,\qquad (V_3)_{ij}=0.95\ (i\neq j) \]

これは「全観測がほぼ同じ誤差（共通成分）を共有している」極端な状況であり，
\(V_3\) は 1 方向（平均方向）に大きな分散を持つ一方で，
それ以外の方向の分散は小さくなる（ほぼランク 1 に近い）。

本例では\(V_3^{-1}\) によってその方向が強く抑制され，\(\hat a\) は Case A よりも大きく縮小されることが期待される。

# -------------------------
# Setup (n = 4)
# -------------------------
n <- 4
Z <- diag(n)
r<-c(1,1,1,1)
sigma_a2 <- 0.7
sigma_e2 <- 0.3  # <- もし Sigma=I にしたいなら 1-sigma_a2

G <- sigma_a2 * diag(n)  # Var(a)
R <- sigma_e2 * diag(n)  # Var(e)

# -------------------------
# Case A: Sigma = I
# -------------------------
Sigma1 <- Z %*% G %*% t(Z) + R 
Sinv1  <- solve(Sigma1)
a_hat1 <- as.vector(G %*% t(Z) %*% Sinv1 %*% r)  

# -------------------------
# Case B: Sigma = I + off-diagonal correlation (1,4)
# -------------------------
rho <- 0.7
Sigma2 <- Z %*% G %*% t(Z) + R 
Sigma2[1,4] <- rho
Sigma2[4,1] <- rho
Sinv2  <- solve(Sigma2)
a_hat2 <- as.vector(G %*% t(Z) %*% Sinv2 %*% r)

# -------------------------
# Case C: diag = 1, offdiag = 0.99 (almost fully correlated)
# -------------------------
rho_all <- 0.95
Sigma3 <- matrix(rho_all, n, n)
diag(Sigma3) <- 1
Sinv3  <- solve(Sigma3)
a_hat3 <- as.vector(G %*% t(Z) %*% Sinv3 %*% r)

# -------------------------
# Data for plotting
# -------------------------
df <- data.frame(
  sample = factor(1:4),
  r      = r,
  a_I    = a_hat1,
  a_corr = a_hat2,
  a_all  = a_hat3
)

# y-limits (3 panels共通) : r と a_hat の範囲を全部含める
ylim_all <- range(c(df$r, df$a_I, df$a_corr, df$a_all))

# -------------------------
# Helper: simple heatmap (base)
# -------------------------
plot_heat <- function(M, title = ""){
  image(
    x = 1:nrow(M), y = 1:ncol(M), z = t(M)[, nrow(M):1],
    axes = FALSE, xlab = "", ylab = "", main = title
  )
  axis(1, at = 1:nrow(M), labels = 1:nrow(M))
  axis(2, at = 1:ncol(M), labels = rev(1:ncol(M)))
  box()
  # values overlay (見やすいので数値も表示)
  for(i in 1:nrow(M)){
    for(j in 1:ncol(M)){
      text(i, nrow(M)-j+1, labels = sprintf("%.2f", M[i,j]), cex = 0.9)
    }
  }
}

# -------------------------
# Plot: 1st row = barplots (a_hat)
#       2nd row = heatmaps (Sigma)
# -------------------------
par(mfrow = c(2, 3), mar = c(4,4,3,1))

# Panel 1: Sigma = I
barplot(
  height = df$a_I,
  names.arg = df$sample,
  col = "grey70",
  ylim = c(0,1),
  main = expression(Sigma == I),
  ylab = expression(hat(a))
)
abline(h = 0)
abline(h = sigma_a2, lty = "dashed")

# Panel 2: Sigma with (1,4) correlation
barplot(
  height = df$a_corr,
  names.arg = df$sample,
  col = "steelblue",
  ylim = c(0,1),
  main = expression(~Sigma~"(1,4 corr)"),
  ylab = expression(hat(a))
)
abline(h = 0)
abline(h = sigma_a2, lty = "dashed")

# Panel 3: Sigma with offdiag 0.99
barplot(
  height = df$a_all,
  names.arg = df$sample,
  col = "tomato",
  ylim = c(0,1),
  main = expression(Sigma~"diag=1, offdiag=0.99"),
  ylab = expression(hat(a))
)
abline(h = 0)
abline(h = sigma_a2, lty = "dashed")

# -------------------------
# 共通カラースケールを計算
# -------------------------
zlim_all <- range(c(Sigma1, Sigma2, Sigma3))

# -------------------------
# 共通カラースケール
# -------------------------
zlim_all <- range(c(Sigma1, Sigma2, Sigma3))

# -------------------------
# Helper: heatmap (white -> black)
# -------------------------
plot_heat <- function(M, title = "", zlim){
  image(
    x = 1:nrow(M),
    y = 1:ncol(M),
    z = t(M)[, nrow(M):1],
    axes = FALSE,
    xlab = "",
    ylab = "",
    main = title,
    zlim = zlim,
    col = gray.colors(100, start = 1, end = 0)  # white -> black
  )
  axis(1, at = 1:nrow(M), labels = 1:nrow(M))
  axis(2, at = 1:ncol(M), labels = rev(1:ncol(M)))
  box()

  # 数値を重ねる
  for(i in 1:nrow(M)){
    for(j in 1:ncol(M)){
      text(i, nrow(M)-j+1,
           labels = sprintf("%.2f", M[i,j]),
           cex = 0.9,
           col = ifelse(M[i,j] > mean(zlim), "white", "black"))
    }
  }
}

# -------------------------
# Heatmaps (same scale, white -> black)
# -------------------------
plot_heat(Sigma1, title = expression("Heatmap of V: "~Sigma==I), zlim = zlim_all)
plot_heat(Sigma2, title = expression("Heatmap of V: "~Sigma~"(1,4 corr)"), zlim = zlim_all)
plot_heat(Sigma3, title = expression("Heatmap of V: "~Sigma~"(offdiag=0.99)"), zlim = zlim_all)

• 上段（棒グラフ）：各ケースの \(\hat a\) の成分（4サンプル）を表示する
  • 破線：\(\sigma_a^2=0.7\)（Case A で理論的に出る水準）
    
• 下段（ヒートマップ）：各ケースの \(V\) を白→黒で表示する（黒いほど値が大きい）

---

8.7 まとめ（このシミュレーションが示すこと）

この 4 サンプル例では，\(r\) を固定した上で，

• \(V\) が 独立（Case A） なら \(V^{-1}\) は単純で，\(\hat a\) は基準値に揃う
• \(V\) が 部分的に相関（Case B） なら，\(V^{-1}\) は相関した観測を冗長として扱い，\(\hat a\) を調整する
• \(V\) が ほぼ完全相関（Case C） なら，平均方向の情報が強く割り引かれ，\(\hat a\) は大きく縮小される

つまり，\(V^{-1}\) は単なる逆行列ではなく，残差から「独立で信頼できる成分」を抽出するフィルタ（白色化・重み付け）として働き，BLUP の縮小量を決める。

9 まとめ

• LMM = 正規分布における GLS
• \(V^{-1}\) は：
  • 信頼できる情報のみを抽出し
  • ランダム効果に割り当てる
• 相関が強いほど，BLUP は控えめになる

10 参考

Mrode, R. (2014). Linear models for the prediction of animal breeding values. https://www.cabidigitallibrary.org/doi/book/10.1079/9780851990002.0000