Documentation Example 13 for PROC MCMC
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: MCMCEX13 */
/* TITLE: Documentation Example 13 for PROC MCMC */
/* Exponential and Weibull Survival Analysis */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: MCMC */
/* DATA: */
/* */
/* SUPPORT: Fang Chen */
/* REF: PROC MCMC, EXAMPLE 13 */
/* MISC: */
/****************************************************************/
data e1684;
input t t_cen treatment @@;
if t = . then do;
t = t_cen;
v = 0;
end;
else
v = 1;
ifn = treatment - 1;
et = exp(t);
lt = log(t);
drop t_cen;
datalines;
1.57808 0.00000 2 1.48219 0.00000 2 . 7.33425 1
2.23288 0.00000 1 . 9.38356 2 3.27671 0.00000 1
. 9.64384 1 1.66575 0.00000 2 0.94247 0.00000 1
1.68767 0.00000 2 2.34247 0.00000 2 0.89863 0.00000 1
. 9.03288 2 . 9.63014 2 0.52603 0.00000 1
1.82192 0.00000 2 0.93425 0.00000 1 . 8.98630 2
3.35068 0.00000 1 8.67397 0.00000 1 0.41096 0.00000 2
2.78630 0.00000 1 2.56438 0.00000 1 . 8.75342 2
0.56986 0.00000 2 . 8.40000 1 . 7.25205 1
4.38630 0.00000 2 . 8.36712 2 . 8.99178 2
0.86575 0.00000 2 . 4.76986 1 1.15616 0.00000 2
. 7.28767 1 3.13151 0.00000 1 . 8.55068 2
. 8.45753 2 4.59452 0.00000 1 2.88219 0.00000 2
0.89589 0.00000 1 1.76164 0.00000 2 . 7.81370 1
. 8.33425 2 2.62192 0.00000 1 0.16164 0.00000 2
. 8.24658 1 1.52603 0.00000 1 5.30959 0.00000 1
0.87123 0.00000 2 0.41644 0.00000 1 4.24110 0.00000 1
0.13699 0.00000 1 7.07671 0.00000 2 0.13151 0.00000 2
. 8.02740 1 . 6.16164 2 1.29863 0.00000 2
1.29041 0.00000 2 . 7.99726 1 . 8.34795 1
. 7.30137 2 2.32877 0.00000 2 0.56438 0.00000 1
5.62740 0.00000 2 1.23014 0.00000 1 . 7.94521 1
5.06301 0.00000 1 3.27671 0.00000 2 . 0.60822 2
0.65753 0.00000 1 0.84110 0.00000 2 . 8.40000 2
0.18356 0.00000 1 2.62466 0.00000 2 . 7.96438 2
. 7.77808 1 0.22192 0.00000 1 2.33973 0.00000 1
0.52329 0.00000 1 . 8.04110 2 . 7.83288 1
0.64110 0.00000 1 0.38356 0.00000 1 . 7.82192 2
0.51781 0.00000 2 . 8.09863 2 . 8.16712 2
4.42740 0.00000 1 0.88493 0.00000 1 2.78356 0.00000 1
2.64658 0.00000 2 . 8.21370 2 . 7.41918 2
0.99726 0.00000 1 5.88493 0.00000 2 0.41644 0.00000 1
3.53699 0.00000 1 . 7.56164 1 . 7.53151 1
0.27671 0.00000 1 0.76986 0.00000 2 . 7.62192 2
. 7.79726 1 0.64110 0.00000 1 1.14521 0.00000 2
2.01644 0.00000 1 2.84384 0.00000 1 . 7.00000 2
1.27397 0.00000 2 . 7.09589 1 2.04110 0.00000 1
0.83562 0.00000 1 0.92329 0.00000 1 0.07397 0.00000 1
. 7.30685 2 2.07671 0.00000 2 . 7.70959 2
. 6.15890 1 . 6.89315 2 3.30685 0.00000 1
0.36164 0.00000 1 1.97808 0.00000 2 1.23836 0.00000 2
0.10685 0.00000 1 . 7.63836 1 2.06301 0.00000 1
. 7.42466 2 0.50959 0.00000 1 0.65753 0.00000 1
. 6.93151 1 . 7.23288 2 6.01096 0.00000 1
0.33699 0.00000 1 . 6.47123 2 0.94795 0.00000 1
2.91781 0.00000 2 1.59726 0.00000 2 0.84932 0.00000 2
1.38356 0.00000 1 3.81644 0.00000 2 . 7.06849 1
. 7.04110 2 1.00274 0.00000 2 . 6.34795 2
1.18082 0.00000 1 0.97534 0.00000 1 2.16712 0.00000 1
. 6.85479 1 1.38356 0.00000 1 1.71507 0.00000 2
0.79452 0.00000 2 . 6.86301 2 . 6.50411 1
0.42466 0.00000 2 0.98630 0.00000 1 . 6.13699 2
3.80000 0.00000 2 . 6.48493 1 . 6.96438 2
. 6.78082 2 0.56164 0.00000 2 2.67123 0.00000 1
1.56712 0.00000 2 2.07397 0.00000 2 0.33973 0.00000 1
3.37808 0.00000 2 3.15068 0.00000 1 . 6.81096 2
3.20822 0.00000 2 0.62740 0.00000 1 1.64384 0.00000 1
1.40822 0.00000 1 . 6.06575 1 1.66301 0.00000 2
1.36986 0.00000 2 5.46849 0.00000 1 0.42740 0.00000 1
1.13973 0.00000 2 1.73699 0.00000 2 . 5.54521 2
0.85205 0.00000 1 0.43014 0.00000 1 1.20822 0.00000 2
4.36164 0.00000 1 0.52877 0.00000 2 . 6.51507 1
2.89863 0.00000 2 . 6.20274 2 1.21644 0.00000 2
. 6.00000 2 . 6.25479 1 . 6.49863 1
1.13699 0.00000 2 1.69589 0.00000 1 . 6.41096 2
. 6.02192 1 3.04932 0.00000 2 . 5.62740 2
0.72603 0.00000 1 0.73425 0.00000 2 1.47945 0.00000 2
0.37808 0.00000 2 . 5.75890 2 1.48219 0.00000 2
. 5.88493 1 . 1.80274 1 1.40548 0.00000 2
. 4.74795 1 . 5.24658 1 0.29041 0.00000 1
. 5.83836 1 . 5.32055 2 5.16712 0.00000 2
. 5.59178 2 . 5.77808 1 0.53425 0.00000 2
. 2.22466 1 3.59726 0.00000 1 . 5.32329 1
1.78630 0.00000 2 0.70411 0.00000 2 . 4.94795 2
. 5.45479 2 4.32877 0.00000 1 1.16164 0.00000 2
. 5.20274 2 . 4.40822 1 1.41096 0.00000 1
. 4.92877 2 . 5.42192 2 0.98904 0.00000 1
0.36438 0.00000 1 . 4.38082 1 0.77260 0.00000 2
4.90959 0.00000 2 1.26849 0.00000 1 0.58082 0.00000 2
. 4.95616 1 . 5.12329 1 . 4.74795 1
. 4.90685 2 1.41918 0.00000 1 0.44110 0.00000 1
. 4.29863 2 . 4.63836 2 . 4.81370 1
4.50137 0.00000 2 3.92329 0.00000 2 . 4.86027 2
0.52603 0.00000 1 2.10685 0.00000 2 . 4.24384 1
3.39178 0.00000 1 . 4.36164 2 . 4.81918 2
;
/****************************************************************/
/* Exponential Survival Model */
/****************************************************************/
title 'Exponential Survival Model';
ods graphics on;
proc mcmc data=e1684 outpost=expsurvout nmc=10000 seed=4861
diag=(mcse ess);
ods select PostSumInt TADpanel
ess mcse;
parms (beta0 beta1) 0;
prior beta: ~ normal(0, sd = 10000);
/*****************************************************/
/* (1) the logpdf and logsdf functions are not used */
/*****************************************************/
/* nu = 1/exp(beta0 + beta1*ifn);
llike = v*logpdf("exponential", t, nu) +
(1-v)*logsdf("exponential", t, nu);
*/
/****************************************************/
/* (2) the simplified likelihood formula is used */
/****************************************************/
l_h = beta0 + beta1*ifn;
llike = v*(l_h) - t*exp(l_h);
model general(llike);
run;
/****************************************************************/
/* Weibull Survival Model */
/****************************************************************/
title 'Weibull Survival Model';
proc mcmc data=e1684 outpost=weisurvout nmc=10000 seed=1234
monitor=(_parms_ surv_ifn surv_noifn) stats=(summary intervals);
ods select PostSummaries;
ods output PostSummaries=ds PostIntervals=is;
array surv_ifn[10];
array surv_noifn[10];
parms alpha 1 (beta0 beta1) 0;
prior beta: ~ normal(0, var=10000);
prior alpha ~ gamma(0.001,is=0.001);
beginnodata;
do t1 = 1 to 10;
surv_ifn[t1] = exp(-exp(beta0+beta1)*t1**alpha);
surv_noifn[t1] = exp(-exp(beta0)*t1**alpha);
end;
endnodata;
lambda = beta0 + beta1*ifn;
/*****************************************************/
/* (1) the logpdf and logsdf functions are not used */
/*****************************************************/
/* gamma = exp(-lambda /alpha);
llike = v*logpdf('weibull', t, alpha, gamma) +
(1-v)*logsdf('weibull', t, alpha, gamma);
*/
/****************************************************/
/* (2) the simplified likelihood formula is used */
/****************************************************/
llike = v*(log(alpha) + (alpha-1)*log(t) + lambda) -
exp(lambda)*(t**alpha);
model general(llike);
run;
proc format;
value alphafmt low-<1 = 'alpha < 1' 1-high = 'alpha >= 1';
run;
proc freq data=weisurvout;
tables alpha /nocum;
format alpha alphafmt.;
run;
/* define macro stackdata */
%macro StackData(dataset,output,vars);
data &output;
length var $ 32;
if 0 then set &dataset nobs=nnn;
array lll[*] &vars;
do jjj=1 to dim(lll);
do iii=1 to nnn;
set &dataset point=iii;
value = lll[jjj];
call vname(lll[jjj],var);
output;
end;
end;
stop;
keep var value;
run;
%mend;
/* stack the surv_ifn variables and saved them to survifn. */
%StackData(weisurvout, survifn, surv_ifn1-surv_ifn10);
proc sgplot data=survifn;
yaxis label='Survival Probability' values=(0 to 1 by 0.2);
xaxis label='Time' discreteorder=data;
vbox value / category=var;
run;
data surv;
set ds;
if _n_ >= 4 then do;
set is point=_n_;
group = 'with interferon ';
time = _n_ - 3;
if time > 10 then do;
time = time - 10;
group = 'without interferon';
end;
output;
end;
keep time group mean hpdlower hpdupper;
run;
proc sgplot data=surv;
yaxis label="Survival Probability" values=(0 to 1 by 0.2);
series x=time y=mean / group = group name='i';
band x=time lower=hpdlower upper=hpdupper / group = group transparency=0.7;
keylegend 'i';
run;
ods graphics off;
/****************************************************************/
/* Model Comparison between Weibull and Exponential */
/****************************************************************/
title 'Model Comparison between Weibull and Exponential';
proc mcmc data=e1684 outpost=weisurvout nmc=10000 seed=4861 dic;
ods select dic;
parms alpha 1 (beta0 beta1) 0;
prior beta: ~ normal(0, var=10000);
prior alpha ~ gamma(0.001,is=0.001);
lambda = beta0 + beta1*ifn;
llike = v*(log(alpha) + (alpha-1)*log(t) + lambda) -
exp(lambda)*(t**alpha);
model general(llike);
run;
proc mcmc data=e1684 outpost=expsurvout nmc=10000 seed=4861 dic;
ods select dic;
parms beta0 beta1 0;
prior beta: ~ normal(0, var=10000);
begincnst;
alpha = 1;
endcnst;
lambda = beta0 + beta1*ifn;
llike = v*(log(alpha) + (alpha-1)*log(t) + lambda) -
exp(lambda)*(t**alpha);
model general(llike);
run;
proc mcmc data=e1684 outpost=expsurvout1 nmc=10000 seed=4861 dic;
ods select none;
parms (beta0 beta1) 0;
prior beta: ~ normal(0, sd = 10000);
l_h = beta0 + beta1*ifn;
llike = v*(l_h) - t*exp(l_h);
model general(llike);
run;
proc compare data=expsurvout compare=expsurvout1;
var beta0 beta1;
run;
