医療統計学を学ぶ大学生のブログ

医療統計学、因果推論を専攻しています。R, SASユーザーです。

SASによる解析:データのシミュレート(線形回帰)

f:id:NorihiroSuzuki:20210701161721p:plain


 ある仮定を持つデータのシミュレートは、手法の頑健性などを調べる際にはよく行います。 この記事では、SAS上での線形関係があるデータのシミュレートと、実際にそれを用いた重回帰分析を行います。

 実行環境はSAS Studioですが、SAS9.4 でも同様に実行可能です。最後にSAS9.4での実行用のコードも載せておきます。

 

 

SAS Studioでのシミュレート

 SAS Studioでは、データのシミュレートのガイドとして、いくつかのスニペットが準備されています。今回はこのうち「線形回帰データのシミュレート」を利用しつつ、線形な関係を持つ変数データの作成、またそれを用いての重回帰分析をやってみます。

 

まずはSAS Studioにログインして、SAS Studioを起動します。

次に、画面左にあるスニペットを展開し、データ、線形回帰分析データのシミュレートを選択します。

f:id:NorihiroSuzuki:20210701143834p:plain

 

すると、以下のように自動的にコードが生成されます。このコードの内容について解説をまず行います。

 

f:id:NorihiroSuzuki:20210701144505p:plain

 

 今回のシミュレートはDATAステップを用いて行われているため、初めにデータセット名である"regdata"が指定されています。

 

 次に、call streaminit(数値)で乱数シードを指定しています。この数値が0でない場合、このコードで作成される乱数は同一になります。それに対して、数値を0とした場合には毎回異なる乱数が発生します。

 

 この発生させる乱数をデータセットに含める必要はないため、その次の行でkeep statementを用いて、x1~x3, yのみをデータセットに残すように指定がされています。

 

 次の行のコードでは、目的変数yに対しての重回帰式に含まれる各パラメータを指定しています。今回シミュレートで考える関係式は以下の通りです。

y = b0 + b1*x1 + b2*x2 + b3*x3 + 誤差項

先ほどの行のパラメータ指定によって、上記の切片、偏回帰係数の値を

b0 = 4,  b1 = 0.8, b2 = 1.2, b3 = 2.4

としているわけです。

 

シミュレートで作成するデータのサイズはその次の行のdoステートメントによって指定がされています。後述するrand関数でx1~x3、誤差項の分布を指定していますが、それをdoステートメントによって何回繰り返すか決めています。(今回は100回)

 

次に各説明変数、誤差項が従う分布を指定します。デフォルトで作成された今回のコードでは次のような分布に各変数は従っています。

x1:平均5、標準偏差0.5の正規分布

x2:平均8、標準偏差0.3の正規分布

x3:平均6、標準偏差0.1の正規分布

epsilon(誤差項):平均0、標準偏差0.8の正規分布

※rand関数で指定できる分布の詳細いについては、この記事の後半にまとめています。 

 

回帰係数、説明変数、誤差項の指定が終わったのでここで目的変数yの計算が行われています。すなわちそれぞれのiについて

y = b0 + b1*x1 + b2*x2 + b3*x3 + 誤差項

という計算が行われています。

 

そしてdoループの最後として、outputステートメントによって、生成したすべての変数を新規オブザベーションとしてregdataに格納しています。

 

作成したー他の一部はこんな感じ。

f:id:NorihiroSuzuki:20210701151237p:plain

 

 

シミュレートデータを用いた重回帰分析

 シミュレートしたデータが、想定していた通りに作られたかどうかを確認するために作成したデータを用いて重回帰分析を行ってみます。

つまり、

b0 = 4,  b1 = 0.8, b2 = 1.2, b3 = 2.4

y = b0 + b1*x1 + b2*x2 + b3*x3 + 誤差項

としてデータを生成したわけですから、これらのパラメータに近い値が推定されればいいわけです。

 

というわけで重回帰分析を行ってみた結果がこちら。

f:id:NorihiroSuzuki:20210701153349p:plain

 

 切片の値が約-1.7と生成した際の4とは大きく離れていますが、これは説明変数の中心化を行わなかったためです。

 偏回帰係数はそれぞれ0.8, 1, 3.6と、 データ作成時のb1 = 0.8, b2 = 1.2, b3 = 2.4 と近い値が出ています。x3の偏回帰係数はやや推定結果のほうが大きく出ていますが、これはサンプルサイズをより大きくすれば解決します。

 

 サンプルサイズを1000に増やした場合の結果はこちら。

f:id:NorihiroSuzuki:20210701154010p:plain

より近い値が出ていることがわかるかと思います。

 

 

rand関数で指定可能な分布の詳細

 関数rand()では以下のような20種類の分布が指定できます。

 

f:id:NorihiroSuzuki:20210701160929p:plain

 

SASコード

今回SAS Studioで行ったシミュレート、解析のコードです。SAS9.4でも同様の結果が得られることを確認済みです。コピペ用にどうぞ。

 

data regdata;
call streaminit(112358);
keep x1-x3 y;
b0 = 4; b1 = 0.8; b2 = 1.2; b3 = 2.4;
do i=1 to 100;
x1 = rand( 'NORMAL',5,0.5 );
x2 = rand( 'NORMAL',8,0.3 );
x3 = rand( 'NORMAL',6,0.1 );
epsilon = rand( 'NORMAL',0,0.8 );
y = b0 + b1*x1 + b2*x2 + b3*x3 + epsilon;
output;
end;
run;

proc reg data=WORK.REGDATA alpha=0.05 plots(only)=(diagnostics residuals
observedbypredicted);
model y=x1 x2 x3 /;
run;
quit;

 

参考

Simulating Linear Regression Data

video.sas.com

support.sas.com