EZR で多重代入法を行う方法

多重代入法を EZR で行う方法はない

途中から、R スクリプトで行っていく

コピペでなんとかなるので、がんばってほしい

はじめに

欠損値（欠測値と同じ）があるデータセットにおいて、推定値にバイアスがかかると言われており、欠損値を補完して解析をしたいという要望は多い

欠損値を補完する方法として、現在主流なのが、多重代入法である

EZR で行う方法を知りたい人は多い

しかし、2024 年 6 月 1 日現在、EZR のメニューで実行はできない

本記事では、EZR の R スクリプト窓に、R スクリプトを書くことで、EZR で多重代入法を行う方法を紹介する

できる限りコピペで済むように、シンプルで必須の内容に留めて紹介する

EZR で多重代入法を行う準備

EZR で多重代入法を行うためには、まず準備が必要である

EZR を起動した後、以下の一行を R スクリプト窓に書く（もしくは、以下をコピペする）

install.packages("mice")

mice というのが、R の多重代入法のパッケージ名である

その mice パッケージのインストールを行わないといけない

これは、1 回だけでよい

この 1 行「install.packages(“mice”)」を書いた後、その行にカーソルをおいたままにして、実行ボタンをクリックする（あるいはキーボードの Ctrl と R を押す）と実行される

すると、以下のような窓が出現する

Japan (Yonezawa) を選択して、OK をクリックする

すると、mice パッケージが PC にインストールされる

これで、インストールが完了する

次に、以下の 1 行を書いて、実行する

library(mice)

これは、mice パッケージを呼び出す指令である

この library は、EZR を起動して mice を使いたいときは毎回、mice を使う前に、1 回だけ実行する

EZR が起動している間は、以後は実行する必要はない、そういう指令である

これで、多重代入法の準備が整った

多重代入で重回帰分析を行う方法

すべての変数が連続データの多変量線形回帰、いわゆる重回帰分析の場合、どのように R スクリプトを書くか紹介する

まず、エクセルなどから、欠損値を含むデータを EZR に取り込む

取り込んだデータセット名が Dataset であるとする

データ型の確認と変更方法

補完の前に、データセット内のデータの型が適切になっていると mice が自動で補完方法を選んでくれる

アクティブデータセット → 変数の操作 → データセット内の変数を一覧するを用いて、データ型を確認する

連続データは num、カテゴリカルデータ（順序なし）は Factor、カテゴリカルデータ（順序あり）は Ord.factor になっていれば、それに応じて適切な補完方法が選ばれる

適切なデータ型になっていない場合、アクティブデータセット → 変数の操作の中の

• 因子あるいは文字列として扱われている数値を連続変数に変換する

• 連続変数を因子に変換する

• 因子水準を再順序化する

を用いて、適切なデータ型に変更する

ただし、連続データに自動で適用される pmm（予測平均マッチング）は、カテゴリカルデータにも適用でき、全データに pmm を適用するのがよい

次節で設定方法を示す

補完方法

そこで、R スクリプト窓に以下の一行を書いて、実行する

imp <- mice(Dataset, m=100, maxit=10, method="pmm", seed=1)

ここで、imp は多重代入した後のデータセットの名前、m は多重代入データセット数、maxit は最大反復回数、method は補完方法（pmm は予測平均マッチングを示す）、seed は再現できるように、ランダム変数のシード（種）を指定したもの

データセット数は、m = 100 とした

100 から 1000 のデータセットが欲しいと言われている

多重代入データセットの数については、こちらも参照のこと

Fraction of Missing Information (FMI) という値の 100 倍を推奨している論文もある

Multiple imputation using chained equations: Issues and guidance for practice

maxit は、最大反復回数の意味だが、10 から 20 回は行うとされる

補完データセット一つにおける代入値を決める計算は、初期値を決めて、計算を繰り返すことで、収束させていくという方法を取るため、反復計算の回数を決める必要がある

ここでは、推奨の最低限の 10 とする（ただし、下記のスクショでは、デフォルトのまま（5 になっている）であるので、その点ご留意を）

method は前節で説明したとおり、連続データのほか、カテゴリカルデータにも適用できる方法で、最適な補完方法である

また、seed は、ここでは、1 としたが、何でもよい

通常は、正の整数にする

出力窓に何やらたくさん出てくるが、多重補完されたことを示している変数名の一覧というだけなので、気にしない

これで、多重補完した 100 のデータセットからなる多重代入データセット imp ができた

例えば、3 番目の補完データセットの先頭部分のみ少し見てみる場合、以下の一行を書くと見られる

head(complete(imp, 3))

ここで、head は先頭の 6 行を表示する指令、complete は、補完したデータを取り出す指令、 3 は、3 番目の補完データセットの指定という意味である

もともと、欠損値があった最初のデータセットと見比べてみると、確かに欠損値が補完されているのがわかる

3 番目ではなく、1 とか、 5 とかの数字にしてみると、ほかの数字で補完されていることがわかるので、見てみてほしい

では、欠損値が補完された 100 個のデータセットを使って、どうやって計算するか

計算方法

計算方法は、以下のように行う

2 行目の最後のカッコまで選択して、実行する

以下の、スクリプトを活用して、変数名を置き換えて活用してほしい

lm_imp <- with(imp, lm(Sales ~ Income + Population +
Price + Age + Education))

ここで、

• lm_imp は、計算結果につける名前なので、何でもよい（ただし、半角英字先頭の半角英数のみがおすすめ。文字は _ と . が使える）
• with は 次に書いてある、ここでは imp という名前の多重補完したデータセットを使う関数
• lm は、単回帰、重回帰を行う関数
• Sales 以降は、回帰のモデル式で、~ の左側が目的変数、右側が説明変数で、+ でつなぐ（記号はすべて半角）

である

with と lm という関数名は変更せず、半角カッコや半角カンマの構造はそのままに、変数名は、ご自分の変数名に置き換えて、試してみてほしい

統合方法

計算結果は、100 個の補完データセットで計算された、100 個の結果が格納されている

これを統合して、結果を表示させるには、以下のように書いて実行する

summary(pool(lm_imp))

すると、以下の赤枠のように、統合された一つの計算結果が出てくる

estimate （偏回帰係数）は、term（説明変数）が、1 上昇したときの、目的変数の推定値の変化を示している

右端の p.value が P 値で、これが、 0.05 未満であれば、estimate が 0 という帰無仮説が棄却されて、説明変数が 1 上昇したときの変化がゼロではなく統計学的に意味があると言えることになる

ちなみに、e-01, e-02 というのは、× 10 の －1乗 すなわち 0.1 倍、× 10 の － 2 乗 すなわち 0.01 倍という意味の科学的表現である

この例の場合、Age と Price と Income が 0.05 未満であり、統計学的有意である

Intercept（切片）は、ほかの説明変数が 0 のときの平均値で、この P 値には、出力されているが注目・解釈する必要はない

ちなみに、元のデータの結果と比べると、同様の結果であることが見て取れる

なので、今回の例では、欠損値がある元データの解析をメインにして、多重代入の結果をサブ（感度解析）というふうに示していくのがよいと思う

異なる場合は、多重代入することで、バイアスが取り除かれたと考えて進めるのもよいだろう

多重代入を行い共分散分析を行う場合

群の比較を背景因子で調整した多変量線形回帰、いわゆる共分散分析（カテゴリカルデータをダミー変数にした重回帰分析と同じ）の場合、どのように R スクリプトを書くか紹介する

補完方法までは、上記と同一である

計算方法と統合方法

計算方法は、以下の一行を R スクリプトに書く

変数名は、各自のデータセット内の変数名に変更してほしい

ancova_imp <- with(imp, lm(Sales ~ factor(US1) +
Price + Age + Education))

Sales に関して、US1（二値のカテゴリカルデータ）の群間比較を、Price、Age、Education の共変量で調整した共分散分析という式である

統合方法も、上記と同様である

summary(pool(ancova_imp))

そうすると、以下のように計算結果が得られる

factor(US1)[T.1] というのは、US1 の 0 を参照カテゴリとして、1 の推定値が計算されていると読む（US1 が 0 と 1 のカテゴリカルデータの場合）

Price、Age、Education を調整しても、US とそれ以外で Sales の平均値の差は、統計学的有意と解釈できる

多重代入でロジスティック回帰分析を行う場合

ロジスティック回帰分析を行う場合は、glm 関数を使う

多重代入までは、上記と同じである

例えば、以下のように書いて、実行する

imp_lr <- mice(Dataset, m=100, maxit=10, method="pmm", seed=1)

文字データのカテゴリカルデータは、数値（例えば、0, 1）に置き換えておくほうがスムーズと思う

解析は、例えば、以下のように書いて実行する

glm_imp <- with(imp_lr, glm(event ~ exposure + age + sex,
family=binomial))

~ の左側に目的変数、右側に説明変数を置く

family=binomial がロジスティック回帰を指定している

最後に、以下のように書いて実行すると、統合結果が得られる

pooled_glm <- pool(glm_imp)
summary(pooled_glm, exponentiate=TRUE, conf.int=0.95)

estimate が オッズ比で、2.5 %, 97.5 % が 95 % 信頼区間の下限と上限である

もともとの結果は以下のとおりで、方向性は同じだが、少々違うことがわかる

多重代入を行い Cox 回帰を行う場合

Cox 回帰（Cox の比例ハザードモデルと同じ）の場合も、補完データ作成は同じ

例えば、以下のように R スクリプト窓に書いて、実行する

imp_cox <- mice(Dataset, m=100, maxit=10, method="pmm", seed=1)

次に、with と coxph を使って、解析を実行する

library(survival)
coxph_imp <- with(imp_cox, coxph(Surv(time, status01) ~
age + sex + ph.ecog012 + meal.cal + wt.loss))

最後に、以下のように統合・結果出力を行う

pooled_cox <- pool(coxph_imp)
summary(pooled_cox, exponentiate=TRUE, conf.int=0.95)

オリジナルの欠測があるデータセットの場合は、以下のようになる

傾向は同じだが、少しだけ数値が異なるのがわかる

t 検定とカイ二乗検定はどうやるか

2 群や 3 群以上の平均値や割合の比較はどのようにすればよいか

平均値であれば、カテゴリカル変数一つのみを説明変数にした、回帰分析を行えばよい

割合であれば、カテゴリカル変数一つのみを説明変数にした、ロジスティック回帰分析を行えばよい

例えば、平均値の場合は、以下のようにする

ttest_imp <- with(imp, lm(Sales ~ factor(US1)))
summary(pool(ttest_imp))

結果が出力されたら、カテゴリカル変数の P 値を確認する

3 カテゴリでも参照カテゴリとの差が検討できる（この例の場合は、ShelveLoc = 1 が参照カテゴリ）

multcomp_imp <- with(imp, lm(Sales ~ factor(ShelveLoc1)))
summary(pool(multcomp_imp))

割合の場合は、以下のように書いて実行する

chisqtest_imp <- with(imp_lr, glm(event ~ exposure,
family=binomial))
pooled_chisqtest <- pool(chisqtest_imp)
summary(pooled_chisqtest)

3 群以上でも解析できて、その場合は、参照カテゴリとの比較結果になる

補完値はどうやって生成しているのか

連続データでも、カテゴリカルデータでも、予測平均マッチング（predictive mean matching, PMM）という方法を使うという話をしてきた

PMM は、線形回帰（重回帰）モデルによる代入方法を基礎にしている

線形回帰モデルによる代入方法は、欠損値を除いた完全データで、欠損値を含む連続データを予測する重回帰式を基にして、補完する方法である

PMM は、その線形回帰モデルで計算される「欠測データの予測値候補」と、実際に観測されている「観測データの完全データからの予測値」の差（距離）が近い K 人を選び、その中からランダムに選んだ人の観測データを補完値とする方法である

K は、5 から 10 とか、3 から 10 とか、が良いと言われており、mice は、5 に設定されていて、変更は不要である

詳細は、野間 2017 参照

以下も、よければ参照

Tuning multiple imputation by predictive mean matching and local residual draws

まとめ

多重代入法を EZR で実施する方法を解説した

多重代入法は、欠損値の多重代入・解析・統合の 3 ステップが一つになった方法である

mice パッケージの mice と pool を使って実施する

重回帰、共分散分析（群間の背景を調整することが目的の重回帰）、ロジスティック回帰、Cox 回帰の方法について解説した

カテゴリカルデータを説明変数にした単回帰、ロジスティック回帰を使えば、t 検定やカイ二乗検定の統合の代わりになることも説明した

何らか参考になれば幸い

参考文献

連鎖方程式による多重代入法

mice: Multivariate Imputation by Chained Equations in R | Journal of Statistical Software

Tuning multiple imputation by predictive mean matching and local residual draws

Multiple imputation using chained equations: Issues and guidance for practice