欠損値があるデータを扱うＭＩプ
ロシジャについて小野裕亮
テクニカルサポートグループ（株）ＳＡＳインスティチュートジャパン
jpnyco@jpn.sas.com

Copyright ©2000, SAS Instit ute Inc. All r ights r eser ved.

内容：
multiple imputation （多重代入；以後ＭＩ）を、ＳＡＳで行なう方法を紹介→プログラミングが中心。
×理論 ×実務 （？？？）■ 1） ＭＩを行なうためのソフトウェア必要要件■ 2） Version6では？
■ 3) MIの一般論■ 4） PROC MIについて■ 5） PROC MIANALYZEについて1

ソフトウェア必要要件＊MI
MIを行なうプロシジャ＊を行なうプロシジャ＊＊バージョン8.1～＊バージョン ～＊SAS/STATプロダクト＊評価版（experimental）
バージョン6には存在しない。
■日本：バージョン8.1 （日本語版）を、■Microsoft Windows版のみリクエストに応じて出荷中→
http://www.sas.com/japan/service/v8/index.html■１１月の現状：マニュアルなどは英語,■
Base SASの一部分メッセージが日本語。
■V6における欠損値（．）の扱いにおける欠損値（．）の扱い■ほとんどのプロシジャはリストワイズで削除

x１
10.5
14.2．
12.8
10.4x2x3
19.1 30.9
17.5 35.5
20.1 35.4
20.4 35.2
15.2 33.82

SAS V6における現状における現状A)単一代入（single imputation)■• 何らかの値で、事前にとりあえず穴を埋める努力• 1)STANDARDプロシジャ（Base SAS） … 平均。Easy!
• 2)ユーザ自身のマクロプログラム、データステップなどによる涙ぐましい努力– 例）Hot-Deck Imputation … 米国SUGI 1998, USA Census

B)CORRプロシジャ (Base SAS)・デフォルトは、ペアワイズで削除C)MIXEDプロシジャ（SAS/STAT）…Observed-data likelihood■

SAS V6における現状における現状D)打ち切りデータ（生存時間モデル・トービットモデル）

Ｅ）ＩＭＬによって自分自身でプログラミング。。。。。。
例）ペンシルバニア大学のPaul Allison...

3

余談：ＳＡＳデータセットの欠損値（「欠損」といってもいろいろあるが）ＳＡＳデータセットの数値変数における欠損値
基本 .
（ピリオド）
特殊欠損値 .A .B ... （ピリオド＋アルファベット）
分析において、この違いを積極的に使うものは特にない（例外）
MEANSのCLASSステートメント,FREQの単純集計
UNIVARIATEの欠損値数出力
TRANSREG,PRINQUALプロシジャ

記号：（ＳＡＳマニュアル） Shafer ?
Yobs ...観測されたデータ
Ymis …観測されなかったデータ
Ｒ … 欠損値の位置を表すインデックス
Θ … 完全データに対するモデルパラメータ
ξ … 欠損メカニズムを決定するパラメータ例） X２＝β0 + β1 X1 + εPr(X1=.)=p1
Pr(X2=.) =p2
モデルパラメータ θ = ( β0, β1, σ)
欠損メカニズムのパラメータ ξ=(p1 p2)

4

欠損値いろいろ１）Missing Completeｌｙ At Random
Yobs, Ymisに依存せず、欠損がランダムに生じている

２）Missing At Random （ランダムな欠損）
Yobsだけに依存して欠損がランダムに生じている。
Pr（Ｒ｜Yobs, Ymis, ξ） = Pr（Ｒ｜Yobs,ξ）
条件や方法によっては、欠損が無視できる… ...

３）nonignorable （欠損がinformativeな場合）
欠損のメカニズムをモデル化する必要

Multiple Imputation, Observed-data Likelihoodの前提：の前提：
欠損を無視できる？
Ｌ（θ| Yobs)が最大になるようなθ｜｜Ｌ（θ,ξ｜Ｙobs , Ｒ） が最大になるようなθ,ξ
１） Missing At Random (or MCAR）
２） θとξがdistinct （ベイズ流：θとξが独立）。
・分析者は、ξには特に関心ない。
・分析者は、完全データ上のモデルパラメータ θに関心

5

喩え話：
２正規変量X1
において、２正規変量 と X2において、X2にに欠損値が存在・A) MCAR
サイコロを振って、Ｘ１,Ｘ２の一部分を観測しない。
・B) MAR （Yobsのみに依存）。
Ｘ１＞１以上の場合には、サイコロを振ってＸ２の値を観測しない。
・C) Ymisに依存。
Ｘ２ ＞ １の場合には、サイコロを振ってＸ２の値を観測しない。

平均と分散共分散行列の計算研究者の興味は、μ研究者の興味は、μとΣA）
リストワイズの削除でもＯＫ。
B)
リストワイズの削除→ ×Ｘ２の平均および分散, 共分散, 相関×C)
リストワイズの削除→ × Ｘ２の平均および分散, 共分散, 相関×

6補足：
ＲＥＧ：Ｘ２＝β＋ε：
ＲＥＧ：Ｘ２＝β0 + β1 X1＋ε：
研究者の興味はβ
研究者の興味はβ0, β1, σ- A)
リストワイズの削除でもＯＫ （推定効率,検出力 ↓）- B)
リストワイズの削除でもＯＫ推定効率,検出力 ↓。標準化偏回帰係数, R2乗値は×- C)
リストワイズ削除： β0,β1,σの推定値とも×。
打ち切りを考慮して、尤度関数を定義＝LIFEREG補足：
ＲＥＧ：Ｘ１＝β２＋ε：ＲＥＧ：Ｘ１＝β0 + β1 X２＋ε：
研究者の興味はβ
研究者の興味はβ0, β1, σ-A)
リストワイズの削除でもＯＫ（推定効率,検出力 ダウン）。
-B)
リストワイズ削除： β0,β1,σの推定値とも×。
-C)
リストワイズ削除： β0,β1,σの推定値ともＯＫ（推定効率,検出力 ダウン、標準化偏回帰係数、Ｒ２乗値は×）7

V6の機能：の機能：MIXEDプロシジャ
プロシジャの機能：“Observed-data likelihood”をを最大化“Observed-data Likelihood ” を最大化 →Ｌ（θ,ξ｜Ｙobs , Ｒ） ではなく、Ｌ（θ| Yobs)を最大化。
? max( Ｌ（θ,ξ｜Ｙobs , Ｒ）)→ max ( Ｌ（θ| Yobs)) ?
本当は、max（ Ｌ（θ| Ymis, Yobs)）がベスト。

V6の機能：の機能： MIXEDプロシジャ
プロシジャ
多変量分析が目的ではないので…
多変量分析が目的ではないので 。
proc transpose data=data1 name=var out=out1;
by id;
run;proc mixed data=out1;
class id var;model col1=var /ddfm=kenwardroger solution ;repeated var /subject=id type=un r rcorr;run;（注：小標本のときの振舞い）

8

ＭI推定を行なう手順
推定を行なう手順（２つに分類される）ＭＩ実行者の任務… “モデルＡ”に従う乱数で穴埋めを実行。この処理を複数回実行し、穴埋めされた複数個のデータセットを作成。
分析者の任務… 各々の完全化されたデータセットを、“モデルＢ”にあてはめる。穴埋めによるバラツキも考慮して、複数の結果をまとめる。

MI推定のＳＡＳプログラム
推定のＳＡＳプログラム… ３ステップに分けられる。
ＭＩ実行者
１）PROC MIを用いる。（重要任務：モデルＡの選択）
分析者
２）既存プロシジャでモデルＢをあてはめる。
3) 結果をまとめるためPROC MIANALYZEを使用。

9

ステップ１：
ＭＩプロシジャを実行する。
proc mi data=data1 out=out1 nimpu=70
seed=654321;
multinormal method=MCMC( prior=jeffreys
initial=EM(MLE)
biter=50
chain=multiple);var x1 x2 x3;
run;

ステップ２：
分析者は既存プロシジャを実行する。
PROC REG DATA=WORK . OUT1OUTEST=WORK.OUT2 COVOUT ;BY _IMPUTATION_;MODEL X3=X1 X2 ；
ＲＵＮ；

10

ステップ３：分析者は、さらにMIANALYZEプロシジャを実行。
プロシジャを実行。

proc mianalyze data=work.out2 edf=97;
var Intercept x1 x2;
run;

ステップ１：PROC
ステップ１：
MI
代入モデルを選ぶ。
PROC MIで用意されている方法は、全部で３つ。
MULTINORMAL METHOD= ＊（1）ロジスティックモデルに基づく方法 PROPENSITY（propensity score method)
×（2） 回帰モデル法 REGRESSION Version8.1×（３）多変量正規をMCMCで発生させる方法 MCMC（制限） 1,2は、欠損が単調(monotone)である時しか適用できない。

11Monotoneなな 欠損構造
T1
X
X
X
X

T2
X
X.
X

T3 T4.
.
X X.
.
X .
＊前が欠損だったら、後ろも欠損。

（１）Propensity（１）
Score Method
METHOD=PROPENSITY（a）「欠損値 or 非欠損値」の２値を従属変数としてロジスティック（b）ロジスティックモデルの予測値が近い観測値を抽出して、欠損値を穴埋め。
－欠損値 or 非欠損値の情報しか用いていない。
付随するオプション：
METHOD=PROPENSITY（ GROUP= n ）－何グループにするか？

12

（2）回帰モデルに基づく方法）回帰モデルに基づく方法
METHOD=REGRESSION
回帰モデルに基づき、穴埋めを行なっていく。
欠損構造がmonotoneの時のみ。
でも、ＭＣＭＣ法よりも高速。
基本的に、多変量正規分布。
・非欠損値の部分から、βおよびσの推定値→β* , σ*を乱数で生成→ Ｙ＝Ｘβ* + σ* z で欠損値部分を埋める。

（2）回帰モデルに基づく方法）回帰モデルに基づく方法
METHOD=REGRESSION；；

重要：重要：Version8.1では、METHOD=REGRESSIONは常に間違った結果になっています。Version8.2で、このバグは修正されています。

13(3)MCMC法
法 METHOD=MCMC
多変量正規分布“Impute”-step と“Posterior”-stepを交互に行なう。
（Ｙmis｜μ,Σ,Yobs) →（μ,Σ｜Y）→（μ,Σ|Ｙ）→ ….

→ （Ｙmis｜μ,Σ,Yobs)＜オプション＞
ＩＮＩＴＩＡＬ＝ＥＭ（ BOOTSTRAP = p
MLE (デフォルトは事後分布）
ＣＨＡＩＮ＝ＳＩＧＬＥ｜ＭＵＬＴＩＰＬＥ 単鎖 or 複鎖
PRIOR= 事前分布（デフォルトは、Jeffreys)
BITER＝ 各連鎖のまえに、何回、回しておくか？
ITER＝ 単鎖のときの間隔

ＰＲＯＣ ＭＩステートメントのオプション
幾つのデータ（完全化されたデータ）を作成するか？
NIMPU= k ;乱数系列のシード値 SEED=
代入値の最大, 最小 MAX= MIN=
用いるデータ名 DATA=
作成する穴埋めデータ名 OUT=
結果を表示しない NOPRINT

14ＰＲＯＣ ＭＩＡＮＡＬＹＺＥ推定値および推定値間の分散共分散行列を与えなければいけない。
入力データの形式は２通り。
・ＤＡＴＡ＝_TYPE_=“EST” & _TYPE_=“COV”
・PARMS= 推定値
COVB=
分散共分散行列出力される情報： （別紙）☆MI推定値＝ 各推定値の平均☆欠損値を埋めることに伴う変動
・Within-impute Variance (W)・Between-impute Variance (B) →欠損なしの時０・Total Variance ＝ W + B /m + B☆欠損がＭＩ推定値の変動にどれだけ影響しているか？
－ a)fraction of missing information about Q－ b)relative increase in variance due to nonresponse15

今後の予定（11/28 現在）バージョン8.2で、変更・拡張を行なった。マニュアルも、
V8.1のPROC MIは46ページだが、現在、V8.2は72ペー
ジ。
バージョン8.2でも、MIおよびMIANALYZEは、「評価版」。バージョン8.2の次バージョンでは、プロダクト版となる予定。
なお、カテゴリーデータに対する処理は、今後の取り組むべき課題。

Version8.2にに追加される機能（ＭＩプロシジャ）
1. EMステートメント
MIを行なわずに、EMアルゴリズムで（μ,Σ）の推定だけを行ないたい時。
2. TRANSFORMステートメント データを変換。
3. Monotone MCMC method （monotone missing patternのデータを作成する）
4. ＭＣＭＣ法において、定常になったかどうかをチェックするための自己相関プロット。

16

バージョン8.1のバージョン のドキュメント （英語）- ＭＩプロシジャ （Version8.1)
http://www.sas.com/service/techsup/faq/stat_proc/miproc.html- MIANALYZEプロシジャ (Version8.1)
http://www.sas.com/service/techsup/faq/stat_proc/mianalyzeproc.html（ＰＤＦ形式。Adobe社が無料提供しているAcrobat
Readerが必要です）。

アウトプット例：別紙）

17

18