多重代入法 (Multiple Imputation, Mice) は、欠損値に対処するための強力な統計的手法である。しかし、多重代入法によって作成された複数のデータセットそれぞれに対して統計解析を行った後、それらの結果をどのように統合すればよいか迷うことがある。本記事では、特にCox比例ハザード回帰モデルにおけるWald検定（ANOVA）のカイ二乗値を、rmsパッケージのcph関数とmiceaddsパッケージのmicombine.chisquare関数を用いて統合する方法について解説する。
どんなときに使う方法か
この方法は、以下のような状況で特に役立つ。
- データに欠損値があり、多重代入法を用いて欠損値を補完したケース。
- イベント発生までの時間（生存時間）をアウトカムとするCox比例ハザード回帰モデルを構築したい場合。
- モデルに含まれる特定の変数（カテゴリカル変数など）の全体的な効果を評価するため、Wald検定（ANOVA）のカイ二乗値を利用したい場合。
- 多重代入法によって得られた複数のカイ二乗値を適切に統合し、単一の信頼できるp値を導き出したい場合。
どのような方法か
多重代入法では、m 個の代入済みデータセットが生成される。それぞれのデータセットに対してCox回帰モデルを当てはめ、Wald検定（ANOVA）のカイ二乗値を算出する。これらの m 個のカイ二乗値は、直接平均をとるだけでは適切に統合できない。
ここで登場するのが、Rubinの統合ルールの考え方に基づいたカイ二乗値の統合である。miceaddsパッケージのmicombine.chisquare関数は、複数のカイ二乗値とそれぞれの自由度を用いて、プールされたカイ二乗値とプールされた自由度、そしてp値を算出する。この関数は、特にANOVAタイプのカイ二乗検定の統合に適している。
rmsパッケージのcph関数は、Cox回帰モデルを柔軟に指定できる関数であり、モデルの当てはめ後にanova()関数を用いることで、各変数のWald検定（ANOVA）のカイ二乗値と自由度を得ることができる。
具体例
ここでは、架空のデータセットを用いて具体的な手順を説明する。
シナリオ: ある疾患の患者の生存時間を予測するモデルを構築したいと考えている。患者の年齢、性別、治療法（A, B, C）が生存時間に影響を与えるかを検討する。ただし、データには欠損値が含まれているとする。
目的: 治療法が生存時間に与える全体的な効果を、多重代入法を適用したデータでCox回帰を行い、Wald検定（ANOVA）のカイ二乗値を統合して評価する。
具体例のRでの計算
R 計算例：
# 必要なパッケージのインストール (初回のみ)
# install.packages("mice")
# install.packages("rms")
# install.packages("miceadds")
# install.packages("survival") # cphの内部で使用されるため
# パッケージのロード
library(mice)
library(rms)
library(miceadds)
library(survival)
# データの作成 (例として欠損値を含むデータを作成)
set.seed(123)
n <- 100
df_orig <- data.frame(
time = rexp(n, rate = 0.01),
status = sample(0:1, n, replace = TRUE, prob = c(0.3, 0.7)),
age = rnorm(n, 60, 10),
sex = factor(sample(c("Male", "Female"), n, replace = TRUE)),
treatment = factor(sample(c("A", "B", "C"), n, replace = TRUE, prob = c(0.4, 0.3, 0.3)))
)
# 欠損値をランダムに導入
df_orig$age[sample(1:n, 10)] <- NA
df_orig$sex[sample(1:n, 5)] <- NA
df_orig$treatment[sample(1:n, 8)] <- NA
# 多重代入法の実行
# method = "pmm" (Predictive Mean Matching) は連続変数、"polyreg" はカテゴリカル変数に適用
imp <- mice(df_orig, m = 5, seed = 456, printFlag = FALSE)
# 各代入済みデータセットでCox回帰を実行し、Wald検定（ANOVA）のカイ二乗値と自由度を抽出
chi_squares <- list()
dfs <- list()
for (i in 1:imp$m) {
# i番目の代入済みデータセットを取得
df_complete <- complete(imp, i)
# Cox回帰モデルの構築 (rms::cphを使用)
# survオブジェクトを作成
surv_obj <- Surv(df_complete$time, df_complete$status)
# cphモデルの当てはめ
# x=TRUE, y=TRUE は後でanova()を使うために必要
fit_cph <- cph(surv_obj ~ age + sex + treatment, data = df_complete, x=TRUE, y=TRUE)
# treatment変数のANOVA検定結果を取得
# anova(fit_cph)はデータフレームを返す
anova_result <- anova(fit_cph)
# treatmentの行を取得 (変数名 'treatment' が含まれる行)
# anova_resultの行名は変数名になっていることが多い
# exactなマッチング 'treatment' を行う
treatment_row <- anova_result[rownames(anova_result) == "treatment", ]
# カイ二乗値と自由度を抽出
# カイ二乗値は 'Chi-Square' 列、自由度は 'd.f.' 列に格納されている
chi_squares[[i]] <- treatment_row["Chi-Square"]
dfs[[i]] <- treatment_row["d.f."]
}
# リストからベクトルに変換
chi_squares_vec <- unlist(chi_squares)
dfs_vec <- unlist(dfs)
# 統合されたカイ二乗値を計算
# micombine.chisquare関数を使用
combined_result <- micombine.chisquare(
dk = chi_squares_vec,
df = dfs_vec[1]
)実行結果と解釈
上記のRコードを実行すると、combined_resultオブジェクトに統合されたカイ二乗検定の結果が格納される。出力は以下のようになるだろう。
> combined_result <- micombine.chisquare(
+ dk = chi_squares_vec,
+ df = dfs_vec[1]
+ )
Combination of Chi Square Statistics for Multiply Imputed Data
Using 5 Imputed Data Sets
F(2, 46.03)=1.969 p=0.1512
> F(2, 46.03): 統合されたカイ二乗統計量である。カッコ内の数字は自由度である。p-value: 統合されたカイ二乗統計量と自由度に基づいたp値である。このp値は、treatmentという変数全体が、他の変数を調整した後でも生存時間に有意な影響を与えるかどうかを示す。
この例では、p-valueが0.05を下回っていないので、統計的に有意な結果とは言えない。これは、治療法が生存時間に有意な影響を与えるとは言えないことを示唆している。
まとめ
多重代入法後のCox回帰において、特定の変数（特にカテゴリカル変数など複数の自由度を持つ変数）の全体的な効果を評価する際には、Wald検定（ANOVA）のカイ二乗値を統合する必要がある。
本記事で紹介したrmsパッケージのcph関数とmiceaddsパッケージのmicombine.chisquare関数を組み合わせることで、このプロセスを効果的に実行できる。これにより、多重代入法による欠損値補完の恩恵を受けつつ、統計的に適切な方法でモデルの評価を行うことが可能となる。
この方法は、臨床研究や疫学研究など、欠損値が頻繁に発生する分野でCox回帰モデルを扱う研究者にとって非常に有用である。
コメント