*------------------------------------------------------------------; * 昭和大学 第13回実践臨床統計学セミナー・ハンズオン2 ; * 多重代入法による欠測データの統計解析 ; * 2021年1月18日   ; *------------------------------------------------------------------; *------------------------------------------------------------------; * 1. CSVデータファイルの読み込み               ; *------------------------------------------------------------------; proc import out=work.ovc datafile = "C:\SAS2\ovariancancer.csv" dbms = csv replace; getnames = yes; datarow = 2; run; *-- datafile=... のフォルダ・ファイルのパスは適当に変更してください。; *-- Cドライブ直下に "SAS2" というフォルダを作っていただくと、以後の実習でも、プログラムの改変なしに、配布したコードでデモを実行することができます。; *------------------------------------------------------------------; * 2. MEAN, FREQ プロシジャ: 欠測値の集計           ; *------------------------------------------------------------------; proc means data=ovc; var age log_ca125 log_alp; run; *-- 通常のMEANプロシジャです。連続変数の要約統計量をまとめて計算してくれます。 proc means data=ovc nmiss; var age log_ca125 log_alp; run; *-- MEANプロシジャで、nmissというオプションをつけると、欠測データがどの程度の頻度であったかの集計をとってくれます。; *-- 雑誌によっては、欠測値の割合についての集計結果が求められることもありますが、この機能を使えば、簡単に計算することができます。; proc freq data=ovc; tables figo grade histol ascites ps resdis / missing; run; *-- FREQプロシジャでも、missingというオプションをつけると、欠測データがどの程度の頻度であったかの集計をとってくれます。; *------------------------------------------------------------------; * 3. PHREGプロシジャ: Cox回帰分析での多変量解析を行います          ; *------------------------------------------------------------------; proc phreg data=ovc; class figo(ref="I") grade(ref="well differentiated") histol(ref="Serous papillary") ascites(ref="absent") ps(ref="0") resdis(ref=">5 cm"); model t*d(0) = age figo grade histol ascites ps resdis log_ca125 log_alp/ rl; run; *-- PHREGプロシジャは、Cox回帰分析のためのプロシジャです。; *-- 欠測を含むデータに対して、特に、何の操作も行わずに、PHREGをかけると、「結果変数・説明変数のいずれか1つでも欠測した対象者を除外した解析(Complete Case Analysis)」が行われます。; *-- MODELステートメントで、「結果変数=説明変数」の構文で、回帰式の指定を行います。; *-- 結果変数には、「イベントまでの時間」と「イベント/打ち切りの指示変数」の2つの変数の指定を行う必要があります。; *-- rlというオプションを加えることで、ハザード比の信頼区間を計算することができます。 *------------------------------------------------------------------; * 4. MIプロシジャ: 多重代入法による欠測値の補完を行います ; *------------------------------------------------------------------; proc mi data= ovc nimpute=20 out=mi_ovc; class figo grade histol ascites ps resdis; var age figo grade histol ascites ps resdis log_ca125 log_alp; fcs regpmm(log_ca125) regpmm(log_alp) logistic(figo) logistic(grade) logistic(ascites) logistic(ps) logistic(resdis) ; run; *-- MIプロシジャは、多重代入法による欠測データの補完値を生成するプロシジャです。; *-- nimputeは、補完値の組の数を指定する引数です。通常は100以上に設定されますが、計算時間がかかるので、ここでは、実習のために20と設定しています。任意に設定可能です。; *-- VARステートメントで、今回の解析に使う変数すべてを指定します。互いの変数の情報を利用して、補完値を生成することになります。; *-- FCSステートメントで、MICE (multiple imputation by chaned equation) による解析を行うことを指定します。; *-- regpmm, logisticは、FCSによる補完値の生成において、個々の変数の補完値をどのような方法で生成するかを指定する引数です。; *-- regpmmは、連続変数の補完方法で、予測平均マッチング法という方法になります。; *-- logisticは、2値変数の補完方法で、ロジスティック回帰モデルによる生成を行う方法になります。; proc print data=ovc (obs=20); var age figo grade histol ascites ps resdis log_ca125 log_alp; run; *-- 元のデータセットのはじめの20名分のデータを出力します。これだけでも、かなりの割合で欠測が起こっていることがわかります。; proc print data=mi_ovc (obs=20); var _imputation_ age figo grade histol ascites ps resdis log_ca125 log_alp; run; *-- 1組目の補完後のデータセットのはじめの20名分のデータを出力します。すべての欠測データが埋められていることがわかります。; proc print data=mi_ovc (firstobs=1190 obs=1209); var _imputation_ age figo grade histol ascites ps resdis log_ca125 log_alp; run; *-- 2組目の補完後のデータセットのはじめの20名分のデータを出力します。1組目のデータセットとは異なる補完値が埋められていることがわかります。; *-- このように、補完値の予測の不確実を評価するため、100-1000組ほどの補完値を生成し、そのばらつきを最終的な解析で考慮することになります。; *------------------------------------------------------------------; * 5. PHREGプロシジャ: 補完後のデータの解析を行います ; *------------------------------------------------------------------; proc phreg data=mi_ovc outest=mi_phr covout; class figo(ref="I") grade(ref="well differentiated") histol(ref="Serous papillary") ascites(ref="absent") ps(ref="0") resdis(ref=">5 cm"); model t*d(0) = age figo grade histol ascites ps resdis log_ca125 log_alp/ rl; by _imputation_; run; *-- 補完後の20組のデータセットを、PHREGプロシジャで解析します。; *-- 4行目の by _imputation_ というコマンドで、「補完値の組ごとに、それぞれここで指定したPHREGでの解析を行う」ことを指定することになります。; *-- 1行目の outest=mi_phr というオプションで、「mi_phrという名前のデータセットに、20組のCox回帰の結果を出力する」という指定をしています。; *-- 1行目の covout も必須です。標準誤差を計算するために必要になる分散の推定値を、mi_phrに一緒に出力するという指定になります。; *-- それぞれの補完後のデータセットに対するCox回帰分析の結果も出力されます。; *-- 補完値のばらつきの分だけ、データセットが異なりますので、Cox回帰分析の結果も異なります。; *-- この推定値のばらつきが、「欠測データの予測の不確実性」を反映したものとなります。; *-- 最終的には、この推定値のばらつきを、信頼区間やP値の計算に反映する必要があります。; *------------------------------------------------------------------; * 6. MIANALYZEプロシジャ: 補完後のデータの解析結果の統合を行います ; *------------------------------------------------------------------; proc mianalyze data=mi_phr; modeleffect age figoII figoIII figoIV grademoderately_differentiated gradepoorly_differentiated histolAdenocarcinoma histolEndometrioid histolMesonephroid__clear_cell_ histolMixed_mesodermal histolMucinous histolUndifferentiated ascitespresent ps1 ps2 ps3_or_4 resdis2_5_cm resdis_2_cm log_ca125 log_alp ; run; *-- MIANALYZEプロシジャは、多重代入法による補完後のデータの解析結果を統合するためのプロシジャです。; *-- dataのところで、先ほどの擬似的な完全データの解析によって出力したデータセット(mi_phr)を指定します。; *-- modeleffect で指定した変数について、統合後の回帰係数の推定値・信頼区間・P値を計算します。; *-- 連続変数については、元のデータセットでの変数名をそのまま指定すればOKです。; *-- カテゴリカル変数については、カテゴリごとに変数が定義し直される(カテゴリ間の比較を行うためのダミー変数が生成される)ため、それぞれのカテゴリごとに変数名が再定義されています。; *-- カテゴリごとの変数名は、先ほどのPROC PRINTでのCox回帰分析の結果で出力される表に出されています。少し面倒ですが、推定値・信頼区間が欲しい変数については、表中の変数名をここにコピー&ペーストしてください。; *------------------------------------------------------------------; * 6a. MIANALYZEプロシジャ: 回帰パラメータの推定値をEXP変換して出力する場合 ; *------------------------------------------------------------------; proc mianalyze data=mi_phr; modeleffect age figoII figoIII figoIV grademoderately_differentiated gradepoorly_differentiated histolAdenocarcinoma histolEndometrioid histolMesonephroid__clear_cell_ histolMixed_mesodermal histolMucinous histolUndifferentiated ascitespresent ps1 ps2 ps3_or_4 resdis2_5_cm resdis_2_cm log_ca125 log_alp ; ods output parameterestimates=mi_est; run; *-- MIANALYZEプロシジャは、多重代入法による補完後のデータの解析結果を統合するためのプロシジャです。; *-- dataのところで、先ほどの擬似的な完全データの解析によって出力したデータセット(mi_phr)を指定します。; *-- modeleffect で指定した変数について、統合後の回帰係数の推定値・信頼区間・P値を計算します。; *-- 連続変数については、元のデータセットでの変数名をそのまま指定すればOKです。; *-- カテゴリカル変数については、カテゴリごとに変数が定義し直される(カテゴリ間の比較を行うためのダミー変数が生成される)ため、それぞれのカテゴリごとに変数名が再定義されています。; data HR; set mi_est; HR=exp(estimate); LCL_HR=exp(LCLMean); UCL_HR=exp(UCLMean); pvalue=Probt; run; proc print; var parm HR LCL_HR UCL_HR pvalue; run; *-- 上記の出力した回帰パラメータの推定結果を編集するコマンドを使えば、そのままハザード比の推定値,信頼区間を出力することができます。;