ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rで線形混合モデル3】なぜ個人差を考慮する必要があるのか

生体信号解析、心理実験、行動データ、教育データなど、同じ実験対象者から複数回の測定を行った結果のデータは多くの分野で扱われています。そのようなデータに対して「実験対象者間の個人差はとりあえず無視して、全部まとめて回帰してみよう」と考えるかもしれません。
 しかし、このような安易な簡単化は、統計的には誤った危険な解釈を生むことがあります。そのような過ちを避けるために、線形混合モデル(Linear Mixed Model: LMM)と呼ばれる方法が役に立つかもしれません。線形混合モデルとは、説明変数の効果として全体に共通する成分(固定効果)と、実験対象者ごとに異なる成分(ランダム効果)を同時にモデル化することで、個人差や群内構造を明示的に扱う統計モデルです。
 今回の記事では、Rを使ってサンプルデータを生成し、個人差を無視するとどのような過ちを犯す可能性があるのかを考えてみたいと思います。

1. 個人差を無視すると「解釈を間違える」サンプルデータ

 ここでは、3種類のサンプルデータを生成するRスクリプトを提示します。

1-1. 例1:個人差を無視すると無相関に見えるデータ

 下図のデータは、同じ Subject(実験対象者)から複数回測定されたデータを模擬しています。 左の図では Subject ごとに色分けされており、右の図では Subject の区別を意図的に無視してすべての点をまとめています。下段の図では、個人ごとに回帰直線を当てはめた場合(左)と個人差を無視して回帰直線を当てはめた場合(右)が描いてあります。

個人差を無視すると無相関に見えるデータの例

このデータの重要な特徴は次の 3 点です。

  1. Subject ごとにデータの分布中心(切片)が異なる
     左図を見ると、各 Subject の点群はおおむね直線的な傾向を示していますが、その直線の位置(上下方向)は Subject ごとにずれています。これは「同じ説明変数 x に対しても、Subject によって応答変数 y の平均レベルが異なる」ことを意味します。

  2. 各 Subject 内では正の関係が存在する
     下段左図のように Subject ごとに回帰直線を引くと、どの Subject においても 「x が大きくなると y も増える」 という関係が明確に見て取れます。

  3. Subject 間の違いが、データ全体の構造を支配している
     Subject を区別しない右図では、データは一見するとばらついた雲のように見え、個々の Subject 内で見えていた構造が視覚的に消えてしまいます。

このデータに個人差を無視して単純な線形回帰を適用すると、各 Subject 内では明確に正の傾きが存在するにもかかわらず、傾きがほぼ 0 の水平な直線が推定されています。つまり、「説明変数 x は応答変数 y にほとんど影響しない」という解釈に陥ってしまいます。しかし、これは事実ではありません。各 Subject の内部では、明確に影響が存在しています。
 この例では、線形回帰の代わりに、ランダム切片を仮定した線形混合モデルを適用すると、この誤った解釈を避けることができます。線形混合モデルを実行するRスクリプトは、次節で示します。
 みなさんが、ここに示したデータを生成するために、以下のRスクリプトを実行してください。

############################################################
# Linear Mixed Model (LMM) 学習用データ生成スクリプト(その1)
############################################################
set.seed(7)
# =========================
# モデルの真のパラメタ
# =========================
beta_slope <- 1
beta_intercept <- c(2.4, 1.2, 0, -1.2, -2.4)
x_center <- c(1, 2, 3, 4, 5)

# =========================
# データサイズ
# =========================
n_per_subject <- 30
n_subject <- length(beta_intercept)

# =========================
# データ生成
# =========================
Subject <- integer()
x <- numeric()
y <- numeric()
for (i in 1:n_subject) {
  Subject <- c(Subject, rep(i, n_per_subject))
  x_i <- rnorm(n_per_subject, mean = x_center[i], sd = 0.7)
  y_i <- beta_slope * x_i +
    beta_intercept[i] +
    rnorm(n_per_subject, mean = 0, sd = 0.4)
  x <- c(x, x_i)
  y <- c(y, y_i)
}
dat <- data.frame(
  Subject = factor(Subject),
  x = x,
  y = y
)

# =========================
# 可視化設定
# =========================
cols <- c("blue", "orange", "darkgreen", "red", "purple")
# 点用(半透明)
cols_pt <- adjustcolor(cols, alpha.f = 0.6)
# 2枚横並び
op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3.2, 1.0), las = 1)
# ============================================================
# 左:Subject別
# ============================================================
plot(
  dat$x, dat$y,
  col = cols_pt[dat$Subject],
  pch = 16,
  cex = 1.05,
  xlab = "説明変数 x",
  ylab = "応答変数 y",
  main = "Subject 別"
)
legend(
  "bottomright",
  legend = paste("Subject", 1:n_subject),
  col = cols_pt,
  pch = 16,
  bty = "n",
  cex = 1.1
)

# ============================================================
# 右:Subject を区別しない場合
# ============================================================
plot(
  dat$x, dat$y,
  col = adjustcolor("black", alpha.f = 0.35),
  pch = 16,
  cex = 1.0,
  xlab = "説明変数 x",
  ylab = "応答変数 y",
  main = "Subject の区別なし"
)
###
par(op)  # 描画設定を戻す

1-2. 例2:傾きの個人差とサンプル数の偏りにより、全体回帰が誤解を招くデータ

 下図のデータは、同じ Subject(実験対象者)から複数回測定されたデータを模擬しています。 左の図では Subject ごとに色分けされており、右の図では Subject の区別を意図的に無視して、すべての点をまとめています。下段の図では、各 Subject ごとに回帰直線を当てはめた場合(左)と、Subject を区別せずに全体で線形回帰を当てはめた場合(右)が描かれています。
 この例では、切片(平均レベル)はすべての Subject で同じに設定されていますが、傾き(x に対する反応の強さ)とサンプル数が Subject ごとに異なる点が重要です。

傾きの個人差とサンプル数の偏りにより、全体回帰が誤解を招く例

このデータの重要な特徴は、次の 3 点です。

  1. Subject ごとに傾きが大きく異なる(反応の強さに個人差がある)
     下段左図を見ると、各 Subject における回帰直線の傾きが大きく異なっていることが分かります。ある Subject では x に対する y の変化がほとんど見られない一方で、別の Subject では強い正の関係が現れています。これは、説明変数 x に対する反応の仕方そのものが個人ごとに異なることを意味しており、ランダム傾きを考慮すべき典型的な状況です。

  2. 切片は同じでも、傾き差だけでデータ構造は大きく変わる
     上段左図では、各 Subject の点群は同程度の y の水準から始まっていますが、x が大きくなるにつれて Subject ごとの差が拡大しています。つまりこのデータでは、平均レベル(切片)の個人差は存在しないにもかかわらず、傾きの個人差だけでデータの階層構造が生じていることが分かります。

  3. サンプル数の偏りにより、全体回帰が一部の Subject に強く引っ張られる
     凡例に示されている通り、Subject 1 はサンプル数が多く、他の Subject は比較的少なくなっています。そのため、Subject を区別しない右図および下段右図の全体回帰では、データ数の多い Subject の傾きが過大に反映されやすくなります。結果として、全体の回帰直線は一見もっともらしく見えるものの、どの Subject の挙動も正確には代表していない直線になっています。

このデータに対して Subject を無視して単純な線形回帰を適用すると、下段右図のように、全体としては緩やかな正の傾きが推定されます。しかし、この回帰直線を「個人内で x を変化させたときの平均的な反応」と解釈するのは危険です。実際には、下段左図が示すように Subject ごとに傾きが大きく異なるため、傾きの個人差を考慮する必要があります。
 このような場合には、切片は共通としつつ、傾きの個人差をランダム効果として扱う線形混合モデルを適用するのが自然です。とはいえ、実際のデータ分析では、推定が安定するように、切片と傾きの両者にランダム効果を仮定した線形混合モデルを使った方が良いこともあります。そのような線形混合モデルの当てはめを実行する R スクリプトは、次節で示します。
 みなさんが、ここに示したデータを生成するために、以下の R スクリプトを実行してください。

############################################################
# Linear Mixed Model (LMM) 学習用データ生成スクリプト(その2)
############################################################
set.seed(7)
# ============================================================
# 1) 真のパラメタ(ユーザ設定)
# ============================================================
beta_slope     <- c(0.1, 0.2, 0.75, 1.0, 1.3)  # Subject ごとの傾き
beta_intercept <- -0.5                         # 共通切片

n_sub    <- c(100, 50, 20, 20, 20)             # Subject ごとのサンプル数
x_center <- rep(2.5, length(beta_slope))       # Subject ごとの x 平均

n_subject <- length(beta_slope)

# ============================================================
# 2) データ生成
# ============================================================
dat_list <- vector("list", n_subject)

for (i in seq_len(n_subject)) {
  x_i <- rnorm(n_sub[i], mean = x_center[i], sd = 0.8)
  y_i <- beta_slope[i] * x_i +
    beta_intercept +
    rnorm(n_sub[i], mean = 0, sd = 0.4)

  dat_list[[i]] <- data.frame(
    Subject = factor(i, levels = seq_len(n_subject)),
    x = x_i,
    y = y_i
  )
}

dat <- do.call(rbind, dat_list)

# ============================================================
# 3) 可視化
# ============================================================
cols <- c("blue", "orange", "darkgreen", "red", "purple")
cols_pt <- adjustcolor(cols, alpha.f = 0.6)

op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3.2, 1.0), las = 1)

# --- 左:Subject 別 ---
plot(
  dat$x, dat$y,
  col = cols_pt[as.integer(dat$Subject)],
  pch = 16, cex = 1.05,
  xlab = "説明変数 x",
  ylab = "応答変数 y",
  main = "Subject 別"
)

legend_labels <- sprintf("Subject %d (N=%d)", seq_len(n_subject), n_sub)
legend(
  "topleft",
  legend = legend_labels,
  col = cols_pt,
  pch = 16,
  bty = "n",
  cex = 1.05
)

# --- 右:Subject を区別しない ---
plot(
  dat$x, dat$y,
  col = adjustcolor("black", alpha.f = 0.35),
  pch = 16, cex = 1.0,
  xlab = "説明変数 x",
  ylab = "応答変数 y",
  main = "Subject の区別なし"
)

par(op)

1-3. 例3:個人差を無視すると強い相関に見えるデータ

 下図のデータは、同じ Subject(実験対象者)から複数回測定されたデータを模擬しています。左の図では Subject ごとに色分けされており、右の図では Subject の区別を意図的に無視して、すべての点をまとめています。下段の図では、Subject 別の散布図(左)と、Subject を区別せずに線形回帰を当てはめた結果(右)が描かれています。

個人差を無視すると「強い相関」に見えるデータの例

このデータの重要な特徴は、次の 3 点です。

  1. Subject ごとに (x, y) の平均(中心)がそろって移動している
     左図を見ると、Subject 1 から Subject 5 に向かうにつれて、点群全体が右上方向へ移動しています。つまり、Subject 間で 説明変数 x の平均も応答変数 y の平均も同時に増加しており、Subject の違いそのものがデータ全体の配置を決めています。

  2. 各 Subject 内で観測される変動幅は相対的に狭い
     各 Subject の点群はそれぞれまとまりを持っており、Subject 内では x の取りうる範囲が比較的狭くなっています。このとき、データ全体の「右上がり」の見かけは、Subject 内の関係というよりも、Subject 間の平均値の違いがメインの構造として現れやすくなります。

  3. Subject を混ぜると、Subject 間の差が「x と y の関係」に見えてしまう
     Subject を区別しない右図では、点群全体がはっきりと右上がりに見えます。さらに下段右図では、単純な線形回帰の回帰直線も強い正の傾きを示しています。しかし、この回帰直線が主に表しているのは、個人内(Subject 内)での (x) の効果ではなく、Subject 間で平均が同時に変化している(個人差の構造)である可能性があります。

このデータに対して Subject を無視して単純な線形回帰を適用すると、「x が増えると y が増える」という強い関係が推定されます。ところが、この結論をそのまま「個人内で x を変化させれば y が変化する」と解釈してしまうと、誤りになることがあります。言い換えると、単純な線形回帰は Subject 間の違い(個人差)Subject 内の関係(個人内効果)を区別できないため、「何の効果を見ているのか」が曖昧になりやすいのです。
 このような場合には、線形混合モデル(LMM)を用いて Subject ごとの差をモデル化したり、さらに踏み込んで Subject 内変動と Subject 間変動を分離することで、「個人内の効果」をより適切に評価できるようになります。そのような線形混合モデルの当てはめを実行する R スクリプトは、次節で示します。
 みなさんが、ここに示したデータを生成するために、以下のRスクリプトを実行してください。

############################################################
# Linear Mixed Model (LMM) 学習用データ生成スクリプト(その3)
############################################################
set.seed(7)
# ============================================================
# 1) パラメタ(ユーザ設定)
# ============================================================
x_center <- c(1, 2, 3, 4, 5)  # Subject ごとの x 平均
y_center <- c(1, 2, 3, 4, 5)  # Subject ごとの y 平均

n_per_subject <- 50
sd_within     <- 0.6

n_subject <- length(x_center)

# ============================================================
# 2) データ生成
# ============================================================
dat_list <- vector("list", n_subject)
for (i in seq_len(n_subject)) {
  x_i <- rnorm(n_per_subject, mean = x_center[i], sd = sd_within)
  y_i <- rnorm(n_per_subject, mean = y_center[i], sd = sd_within)
  dat_list[[i]] <- data.frame(
    Subject = factor(i, levels = seq_len(n_subject)),
    x = x_i,
    y = y_i
  )
}
dat <- do.call(rbind, dat_list)

# ============================================================
# 3) 可視化
# ============================================================
cols    <- c("blue", "orange", "darkgreen", "red", "purple")
cols_pt <- adjustcolor(cols, alpha.f = 0.6)
op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3.2, 1.0), las = 1)

# --- 左:Subject 別 ---
plot(
  dat$x, dat$y,
  col = cols_pt[as.integer(dat$Subject)],
  pch = 16, cex = 1.05,
  xlab = "説明変数 x",
  ylab = "応答変数 y",
  main = "Subject 別"
)
legend(
  "topleft",
  legend = sprintf("Subject %d (N=%d)", seq_len(n_subject), n_per_subject),
  col = cols_pt,
  pch = 16,
  bty = "n",
  cex = 1.05
)

# --- 右:Subject を区別しない ---
plot(
  dat$x, dat$y,
  col = adjustcolor("black", alpha.f = 0.35),
  pch = 16, cex = 1.0,
  xlab = "説明変数 x",
  ylab = "応答変数 y",
  main = "Subject の区別なし"
)

par(op)

2. 線形混合モデル当てはめの実際

 ここでは、データに対して回帰直線を当てはめる方法として、個人差を無視した単純な線形回帰モデルに加えて、次の 3 種類の線形混合モデルを取り上げます。

  1. ランダム切片モデル 切片のみに個人差があると仮定するモデル

  2. ランダム傾きモデル 傾きのみに個人差があると仮定するモデル

  3. ランダム切片・ランダム傾きモデル 切片と傾きの両方に個人差があると仮定するモデル

それぞれについて、R を用いた当てはめの例を示します。 これらのスクリプトを試すためのサンプルデータには、前節で示したデータを用いてください。

2-0. 個人差を無視した単純な線形回帰モデル

 まずは線形混合モデル(LMM)ではなく、Subject(実験対象者)を区別せず、すべての観測を独立同分布とみなして当てはめる 単純な線形回帰モデル(LM)を実行します。以下に、そのための R スクリプト例を示します。

############################################################
# Linear Model (LM) 解析スクリプト
#
# 注意:
#   この解析では「Subject を区別しない」。
############################################################

# =========================
# 1) モデル当てはめ(Subject を区別しない LM)
# =========================

fit_lm <- lm(y ~ x, data = dat)

# =========================
# 2) 出力を文字列として回収
# =========================

out_text <- c(
  "======================================================",
  " Linear Model (LM) 解析結果まとめ",
  " (Subject を区別しない単純回帰)",
  "======================================================",
  "",
  "【解析の前提】",
  "・Subject(被験者)を区別しない",
  "・すべての観測点を独立同分布と仮定",
  "・被験者ごとの切片差・傾き差は誤差に含めて扱う",
  "",
  "【モデル式】",
  capture.output({
    cat("LM : ", deparse(formula(fit_lm)), "\n")
  }),
  "",
  "【推定された回帰係数】",
  capture.output(print(coef(fit_lm))),
  "",
  "【詳細 summary】",
  capture.output(summary(fit_lm))
)
# 3) 結果の表示
cat(paste(out_text, collapse = "\n"))
■ 実行例

 例1のデータに対して、Subject(実験対象者)を区別せず、すべての観測点をまとめて単純な線形回帰モデル(LM)を当てはめると、下に示したような結果が得られます。

【モデル式】
LM :  y ~ x 

【推定された回帰係数】
(Intercept)           x 
3.086427190 0.004016736 

 まず、推定された回帰式は

\displaystyle{
y \approx 3.09 + 0.004,x
}

となっており、傾きはほぼ 0、すなわち x が増加しても y はほとんど変化しない、という結果になっています。この結果だけを見ると「説明変数 x は応答変数 y にほとんど影響しない」と結論づけてしまいそうになります。

2-1. ランダム切片を仮定した線形混合モデル

 以下は、説明変数 x の効果(傾き)は全 Subject で共通としつつ、応答変数 y の平均レベル(切片)に Subject ごとの個人差がある(ランダム切片)と仮定して、線形混合モデルを当てはめる R スクリプト例です。

############################################################
# Linear Mixed Model (LMM) 解析スクリプト(ランダム切片のみ)
#
# 解析の前提:
#   - Subject(被験者)を区別して扱う
#   - 被験者ごとの「切片のずれ」をランダム効果としてモデル化する(ランダム切片)
#   - 傾き x の効果は全被験者で共通(固定効果)
############################################################

# =========================
# 0) パッケージ読み込み
# =========================
pkg_needed <- c("lme4", "lmerTest")
for (p in pkg_needed) {
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
}
library(lme4)
library(lmerTest)

# =========================
# 1) モデル当てはめ(ランダム切片のみ)
# =========================
# モデル: y = (固定切片 + 固定傾き*x) + (Subjectごとの切片偏差) + 誤差
fit_lmm_ri <- lmer(y ~ x + (1 | Subject), data = dat, REML = TRUE)

# =========================
# 2) ICC(被験者内相関の目安)を計算
# =========================
vc <- as.data.frame(VarCorr(fit_lmm_ri))
var_subj <- vc[vc$grp == "Subject", "vcov"]
var_res  <- vc[vc$grp == "Residual", "vcov"]
icc <- var_subj / (var_subj + var_res)

# =========================
# 3) 出力を文字列として回収
# =========================
out_text <- c(
  "======================================================",
  " Linear Mixed Model (LMM) 解析結果まとめ",
  " (ランダム切片モデル:Subject を区別して切片差を表現)",
  "======================================================",
  "",
  "【解析の前提】",
  "・Subject(被験者)を区別して扱う",
  "・被験者ごとの平均レベルの差(切片差)をランダム効果で表現する",
  "・x の効果(傾き)は全被験者で共通(固定効果)",
  "",
  "【モデル式】",
  capture.output({
    cat("LMM(ランダム切片) : ", deparse(formula(fit_lmm_ri)), "\n")
  }),
  "",
  "【固定効果(Fixed effects)】",
  capture.output(print(fixef(fit_lmm_ri))),
  "",
  "【ランダム効果(Subject 別の切片偏差)】",
  capture.output(print(ranef(fit_lmm_ri)$Subject)),
  "",
  "【分散成分(Variance components)】",
  capture.output(print(VarCorr(fit_lmm_ri))),
  "",
  "【ICC(被験者内相関の目安)】",
  capture.output(print(icc)),
  "",
  "【詳細 summary】",
  capture.output(summary(fit_lmm_ri)),
  "\n"
)

# 4) 結果を表示
cat(paste(out_text, collapse = "\n"))
■ 実行例(ランダム切片モデル)

 次に、例1のデータに対して、Subject(実験対象者)ごとに切片が異なることを許したランダム切片モデル

\displaystyle{
y \sim x + (1 \mid \mathrm{Subject})
}

を当てはめた結果を見てみましょう。解析結果は以下の通りです。

【モデル式】
LMM(ランダム切片) :  y ~ x + (1 | Subject) 

【固定効果(Fixed effects)】
(Intercept)           x 
  0.1707446   0.9491374 

まず注目すべき点は、固定効果として推定された傾きが約 0.95 と、単純な線形回帰(LM)とは大きく異なる値になっていることです。 この結果は、「同一 Subject の中で見ると、説明変数 x が 1 増えると、応答変数 y も平均して約 0.95 増加する」という関係が、データ全体として明確に存在していることを示しています。すなわち、例1のデータに本来含まれていた 個人内の正の関係 が、適切に回復されています。
 一方、固定切片の推定値は約 0.17 となっていますが、これは「平均的な Subject」における基準レベルを表すものであり、個々の Subject の切片は次に示すランダム効果によって補正されます。

【ランダム効果(Subject 別の切片偏差)】
  (Intercept)
1  2.30562580
2  1.20555235
3  0.01579101
4 -1.24519837
5 -2.28177078

ここに示されている値は、各 Subject が全体平均の切片からどれだけ上下にずれているかを表しています。たとえば Subject 1 は平均よりも高い水準にあり、Subject 5 は低い水準にあることが分かります。重要なのは、これらの 切片差を「誤差」として捨てるのではなく、モデルの構成要素として明示的に扱っている点です。
 次に、分散成分を見てみましょう。

【分散成分(Variance components)】
 Groups   Name        Std.Dev.
 Subject  (Intercept) 1.84119 
 Residual             0.37109 

ここでは、Subject 間の切片のばらつき(標準偏差)が約 1.84 と、残差誤差(約 0.37)に比べて非常に大きいことが分かります。 このことは、例1のデータにおいて 個人差が主要な変動要因であることを定量的に示しています。
 この点をさらに要約した指標が、以下の ICC(被験者内相関)です。

【ICC(被験者内相関の目安)】
[1] 0.9609647

これは直感的には、「全体のばらつきのうち、どれくらいが個人差によるものか」を示す指標です。ICC が無視できない大きさであれば、Subject を無視した解析は危険だと分かります。ICC が約 0.96 という非常に高い値を取っていることから、データ全体の分散のほとんどが Subject 間の違い(個人差)によって説明されていることが分かります。このような状況で Subject を無視した単純な線形回帰を行うと、先に見たように、個人内の関係が完全に見えなくなってしまうのは自然な結果と言えます。

2-2. ランダム傾きを仮定した線形混合モデル

 以下は、切片は全 Subject で共通としつつ、説明変数 x に対する反応の強さ(傾き)に Subject ごとの個人差がある(ランダム傾き)と仮定して、線形混合モデルを当てはめる R スクリプト例です。

############################################################
# Linear Mixed Model (LMM) 解析スクリプト(ランダム傾きのみ)
#
# 解析の前提:
#   - Subject(被験者)を区別して扱う
#   - 被験者ごとに「傾きのずれ」のみを許す(ランダム傾き)
#   - 切片は全被験者で共通(固定効果のみ)
############################################################

# =========================
# 0) パッケージ読み込み
# =========================
pkg_needed <- c("lme4", "lmerTest")
for (p in pkg_needed) {
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
}
library(lme4)
library(lmerTest)

# =========================
# 1) モデル当てはめ(ランダム傾きのみ)
# =========================
# モデル: y = (固定切片 + 固定傾き*x)
#            + (Subjectごとの傾き偏差*x)
#            + 誤差
#
# ※ (0 + x | Subject) は「切片のランダム効果なし、傾きのみランダム」
fit_lmm_rs <- lmer(y ~ x + (0 + x | Subject), data = dat, REML = TRUE)

# =========================
# 2) 出力を文字列として回収
# =========================
out_text <- c(
  "======================================================",
  " Linear Mixed Model (LMM) 解析結果まとめ",
  " (ランダム傾きモデル:Subject ごとに傾きが異なる)",
  "======================================================",
  "",
  "【解析の前提】",
  "・Subject(被験者)を区別して扱う",
  "・切片は全被験者で共通(固定効果)",
  "・x に対する反応の個人差(傾き差)のみをランダム効果で表現する",
  "・その結果、各 Subject は固有の傾き(反応の強さ)をもつ",
  "",
  "【モデル式】",
  capture.output({
    cat("LMM(ランダム傾き) : ", deparse(formula(fit_lmm_rs)), "\n")
  }),
  "",
  "【固定効果(Fixed effects)】",
  capture.output(print(fixef(fit_lmm_rs))),
  "",
  "【ランダム効果(Subject 別の傾き偏差)】",
  capture.output(print(ranef(fit_lmm_rs)$Subject)),
  "",
  "【分散成分(Variance components)】",
  capture.output(print(VarCorr(fit_lmm_rs))),
  "",
  "(注)このモデルでは切片のランダム効果を入れていないため、",
  "     VarCorr には主に「傾き分散」と「残差分散」が出力されます。",
  "",
  "【詳細 summary】",
  capture.output(summary(fit_lmm_rs)),
  "\n"
)

# 3) 結果を表示
cat(paste(out_text, collapse = "\n"))
■ 実行例(ランダム傾きモデル)

 次に、例2のデータに対して、切片は全 Subject で共通とし、説明変数 x に対する反応の強さ(傾き)のみが Subject ごとに異なると仮定したランダム傾きモデル

\displaystyle{
y \sim x + (0 + x \mid \mathrm{Subject})
}

を当てはめた結果を解説します。解析結果は以下の通りです。

【固定効果(Fixed effects)】
(Intercept)           x 
 -0.4625385   0.6623582 

まず、固定効果として推定された傾きは約 0.66 となっています。 これは、「全 Subject を平均した意味で見ると、x が 1 増えると y は平均して約 0.66 増加する」という平均的な関係を表しています。ただし、この値はすべての Subject に共通する反応の強さを意味するものではありません。
 固定切片の推定値は約 −0.46 ですが、このモデルでは切片にランダム効果を入れていないため、これはあくまで「共通の基準点」を表す補助的な量です。本モデルの本質は、次に示す 傾きの個人差にあります。

【ランダム効果(Subject 別の傾き偏差)】
            x
1 -0.55630319
2 -0.47405934
3  0.07173252
4  0.32053223
5  0.63809777

ここに示されている値は、各 Subject の傾きが、全体平均の傾き(固定効果)からどれだけずれているかを表しています。 たとえば、

  • Subject 1 では \displaystyle{0.66 - 0.56 \approx 0.10} と、x の効果は非常に弱く、
  • Subject 5 では \displaystyle{0.66 + 0.64 \approx 1.30} と、非常に強い正の関係が現れています。

この結果は、例2のデータにおいて「x と y の関係の強さが Subject によって大きく異なる」 ことを明確に示しています。
 次に、分散成分を確認します。

【分散成分(Variance components)】
 Groups   Name Std.Dev.
 Subject  x    0.51306 
 Residual      0.40771 

ここでは、Subject ごとの傾きのばらつき(標準偏差)が約 0.51 と、残差誤差(約 0.41)よりも大きい値を取っています。 これは、例2のデータにおいて 傾きの個人差が主要な変動要因であり、無視できない構造を持っていることを意味します。
 以下では、最後の段落を、指定された単純な線形回帰(LM)の結果を踏まえて修正した形を示します。 それ以前の説明との連続性を保ちつつ、「なぜ LM では誤解が生じるか」が明確になるように解釈を加えています。
 結果の違いを見るため、個人差を無視した単純な線形回帰モデル

\displaystyle{
y \sim x
}

も適用してみます。その結果は以下の通りです。

【推定された回帰係数】
(Intercept)           x 
 -0.2685711   0.2929902 

線形混合モデルの固定効果として推定された傾きは約 0.66でしたが、このモデルでは、傾きは約 0.29 と小さく推定されています。この値は 各 Subject 内での反応の強さを反映したものではなく、Subject 間の違いをすべて平均化した結果にすぎません。
 実際には、ランダム傾きモデルで見たように、x に対する反応の強さは Subject ごとに大きく異なっており、強い正の関係を示す Subject と、ほとんど関係を示さない Subject が混在しています。 単純な線形回帰では、これら異なる傾きを持つデータを一つの直線で近似するため、結果として どの Subject にとっても代表的とは言えない「中途半端な傾き」が推定されます。
 したがって、この例2のデータにおいて LM の結果をそのまま 「x と y には弱い正の相関がある」 と解釈してしまうと、実際には存在する傾きの個人差という本質的構造を見落とすことになります。この点こそが、線形混合モデル、とくにランダム傾きモデルを用いる意義であり、「個人差を明示的にモデル化すること」の重要性を端的に示しています。

2-3. ランダム切片・ランダム傾きを仮定した線形混合モデル

 以下は、応答変数 y の平均レベル(切片)と、説明変数 x に対する反応の強さ(傾き)の両方に Subject ごとの個人差がある(ランダム切片・ランダム傾き)と仮定して、線形混合モデルを当てはめる R スクリプト例です。

############################################################
# Linear Mixed Model (LMM) 解析スクリプト(ランダム切片+ランダム傾き)
#
# 解析の前提:
#   - Subject(被験者)を区別して扱う
#   - 被験者ごとに「切片のずれ」と「傾きのずれ」を許す
#     (ランダム切片+ランダム傾き)
############################################################

# =========================
# 0) パッケージ読み込み
# =========================
pkg_needed <- c("lme4", "lmerTest")
for (p in pkg_needed) {
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
}
library(lme4)
library(lmerTest)

# =========================
# 1) モデル当てはめ(ランダム切片+ランダム傾き)
# =========================
# モデル: y = (固定切片 + 固定傾き*x)
#            + (Subjectごとの切片偏差 + Subjectごとの傾き偏差*x)
#            + 誤差
fit_lmm_rs <- lmer(y ~ x + (1 + x | Subject), data = dat, REML = TRUE)

# =========================
# 2) 出力を文字列として回収
# =========================
out_text <- c(
  "======================================================",
  " Linear Mixed Model (LMM) 解析結果まとめ",
  " (ランダム切片+ランダム傾きモデル:Subject ごとに直線が異なる)",
  "======================================================",
  "",
  "【解析の前提】",
  "・Subject(被験者)を区別して扱う",
  "・被験者ごとの平均レベルの差(切片差)をランダム効果で表現する",
  "・さらに、x に対する反応の個人差(傾き差)もランダム効果で表現する",
  "・その結果、各 Subject は固有の回帰直線(切片・傾き)をもつ",
  "",
  "【モデル式】",
  capture.output({
    cat("LMM(ランダム切片+傾き) : ", deparse(formula(fit_lmm_rs)), "\n")
  }),
  "",
  "【固定効果(Fixed effects)】",
  capture.output(print(fixef(fit_lmm_rs))),
  "",
  "【ランダム効果(Subject 別の切片偏差・傾き偏差)】",
  capture.output(print(ranef(fit_lmm_rs)$Subject)),
  "",
  "【分散成分(Variance components)】",
  capture.output(print(VarCorr(fit_lmm_rs))),
  "",
  "(注)上の VarCorr には、切片分散・傾き分散・切片と傾きの相関が含まれます。",
  "",
  "【詳細 summary】",
  capture.output(summary(fit_lmm_rs))
)

# 3) 結果を表示
cat(paste(out_text, collapse = "\n"))
■ 実行例(ランダム切片・ランダム傾きモデル)

 次に、例3のデータに対して、切片と傾きの両方に Subject ごとの個人差があると仮定した、ランダム切片・ランダム傾きモデル

\displaystyle{
y \sim x + (1 + x \mid \mathrm{Subject})
}

を当てはめた結果を解説します。 重要な前提として、このデータでは 各 Subject の内部では x と y は無相関になるように生成されています。すなわち、個人内の真の傾きは 0 です。
 まず、固定効果の推定結果を見てみましょう。

【固定効果(Fixed effects)】
(Intercept)           x 
  2.6647633   0.0989607 

固定効果として推定された傾きは約 0.10 と非常に小さく、 「全体平均として見ても、x の効果はほぼ存在しない」 という結果になっています。 これは、各 Subject 内で x と y が無相関であるというデータ生成の前提と整合的です。
 次に、Subject ごとのランダム効果を確認します。

【ランダム効果(Subject 別の切片偏差・傾き偏差)】
  (Intercept)            x
1  -1.6884917 -0.052761085
2  -0.7724311 -0.024136516
3   0.0331044  0.001034428
4   0.8025398  0.025077335
5   1.6252786  0.050785838

一見すると、Subject ごとに傾きの違いがあるようにも見えますが、ここで重要なのは、これらの傾き偏差は「真の効果」ではなく、有限サンプルとノイズによって生じた推定上のばらつきである点です。 真のモデルでは個人内の傾きは 0 であり、ランダム傾きは本来不要です。
 このことは、分散成分を見るとより明確になります。

【分散成分(Variance components)】
 Groups   Name        Std.Dev. Corr 
 Subject  (Intercept) 1.303286      
          x           0.040724 1.000
 Residual             0.619137      

ここで、

  • 切片の標準偏差は約 1.30 と大きい一方で
  • 傾きの標準偏差は約 0.04 と極めて小さい

ことが分かります。 これは、データの主要な構造は切片方向(平均レベルの違い)にあり、傾き方向のばらつきはほとんど存在しないことを示しています。
 また、切片と傾きの相関が 1.000 と推定されている点も重要です。 このような極端な相関は、

  • 傾きの分散がほぼ 0 に近い
  • モデルがデータに対して過剰に複雑である

場合に典型的に現れます。 すなわち、この結果は 「ランダム傾きを含めるだけの情報がデータに含まれていない」 ことを示す警告と解釈できます。
 以上を踏まえると、この例3のデータに対してランダム切片・ランダム傾きモデルを当てはめた結果は、

  • 固定効果・ランダム効果ともに、個人内の x–y 関係が存在しないことを正しく反映している
  • 一方で、ランダム傾きを含めたことで推定上の不安定さ(極端な相関など)が生じている

と言えます。
 したがって、このデータ構造に対しては、

  • ランダム切片モデル
\displaystyle{
y \sim x + (1 \mid Subject)
}

あるいは、目的によっては * 切片のみのモデル(あるいは LM)

が、解釈の明確さと安定性の観点からより適切です。
 この例は、「ランダム切片・ランダム傾きモデルは常に万能ではなく、個人内効果が存在しない場合には過剰なモデルになり得る」という重要な教訓を与えてくれます。

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。

関連記事

chaosmemo.com chaosmemo.com