Rで調整済み平均(EMMean)を算出する
はじめに
データ解析の中で調整済み平均を調べる機会があったので、その概念について調べて、自分が理解できるように噛み砕いたものをまとめています。間違いがあれば教えてください。
> emmeansパッケージ
> emmeansパッケージのCRAN
> emmeansパッケージを解説する動画(yuzaR Data Science)
調整済み平均とは
調整済み平均は、共変量の影響を取り除いた平均値を算出する方法です。共変量の影響を取り除くことで、グループ間の比較をより正確に行うことができます。
調整済み平均の算出方法としては、 最小二乗平均(LSMean / Least Squares Mean)と推定周辺平均(EMMean / Estimated Marginal Mean)が存在します。LSMeanは線形モデルにのみ適用できる手法で、EMMeanはそれ以外のモデルにも適用できる、より一般化された手法ととらえることができます。
最小二乗平均(LSMean)の算出プロセス
-
モデル作成:
目的変数とカテゴリ変数(グループ)および共変量を含む線形モデルを構築する。線形モデルは、各グループと各共変量がそれぞれ目的変数に与える影響を明らかにする。 -
共変量の調整:
全データの共変量の値を、共変量の全体平均に調整する。これにより、グループ間の共変量の違いを極力排除でき、グループの違い(と未知の共変量)による影響のみを目的変数に反映できる。 -
LSMeanの計算:
調整した共変量の値を用いて、改めて各グループにおける目的変数の値を算出する。これらの値が、各グループの最小二乗平均(LSMean)となる。
推定周辺平均(EMMean)のアプローチ
-
モデル作成:
目的変数とカテゴリ変数(グループ)および共変量を含む一般化線形モデル(GLM)または一般化線形混合モデル(GLMM)を構築する。このモデルは、各グループと各共変量がそれぞれ目的変数に与える影響を明らかにする。EMMeanは、より複雑なモデル構造(例:ネストされた要因、交互作用)や非線形関係も扱える。 -
共変量の調整:
EMMeanでは周辺化(marginalization)という手法を用いて、共変量の効果を平均化する。具体的には、共変量の観測値それぞれについてモデルによる予測を行い、これらの予測値の加重平均を取る。この方法により、線形モデル以外のモデルにも対応しつつ、グループ間の比較において共変量の影響を極力排除できる。 -
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