📚

Rで調整済み平均(EMMean)を算出する

2024/08/09に公開

はじめに

データ解析の中で調整済み平均を調べる機会があったので、その概念について調べて、自分が理解できるように噛み砕いたものをまとめています。間違いがあれば教えてください。

> emmeansパッケージ
> emmeansパッケージのCRAN
> emmeansパッケージを解説する動画(yuzaR Data Science)

調整済み平均とは

調整済み平均は、共変量の影響を取り除いた平均値を算出する方法です。共変量の影響を取り除くことで、グループ間の比較をより正確に行うことができます。

調整済み平均の算出方法としては、 最小二乗平均(LSMean / Least Squares Mean)と推定周辺平均(EMMean / Estimated Marginal Mean)が存在します。LSMeanは線形モデルにのみ適用できる手法で、EMMeanはそれ以外のモデルにも適用できる、より一般化された手法ととらえることができます。

最小二乗平均(LSMean)の算出プロセス

  1. モデル作成:
    目的変数とカテゴリ変数(グループ)および共変量を含む線形モデルを構築する。線形モデルは、各グループと各共変量がそれぞれ目的変数に与える影響を明らかにする。

  2. 共変量の調整:
    全データの共変量の値を、共変量の全体平均に調整する。これにより、グループ間の共変量の違いを極力排除でき、グループの違い(と未知の共変量)による影響のみを目的変数に反映できる。

  3. LSMeanの計算:
    調整した共変量の値を用いて、改めて各グループにおける目的変数の値を算出する。これらの値が、各グループの最小二乗平均(LSMean)となる。

推定周辺平均(EMMean)のアプローチ

  1. モデル作成:
    目的変数とカテゴリ変数(グループ)および共変量を含む一般化線形モデル(GLM)または一般化線形混合モデル(GLMM)を構築する。このモデルは、各グループと各共変量がそれぞれ目的変数に与える影響を明らかにする。EMMeanは、より複雑なモデル構造(例:ネストされた要因、交互作用)や非線形関係も扱える。

  2. 共変量の調整:
    EMMeanでは周辺化(marginalization)という手法を用いて、共変量の効果を平均化する。具体的には、共変量の観測値それぞれについてモデルによる予測を行い、これらの予測値の加重平均を取る。この方法により、線形モデル以外のモデルにも対応しつつ、グループ間の比較において共変量の影響を極力排除できる。

  3. EMMeanの計算:
    平均化された共変量の効果を用いて、改めて各グループにおける目的変数の予測値を算出する。
    これらの予測値が、各グループの推定周辺平均(EMMean)となる。

Rでの実装

Rでは推定周辺平均の算出に特化したemmeansパッケージを使用することで、EMMeanを算出することができます。

✓ "emmeans"パッケージのインストール

# パッケージのインストール・読み込み
if (!require(emmeans)) install.packages("emmeans")
library(emmeans)

# あるいはこちらでも

# 読み込みたいパッケージをベクターで列挙
packages <- c("emmeans")

# 未インストールのパッケージを検出してインストール
new_pkgs <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_pkgs)) install.packages(new_pkgs)

# 列挙されたパッケージをまとめて読み込み
for (pkg in packages) {
    suppressMessages(library(pkg, character.only = TRUE))
}

✓ サンプルデータの作成

# サンプルデータの作成
set.seed(123)
n <- 100
data <- data.frame(
  group = factor(rep(c("A", "B", "C"), each = n/3)),
  covariate = rnorm(n, mean = 50, sd = 10),
  response = numeric(n)
)

# 応答変数の生成(グループと共変量の影響を含む)
data$response <- 10 + 
  (data$group == "B") * 5 + 
  (data$group == "C") * 10 + 
  0.5 * data$covariate + 
  rnorm(n, mean = 0, sd = 5)

✓ 線形モデルの作成(他モデルでも可)

lm()で線形モデルを作成し、emmeans()関数でEMMeanを算出します。

# 線形モデルの作成
model <- lm(response ~ group + covariate, data = data)

# EMMeanの算出
emm_result <- emmeans(model, specs = "group")

# 結果の表示
print(emm_result)

観測したいグループ変数や共変量を増やすこともできます。

# 複数のグループ変数と共変量を含む線形モデルの作成
model <- lm(response ~ group1 * group2 + covariate1 + covariate2, data = data)

# EMMeanの算出(複数のグループ変数を指定)
emm_result <- emmeans(model, specs = c("group1", "group2"))

✓ emm_resultの中身と解釈

emm_resultは以下のように出力されます。(あくまでサンプルです。)

group   emmean    SE  df lower.CL upper.CL
 A      35.1234 1.234 97  32.6789  37.5679
 B      40.1234 1.234 97  37.6789  42.5679
 C      45.1234 1.234 97  42.6789  47.5679

Confidence level used: 0.95

group: 各グループの名称。データフレームのグループ変数のカラム名がそのまま引き継がれる。
emmean: 推定周辺平均(Estimated Marginal Mean)。
SE: 標準誤差(Standard Error)。推定値の精度を示す。
df: 自由度(Degrees of Freedom)。t分布や信頼区間の計算に使用される。
lower.CL upper.CL: 95%信頼区間の下限と上限。真の平均値がこの範囲内にある可能性が95%あることを示す。

おわり

調整済み平均を算出したのち、グループ間の比較等の解析も求められます。
それらはもう少し調べて確証を得てから。

Discussion