僕の頁 <SASと臨床試験と雑談と>

徒然なるままにSAS暮らし

----

スポンサーサイト  

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Posted on --/--/-- --. --:-- [edit]

CM: --
TB: --

0308

IMLによる多変量正規乱数の生成  

先日、DATAステップによる2変量正規乱数の生成方法を紹介しましたが、
今回はIMLプロシジャを用いた多変量正規乱数の生成方法を紹介します。
多変量正規乱数の生成方法は色々ありますが、ここではIMLプロシジャの
randnormal関数を使用します。

以下の多変量正規(Multivariate Normal; MVN)分布を想定します。
多変量正規分布
randnormal関数は引数に発生させるデータ数、想定する平均値のベクトルと
分散共分散行列を指定します。
randnormal(N,MEAN,Cov.)

分散共分散行列をそのまま指定してもいいですが、現実的には変数間の
相関行列と標準偏差(または分散)から分散共分散行列を計算して指定
することが多いかなと想像してます。今回のプログラムの流れは、

3変数の相関行列、平均値と標準偏差のベクトルを指定し、分散(_VAR)を算出
⇒分散共分散行列を計算(相関行列の各要素と_VAR`*_VARの各要素を#で掛ける)
⇒randnormal関数で乱数を生成
⇒生成したデータの要約統計量と相関係数を確認する

proc iml ;
call randseed(100) ;
create RES var{X1 X2 X3} ;
*--- 各パラメータの設定 ;
_CORR = {1 0.3 0.5,
0.3 1 0.6,
0.5 0.6 1 } ; *--- 相関行列 ;
_MEAN = {5 3 1} ; *--- 各変数の平均値を格納したベクトル ;
_SD = {3 2 1} ; *--- 各変数の標準偏差を格納したベクトル ;
_VAR = _SD##2 ; *--- 各変数の分散を格納したベクトル ;
_COV = _CORR # sqrt(_VAR` * _VAR) ; *--- 分散共分散行列 ;

*--- 3変量正規乱数の生成 ;
N = 500 ;
_X = randnormal(N,_MEAN,_COV) ; *--- データ数, 平均値ベクトル, 分散共分散行列 ;
print _X ;
*--- 結果をデータセットに格納 ;
append from _X ;
close RES ;
*--- 結果の確認(Pearsonの積率相関係数と要約統計量) ;
submit ;
proc corr data=RES pearson ;
var X: ;
run ;
endsubmit; *--- セミコロンの前にスペースを入れるとなぜかerror ;

quit ;

corrプロシジャの結果を見てみましょう。想定通りの相関係数と
平均値・標準偏差のデータが生成されていることが確認できます。
多変量_IMLの結果
スポンサーサイト

Posted on 2014/03/08 Sat. 17:23 [edit]

CM: 0
TB: 0

0224

IMLプロシジャによる2標本t検定の実行  

IMLプロシジャで2標本t検定を実行してみます。何でこんなことを
する必要があるかと言いますと、何万回もシミュレーションをする
必要がある際にttestプロシジャでは時間がかかるからです。今日は
シミュレーションではなく、単純にデータが与えられたもとで、IML
プロシジャで2標本t検定を実行します。

sashelp.BWeightというデータから、母親の喫煙の有無で乳児の
出生時体重に違いがあるかを2標本t検定で検証してみます。

useとreadで各群のデータをそれぞれ読み込んでX1とX2という行列に
格納します。標本平均およびvar関数で不偏分散を算出し、t統計量と
P値を算出します。

・t統計量(m,nは各群の症例数)
t統計量2

実際のプログラムは以下になります。
proc iml ;
use sashelp.Bweight ;
read all var{WEIGHT} where(MOMSMOKE=0) into X1 ; *--- 喫煙無し群 ;
read all var{WEIGHT} where(MOMSMOKE=1) into X2 ; *--- 喫煙有り群 ;
m = nrow(X1) ; *--- 喫煙無し群の症例数 ;
n = nrow(X2) ; *--- 喫煙有り群の症例数 ;
DF = m + n - 2 ; *--- 自由度 ;
DIFF = X1[:] - X2[:] ; *--- 平均値の差 ;
SE = sqrt( ((m-1)*var(X1) + (n-1)*var(X2))/(m + n - 2)*(1/m + 1/n) ) ; *--- SE ;
T = DIFF/SE ; *--- t統計量 ;
if T < 0 then P = 2*probt(T,DF) ;
else P = 2*(1-probt(T,DF)) ;
*--- 両側P値 ;
print T P ;
quit ;

実行結果は以下になります。
IML結果

ちなみにttestプロシジャの実行結果は以下になります。
ttestプロシジャの結果

後日シミュレーションプログラムも紹介したいと思います。

Posted on 2014/02/24 Mon. 23:56 [edit]

CM: 0
TB: 0

0222

2変量正規乱数の生成(DATAステップ)  

2変量正規乱数の生成方法をメモしておきます。
既に様々な言語・サイトで紹介されていることなので
何も真新しいことはないですが、備忘録として。

今回は2変量なのでDATAステップで実行してみます。
確率変数XとYがそれぞれ以下のように正規分布に
従うとします。
XとYの正規分布
このとき、X=xが得られた場合のYの条件付き分布は
以下になります。言うまでもありませんが、ρが0の
場合は上記Yの分布となりそれぞれ独立に乱数を生成
させた場合と同様になります。
2変量_条件付き正規分布
DATAステップでは、rand関数(正規分布はrand('normal',mean,sd))
でまずXの乱数を生成させ、そのあとに上記の条件付き分布に従った
Yの乱数を生成させます。下記の例では、平均値をそれぞれ-3.5,-2、
標準偏差を14,10、相関係数を0.6として2000組のデータを生成して、
corrプロシジャで各変数の要約統計量と相関係数を確認してみます。
*** 各パラメータ ;
%let _R = 0.6 ; *--- 相関係数 ;
%let _M1 = -3.5 ; *--- 変数1の母平均 ;
%let _M2 = -2 ; *--- 変数2の母平均 ;
%let _SD1 = 14 ; *--- 変数1のSD ;
%let _SD2 = 10 ; *--- 変数2のSD ;
*** 乱数の生成 ;
data _TEST ;
call streaminit(100) ; *--- seed ;
do I = 1 to 2000 ;
_V1 = rand('normal', &_M1, &_SD1) ; *--- 変数1 ;
_V2 = rand('normal', &_M2 + &_R*&_SD2*(_V1-&_M1)/&_SD1, &_SD2*sqrt(1-&_R**2)) ; *--- 変数2 ;
output ;
end ;
run ;
*** 平均値,標準偏差,相関係数を確認 ;
proc corr data=_TEST pearson ;
var _V: ;
run ;

結果は以下になります。相関係数が0.6付近であることが
確認できます。もちろん平均・標準偏差も意図通りの結果
となっています。
2変量_DATAステップの結果
今度は多変量の場合もまとめたいと思います。

Posted on 2014/02/22 Sat. 22:13 [edit]

CM: 0
TB: 0

0203

FREQプロシジャによる割合の差の解析  

FREQプロシジャを使用した割合の差に関する検定・推定方法を紹介します。
異なる薬剤を投与された集団における疾患の改善の割合(改善率)を比較
することを想定します。 割合の差の検定はtableステートメントのriskdiff
オプションで実行できます。

テストデータ:
変数GRP:投与群
変数VAL:カテゴリ(1が疾患の改善、2が非改善)
変数NN:各カテゴリの症例数

以下の解析を実行します。
解析方法FREQプロシジャのオプション
帰無仮説の割合に基づいた検定統計量riskdiff(equal column=1 var=null)
独立性のχ二乗検定chisq
Wald検定(標本割合に基づいた検定統計量)riskdiff(equal column=1 var=sample)
割合の差の正確な信頼区間exact riskdiff

プログラムは以下になります。
*** test data: two proportions ;
data xxx3 ;
input GRP VAL NN @@ ;
cards;
1 1 7 1 2 8 2 1 12 2 2 3
;
run ;
*** Compare two proportions ;
proc freq data=xxx3 ;
weight NN ;
*--- equal,var=null <=> chisq(conditional) ;
table GRP*VAL
/ riskdiff(equal column=1 var=null) chisq nocol nopercent ;
*--- equal,var=sample: Wald type ;
table GRP*VAL
/ riskdiff(equal column=1 var=sample) nocol nopercent ;
*--- exact test ;
exact riskdiff ;
run ;

結果は以下になります。①と②の検定結果が一致することが確認できます。
③のWald検定は統計の教科書に必ず出てくる方法です。④の割合の差の正確な
信頼区間は9.2から使用可能になりましたが、標本数が多い場合は計算の実行に
時間がかかります。
Riskdiff.png 

Posted on 2014/02/03 Mon. 18:56 [edit]

CM: 0
TB: 0

プロフィール

最新記事

最新コメント

最新トラックバック

月別アーカイブ

カテゴリ

訪れた人

▲Page top

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。