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

徒然なるままにSAS暮らし

----

スポンサーサイト  

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

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

CM: --
TB: --

0228

call symputとcall symputx  

後輩と話をしていたときの小ネタです。call symputとcall symputxの違いは、
symputxが前後のスペースを削除した上でマクロ変数に値を格納することですが、
数値をそのままマクロ変数に格納しようとした場合、symputはログに数値を文字値
に変換したことが出力されますが、symputxではそれさえも出力されなくなります。
内部で変換についてどのような違いがあるのかはわかりませんが、とりあえず
スペースも消してくれるし、新しいほうのsymputxを使っておけば問題なさそう
といった感じでしょうか。
data TEST ;
X = 15 ;
call symput("XM",X) ;
run ;
%put *** &XM *** ;

data TEST ;
X = 15 ;
call symputx("XM",X) ;
run ;
%put *** &XM *** ;

ログ画面は以下になります。symputxではNOTEも出力されず、
かつスペースも削除されていることが確認されます。
symputx.png
スポンサーサイト

Posted on 2014/02/28 Fri. 00: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

0223

スタイルテンプレートをまとめて出力  

SASでは、ODSで外部ファイルに出力する際に多くのスタイルテンプレートが
用意されています。約50種類もあるので、私の場合、どれを使用するか、また、
各スタイルテンプレートがどんな見た目なのか毎回確かめつつ使用していました。
ただ、いつも出力して確認する時間が無駄なので、1つのPDFファイルに全スタイル
テンプレートの出力サンプルをまとめておいて使うようにしています。

プログラムの流れは、
①レポート用のテストデータ作成
②辞書テーブルからスタイルテンプレートの名前と数を取得(マクロ変数に格納)
 (dictionary.stylesにスタイルのリストが格納されています)
③取得したスタイルごとにODS PDFでレポート(表とグラフ)を作成
 (レポートはc:\tempに出力されます)
となります。
*** おまじない ;
options papersize=A4 orientation=landscape mprint mlogic nodate nonumber ;
*** テストデータ ;
data TEST ;
do ID = 1 to 100 ;
CITY = rantbl(100, 1/5, 1/5, 1/5, 1/5, 1/5) ;
SALES = round(1000000 + 100000*rannor(101)) ;
GRADE = rantbl(102,1/4,1/4,1/4,1/4) ;
output ;
end ;
run ;
*** スタイル名を格納している辞書テーブルからスタイル名と数をスペース区切りで取得 ;
proc sql noprint ;
select STYLE, count(STYLE) into : _STYLE separated by " ",
: _NN
from dictionary.styles where upcase(style) like 'STYLES%';

quit ;
*** スタイル名と数をログに出力 ;
%put スタイル名: ***** &_STYLE ***** ;
%put スタイル数: ***** &_NN ***** ;
*** 一つの表を全てのスタイルで出力するマクロ ;
%macro _OUT_DEF_STYLE ;
ods listing close ;
ods pdf file="C:\temp\DEF_STYLE.pdf" columns=2 ;
*** スタイルごとに1頁ずつテーブルとグラフを出力 ;
%do I = 1 %to &_NN ;
ods pdf startpage=yes style=%scan(&_STYLE,&I," ") ; *--- 各スタイル名を抽出 ;
title "Style: %substr(%scan(&_STYLE,&I,' '),8)" ;
proc tabulate data=TEST format=yen10. ;
class CITY GRADE ;
var SALES ;
table CITY, GRADE*SALES*(n*f=best. mean) ;
run ;
proc sgplot data=TEST ;
histogram SALES ;
density SALES / type=normal ;
run ;
%end ;
ods pdf close ;
ods listing ;
%mend _OUT_DEF_STYLE ;
*** 実行 ;
%_OUT_DEF_STYLE \;

辞書テーブルは以下のように変数Styleにスタイルテンプレート名を格納しています。
vstyles.png

ログに出力したテンプレート名と数は以下になります。
スペース区切りでスタイル名(51個)が羅列されています。
スタイル名とスタイル数

最後に、出来上がったPDFファイル1ページ目を見てみましょう。
Analysisという名前のスタイルテンプレートです。後ろに残りの
50個が続いて出力されています。
Def_style.png

センスの良いテンプレートがどんどん増えてくれるとカスタマイズする
手間も省けるので期待したいですが最近バージョンが上がってもあまり
追加されていないような気がします。。

Posted on 2014/02/23 Sun. 11:15 [edit]

CM: 2
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

0218

Fringeプロットと密度プロット(GTL)  

GTLによるFringeプロットと密度プロットの作成方法をメモしておきます。
GTLは複雑で私は到底暗記してすぐに記述とかはできないので備忘録的に。
layout latticeとlayout overlayの中で3個のグラフを重ねて描画します。
Fringeプロットは下部に出力されるので、offsetminで少しだけスペースを
与えておきます。
histogram: ヒストグラム
densityplot: 密度プロット
fringeplot: フリンジプロット
*** テストデータの作成 ;
data _TEST ;
do ID = 1 to 100 ;
VAL1 = 10 + 5*rannor(100) ;
VAL2 = 50 + 15*rannor(100) ;
output ;
end ;
run ;
*** テンプレートの作成 ;
proc template;
define statgraph mygraphs.distribution;
begingraph ;
layout lattice / columns=2 rows=1 ; *--- グラフを横に並べる ;
layout overlay / yaxisopts=(offsetmin=.05 label="Percent") ;
histogram VAL1 / scale=percent ;
densityplot VAL1 / normal( ) ;
fringeplot VAL1 ;

endlayout;
layout overlay / yaxisopts=(offsetmin=.05 label="Percent") ;
histogram VAL2 / scale=percent ;
densityplot VAL2 / normal( ) ;
fringeplot VAL2 ;

endlayout;
endlayout;
endgraph ;
end;
run;
*** グラフの描画 ;
title "Histogram, Density and Frnge plot" ;
proc sgrender data=_TEST template=mygraphs.distribution ;
format VAL1 VAL2 8. ;
run ;
title ;


結果を見てみましょう。グラフ下部にFringeプロットが出力されています。
Fringe.png

Posted on 2014/02/18 Tue. 00:29 [edit]

CM: 0
TB: 0

0215

XMLファイルの入出力  

SASによるXMLファイルの入出力方法を紹介します。
といっても使用する機能はlibnameやODSといった
単純なものです。xmlといえば、製薬業界ではCDISCの
Define.xmlを思い浮かべますが、最近SASデータセット
をxmlファイルに変換して提出する方向にFDAやCDISCが
動いているので後々この辺りの機能が注目を集めるかも
しれません。

libnameの出力方法:
libname ライブラリ名 xml "xmlファイルのパス" オプション ;

odsの出力方法:
ods xml file="xmlファイルのパス" オプション ;

以下のようにいくつかのファイルを作ってみます(⑤は読み込み)。
①データの値をそのまま各要素の値として出力
②データの値を属性として出力
③CDISC ODM形式(つまりDefine.xml)で出力(CDISCのデータじゃないのでログにエラーが出ます)
④ODSでprintプロシジャの結果を出力
⑤①で作成したxmlファイルの読み込み

おもしろいのは、③の方法でDefine.xmlもどきが作成できるのですが、
コードリストの出力はformatactive=yesで変数のフォーマットを
出力してくれます。CDISCに馴染みがない方やこれから勉強される方は
「Define.xmlって何ぞや?」という話になるかと思いますが、別の機会に
紹介できればと思います。要はデータの内容を所定のXML形式で米国の
FDAや日本のPMDA等の規制当局に提出が必要な代物です。ただ、ODMの
バージョンが古そうなので開発は置き去りなのかなとも感じております。


*** テストデータの作成 ;
proc format ;
value XF 1 = "A" 2 = "B" 3 = "C" ;
run ;
data _TEST ;
do ID = 1 to 5 ;
X=rantbl(100,0.5,0.25,0.25) ;
output ;
end ;
format X XF. ;
run ;
*** ①libname xmlエンジン(拡張子はxmlでファイル名とデータセット名は別) ;
libname _XMLOUT xml "./TEST.xml" xmldataform=element xmltype=generic ; *--- 設定はデフォルト ;
*** データ出力 ;
data _XMLOUT._XMLTEST ;
set _TEST ;
run ;
*** ②xmldataform=attribute: 属性としてデータを出力 ;
libname _XMLOUT2 xml "./TEST2.xml" xmldataform=attribute xmltype=generic ;
data _XMLOUT2._XMLTEST ;
set _TEST ;
run ;
*** ③CDISC ODM形式でも出力可能(formatactive=yesでコードリスト作成) ;
libname _XMLOUT3 xml "./TEST3.xml" xmltype=cdiscodm formatactive=yes ;
data _XMLOUT3.DM ;
set _TEST ;
run ;
*** ④ODS XMLでアウトプットをXMLファイルに出力 ;
*----- 注意 : ODSで吐き出したXMLファイルは普通には読み込めない ;
ods listing close ;
ods xml file="./_ODSXML.xml" ;
proc print data=_TEST ;
run ;
ods xml close ;
ods listing ;
*** ⑤XMLファイルの読み込み ;
*--- データ読み込み(xmltypeがgeneric(well-formed,つまり整形式)のものはそのまま読み込める) ;
data WK_XMLIN ;
set _XMLOUT._XMLTEST ;
run ;


順に結果を見てみましょう。SASのデータは容量が大きいので、
テキストであるXMLでデータを保存しておく場合等に有効かと
思います。

テストデータ
xml_テストデータ
①の結果
test_xml_1.png
②の結果
test_xml_2.png
③の結果(下のほうにCodelist関連のタグも出力されています)
test_xml_3.png
④の結果
ods_xml.png
⑤の結果
xml_テストデータ_in



Posted on 2014/02/15 Sat. 12:25 [edit]

CM: 2
TB: 0

0212

レポートで値によって動的にセルの背景色を変える方法  

データの値によってレポート内のセルの背景色を変える方法を紹介します。
既に当たり前のように使用されている機能ですが、備忘録的にメモを残して
おきます。有害事象の発現率や症例数設計の検出力、相関係数の高い値を
強調する等、様々な場面で活躍してくれる有用な機能です。

プログラムの流れは以下になります。
・値によって背景色を変えるフォーマットを作成する。
・レポート作成プロシジャのstyleステートメントで作成したフォーマットを適用する。

ここでは、副作用とその発現率のデータが得られたと仮定して、5%以上のセルの背景色を赤に、
5%未満は青色にするフォーマットを作成してREPORTプロシジャで適用しています。

*** 副作用データ ;
data _AEPCT ;
input AE $ PCT ;
cards;
感冒 10.2
下痢 8.3
頭痛 5.3
肺炎 0.5
骨折 0.2
;
run ;
*** 動的に色分けを行うためのフォーマット ;
proc format ;
value _BCGF low -< 5.0 = "blue"
5.0 - high = "red" ;
run ;

*** styleオプションで色分けのフォーマットを変数PCTに適用 ;
ods listing close ;
ods pdf ;
proc report data=_AEPCT nowd ;
column AE PCT ;
define AE / display ;
define PCT / display style(column)={background=_BCGF.} ;
run ;
ods pdf close ;
ods listing ;

結果は以下になります。データが更新されてもプログラムを変更する
必要もないので非常に重宝する機能です。
report_color.jpg

Posted on 2014/02/12 Wed. 23:09 [edit]

CM: 3
TB: 0

0211

逆累積分布曲線(Reverse Cumulative Distribution Plot)  

逆累積分布曲線(Reverse Cumulative Distribution(RCD)プロット)を作成
したのでメモしておきます。RCDプロットは、ワクチンの抗体価の分布を
比較する際などに用いられます。作り方は色々あるかと思いますが、
ここではFREQプロシジャで累積パーセントを算出した後にDATAステップ
で逆累積パーセントに変換し、SGプロシジャのSTEPステートメントで
グラフを作成します。
*** テストデータ(連続量で乱数生成) ;
data _TEST ;
do ID = 1 to 100 ;
GRP = rantbl(100,1/2,1/2) ;
VAL = 20 + 2*GRP + 5*rannor(100) ;
output ;
end ;
run ;
proc sort data=_TEST ; by GRP ; run ;
*** 累積パーセント算出 ;
ods listing close ;
ods output Onewayfreqs=_OUT1(keep=GRP VAL Cumpercent) ;
proc freq data=_TEST ;
by GRP ;
table VAL ;
run ;
ods output close ;
ods listing ;
*** 逆累積パーセント算出 ;
data _OUT1 ;
set _OUT1 ;
by GRP ;
_CDF = 100-Cumpercent ; *--- 逆累積パーセント ;
output ;
*--- 始点から一つ目のデータまでの描画用 ;
if last.GRP then do ;
_CDF = 100 ;
VAL = 0 ;
output ;
end ;

run ;
proc sort data=_OUT1 out=_OUT1_S ;
by descending _CDF ;
run ;
*** 各群のデータを転置 ;
proc transpose data=_OUT1_S out=_OUT1_ST(drop=_NAME_) prefix=_G ;
by descending _CDF ;
var VAL ;
id GRP ;
run ;
*** RCDプロット作成(justify=leftがデフォルト) ;
proc sgplot data=_OUT1_ST ;
step x=_G1 y=_CDF / justify=left name="Treat 1"
legendlabel="Treat 1" ;
step x=_G2 y=_CDF / justify=left name="Treat 2"
legendlabel="Treat 2" ;

keylegend "Treat 1" "Treat 2" ;
xaxis label="Test Value" values=(0 to 40 by 5) ;
yaxis label="Percent" values=(0 to 100 by 10) ;
run ;

・テストデータセット_TEST
_test.png

・描画用データセット_OUT1_ST
_rcd_dt.png

・逆累積分布曲線
RCD.png

LIFETESTプロシジャのカップランマイヤー図でも同じように描けそうですが、
カスタマイズが簡単にはできないのでしばらくこれを使おうかと思います。

Posted on 2014/02/11 Tue. 23:31 [edit]

CM: 2
TB: 0

0210

FCMPプロシジャによる生存時間解析の症例数設計  

FCMPプロシジャで生存時間解析の症例数設計を行います。
Schoenfeldの式と呼ばれる計算式で必要イベント数と症例数を
算出します。

想定する状況:
2つの投与群の生存関数を比較する際に、各群の対数ハザード
(ハザード:瞬間死亡率)が正規分布に従うと仮定しています。
λ1,λ2: 各群のハザード
S(t1), S(t2): 各群の時間tにおける生存率
d: 各群の必要イベント数
N: 各群の必要症例数
CALCS.png

関数の引数は、各群の生存率、時間、AlphaとPower(検出力)となります。
FCMPに「$」を指定することで、戻り値を文字変数とすることができます。
*** 生存関数の比較における各群の症例数算出(Schoenfeldの式) ;
proc fcmp outlib=WORK.FUNCDT.DIFFS ;
function _CALCS(S1,S2,T,Power,Alpha) $ 30 ; *--- 戻り値は文字 ;
H1 = (-1)*log(S1)/T ; *--- A群のハザード ;
H2 = (-1)*log(S2)/T ; *--- B群のハザード ;
HR = H2/H1 ; *--- ハザード比 ;
Za = probit(1-Alpha/2) ; *--- 上側(1-Alpha/2)*100%点 ;
Zb = probit(Power) ; *--- 上側(1-Beta)*100%点 ;
dSchoen = 2*(Za+Zb)**2/(log(HR)**2) ; *--- 各群の必要イベント数 ;
N_S = ceil(2*dSchoen/(2-S1-S2)) ; *--- 各群の必要症例数 ;
*--- 必要イベント数と必要症例数を文字変数で返す ;
D_N = "#Events: "||compress(put(dSchoen,6.))
||' N/group: '||compress(put(N_S,best.)) ;

return(D_N) ;
endsub ;
run ;
*** 関数の実行 ;
options cmplib=WORK.FUNCDT ;
data RES ;
length NN $30. ;
NN = _CALCS(0.6,0.7,4,0.05,0.8) ;
run ;
proc print data=RES ; run ;


結果を見てみましょう。各群の必要イベント数と必要症例数が文字列として
格納されています。
SS2.png


Posted on 2014/02/10 Mon. 23:42 [edit]

CM: 0
TB: 0

0210

これまでの外部発表  

これまでに外部で発表させて頂いた主な内容をまとめておきます。
元気があれば今後も頑張ります。

  • 「SAS ETL Serverを用いた臨床試験データの加工業務への活用事例」SAS Forumユーザー会学術総会2005

  • 「HTMLアプリケーションを用いた簡易SASツールの開発事例 -臨床試験における症例数設計用ツールの構築-」SASユーザー総会2008(共同発表者)

  • 「SGプロシジャとGTLによるグラフの作成とODS PDFによる統合解析帳票の作成~TQT試験における活用事例~」SASユーザー総会2011

  • 「SASとHTMLアプリケーションによるCDISC ADaM形式の解析用データセットを用いた有害事象の解析帳票・グラフ簡易作成ツールの開発事例」SASユーザー総会2012

  • 「ODS LAYOUTによる ワクチンの臨床開発におけるDashboardの作成」SASユーザー総会2013

  • Simple Tool for Creating ADaM Define.xml for Statisticians in Pharmaceutical Companies Using SAS and HTML Application with Excel Metadata File. CDISC Japan Interchange 2013

Posted on 2014/02/10 Mon. 00:12 [edit]

CM: 0
TB: 0

プロフィール

最新記事

最新コメント

最新トラックバック

月別アーカイブ

カテゴリ

訪れた人

▲Page top

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