MIXEDプロシージャー randomステートメント

クロスオーバー試験などのように、被験者が入れ子になっているとき。

random ID(SEQ);

のように、PROC MIXEDでも入れ子表記しなければならない。あれ？　そうだったっけ。と思ったが、そうだった。

変量効果分散分析のように期待平方和をSASに出力させたかった。PROC MIXEDでわざわざTYPE3平方和を指定して実行させてみて、このことを再確認。

最近、線型混合効果モデルと少し離れていたので忘れてしまったようだ。2月から、「MEM in SをRとSASで」みたいな会を企画しているので、復習。

分位点計算

SASによる分位点の計算について、標準的な作業であるが質問の頻度が高い。

ods listing close;
PROC UNIVARIATE data=VPC;
by TIME;
var IPRED;
output pctlpts=2.5 5 50 95 97.5 pctlpre=pct_ out=out;
RUN;

任意の分位点計算をするなら、PROC UNIVARIATEがよい。

90%予測区間の場合は、PROC MEANSでもいい。

ods output summary=out;
ods listing close;
PROC MEAN data=VPC MEAN MEDIAN P5 P90;
by TIME;
var IPRED;
RUN;

経過時間の続き

Longitudinalなデータで、ある時点との差のようにObservationをまたいで計算するテクニック。SASは原則として一度に一つのObservationを処理するから、複数のObservationにまたがった処理をするにはテクニックが必要となる。そのときに、retainを使う以外に、First.ID, Last.IDを用いたりLAG関数を用いることもできる。

DATA demo;
do ID = 1 to 10;
do Time = 0,2,4,8;
output;
end;
end;
RUN;

DATA demo2;
set demo;
by ID;
FID=First.ID;
LID=Last.ID;
INTERVAL = Time - LAG(Time);
if LAG(Time) = . then INTERVAL = 0;
if FID=1 then INTERVAL=0;
RUN;

DATA demo3;
set demo;
by ID;
if First.ID then TIME2 = TIME;
retain TIME2;
TAD =TIME-TIME2;
RUN;

こちらの方が簡単かもしれない。

前の時点からの経過時間を計算する

必要があって、最終投与後時間からの経過時間を計算した。普通は生データはこうなってないから、下に示す計算はあまりやらないと思うのだが。

初回投与後8時間で2度目の投与をし12時間、24時間に採血したとする。プログラムとしては各時点の間隔をINTERVALという変数に格納し、採血した時点についてはINTERVALを積算する。まず、基本的なテクニックは以下の通り。

DATA Example;
retain TIME_OLD;
TIME_OLD = 0;
do TIME = 0, 1, 2, 4, 8, 12, 24;
    INTERVAL = TIME - TIME_OLD;
    output;
  TIME_OLD = TIME;
end;
RUN;

DATA Example2;
set Example;
MDV = 1;
if TIME = 12 or TIME = 24 then MDV = 0;
RUN;

DATA Example3;
set Example2;
sum + INTERVAL ;
if MDV  =0 then sum2 + INTERVAL;
RUN;

こんなかんじ。retainｽﾃｰﾄﾒﾝﾄを使うあたりは慣れないと戸惑うかもしれない。

NONMEMのデータセットが既にできていて、そこから構成する場合は以下の通り。

DATA DAT_Plot;set DAT;RUN;

DATA DAT_Plot2;
set DAT_Plot;
if TIME = 0 then do;TIME_OLD=0;INTERVAL=0;end;
if DV ne 0 then MDV = 0;else MDV=1;
retain TIME_OLD;
INTERVAL = TIME-TIME_OLD;
output;
TIME_OLD = TIME;
keep ID TIME TIME_OLD INTERVAL DV MDV;
RUN;

DATA DAT_Plot3;
set DAT_Plot2;
if MDV = 0 then sum+INTERVAL;
else sum=0;
rename sum=TARD;
RUN;

PROC GPLOT data=DAT_Plot3;
where TARD ne 0;
plot DV * TARD;
symbol value=circle color=blue i=none;
RUN;QUIT;

まぁ、あまりやらないと思う。生データが初回投与後時間になっているということはあまりないし。あくまでNONMEMデータが先に手に入っているとか、それを作った後でプロットを作るというようなとき。

GLIMMIX

今日はモデルについての復習。PROC GENMOD, GLIMMIX, NLMIXEDといったものを整理したPPT資料を作った。
PROC GLIMMIX data=infection2  order=data OR;
class trt clinic;
model cure=trt / dist=bin link = logit solution cl;
random intercept /subject=clinic;
estimate "treatment" trt 1 -1 /exp;
RUN;
比較のためにNLMIXEDを使った。NLMIXEDは小さいデータなら実行速度が気にならない。QPOINTS=20 TECH=NEWRAPとか、重い設定と聞いていたのに。QPOINTSはGauss積分の分位点の数を決めるので当然多いと実行速度が遅くなるのだが。考えてみれば、Iterationの少ない数で収束できる問題については、気にならない速度になるのかもしれない。
あるWebで、「全てをNLMIXEDを使って解析してみよう。そうすれば各統計手法が何をしているのかがわかる」という記述があった。まぁ最近そうかなと感じる。というのもNLMIXEDは結局のところ、尤度を解析者がコードするようなものだし、カテゴリー変数のためのダミー変数も自分で作らないとだめだから、理解は深まるだろうなと。

小さい部品

日頃使う小さいプログラムに手を入れたりしようとしている。
DATA Cp_v;
set Cp_h;
length Cp $ 7;
array CC{13}  C0 C0_25 C0_5 C1 C1_5 C2 C4 C6 C8 C12 C24 C48 C72;
do i = 1 to 13;
Cp = CC{i};
output;
end;
keep id Cp;
RUN;

よくある処理なのだけれど、Cp_hはDDEを使ってExcelから拾ったりするときも無念ながらある。ましてやそのExcelがCSRのPDFから作成されることも。社外から来るデータだからであるが。

IML

ある目的があって、IMLの勉強を続ける。思ったよりも使いやすいかもしれない。SAS dataset に行列の内容を出力させるときに、append も必要ということに気づかず、しばらく迷う。

CWRESの計算ﾌﾟﾛｸﾞﾗﾑをSASでちょっと作ってみたくなった。この手の計算のポイントがどこにあるかというと、残差の分散共分散行列の計算にかかっている。基本的には分散の合成と誤差の伝播則とみなすことができる。

η、εについて線形近似したとして、残差の分散共分散行列の計算は、C=(Gη＋Hε)・t（Gη＋Hε）の期待値を求めることになる。これから、Gηη`G`+Hεε`H`となる。よって、GΩG`+HΣH`という例の式が得られる。IMLのよい練習になりそう。

ｸﾞﾗﾌをまたやる

ｸﾞﾗﾌで視覚的にﾃﾞｰﾀを見るということは統計解析においては非常に重要なﾌﾟﾛｾｽであるにも関わらず、統計解析そのものの数値計算とｸﾞﾗﾌ描画というのはｿﾌﾄにとってはかなり異なる機能である。そのため、各ｿﾌﾄともそれなりに開発に苦労しているさまがうかがえて、ｸﾞﾗﾌ機能が使いやすく優れているｿﾌﾄはやはり人気がある。

最近の趣味（？）は、GOPTIONSｽﾃｰﾄﾒﾝﾄなどのHELPなどみてできるだけきれいなｸﾞﾗﾌがかけるようにとしている。SAS/Graphのﾏﾆｭｱﾙとくびっぴき。今まで感覚で使っているところがわかってきた。GOPTIONS の vsize、hsize、vorigin、horiginなんて基礎機能についてけっこう誤解点があることに気づいて、あせる。

ｸﾞﾗﾌばかり

ｸﾞﾗﾌばかりいじっていた。annotate機能というものがおぼろげながらわかってくる。

annotate機能とはｸﾞﾗﾌに施す様々な修飾をSASﾃﾞｰﾀｾｯﾄにして取り込む機能だ。ふつうこういう風にｿﾌﾄを設計するだろうか？　とちょっと考えたが、これはこれで使えれば便利である。やっていることは簡単で、あるプロット領域にﾗﾍﾞﾙを打つ、線を引くといった作業であり、SASﾃﾞｰﾀｾｯﾄのなかにお約束の変数を指定して、そこに値を入れておくことでｸﾞﾗﾌにいろいろな飾りをつけることができる。

まず初心者である自分の場合、お約束を覚えることが重要である。SASのｸﾞﾗﾌの座標設定は、プロットｴﾘｱ内でＸやＹ変数そのものの値で指定するやり方と、全体を通して振っている座標とがあり、これらを必要に応じて使い分ける。legendの指定にもそういう雰囲気がでていたと思うが、凡例とかﾀｲﾄﾙとかはoverallの座標がわかりやすい一方、あるﾃﾞｰﾀのﾏｰｸなどプロットエリアに打つことが前提のものはプロットｴﾘｱ座標がよい。それを、X、Y、horizontalの3つの方向についてそれぞれ指定できるというのは、なかなかﾏﾆｱｯｸだと思った。また、絶対値とパーセントでの相対指定もある。

ちょっとまじめに勉強もしないと。

ということで、ちょっとだけ数理の本を読むことにする。

今週はいろいろ、研究とは関係ないことで忙しくなる予定。

練習する

ｸﾞﾗﾌ作成の annotate 機能を触っていなかったので、使い方を調べることにする。SASは作りこみさえすればｸﾞﾗﾌはきれいなはず。