Documentation Example 16 for PROC MCMC
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: MCMCEX16 */
/* TITLE: Documentation Example 16 for PROC MCMC */
/* Piecewise exponential frailty models */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: MCMC */
/* DATA: */
/* */
/* SUPPORT: Fang Chen */
/* REF: PROC MCMC, EXAMPLE 16 */
/* MISC: */
/****************************************************************/
title 'Piecewise Exponential Model';
data Blind;
input ID Time Status DiabeticType Treatment @@;
datalines;
5 46.23 0 1 1 5 46.23 0 1 0 14 42.50 0 0 1 14 31.30 1 0 0
16 42.27 0 0 1 16 42.27 0 0 0 25 20.60 0 0 1 25 20.60 0 0 0
29 38.77 0 0 1 29 0.30 1 0 0 46 65.23 0 0 1 46 54.27 1 0 0
49 63.50 0 0 1 49 10.80 1 0 0 56 23.17 0 0 1 56 23.17 0 0 0
61 1.47 0 0 1 61 1.47 0 0 0 71 58.07 0 1 1 71 13.83 1 1 0
100 46.43 1 1 1 100 48.53 0 1 0 112 44.40 0 1 1 112 7.90 1 1 0
120 39.57 0 1 1 120 39.57 0 1 0 127 30.83 1 1 1 127 38.57 1 1 0
133 66.27 0 1 1 133 14.10 1 1 0 150 20.17 1 0 1 150 6.90 1 0 0
167 58.43 0 1 1 167 41.40 1 1 0 176 58.20 0 0 1 176 58.20 0 0 0
185 57.43 0 1 1 185 57.43 0 1 0 190 56.03 0 0 1 190 56.03 0 0 0
202 67.53 0 0 1 202 67.53 0 0 0 214 61.40 0 1 1 214 0.60 1 1 0
220 10.27 1 0 1 220 1.63 1 0 0 243 66.20 0 0 1 243 66.20 0 0 0
255 5.67 1 0 1 255 13.83 1 0 0 264 58.83 0 0 1 264 29.97 1 0 0
266 60.27 0 1 1 266 26.37 1 1 0 284 5.77 1 1 1 284 1.33 1 1 0
295 5.90 1 0 1 295 35.53 1 0 0 300 25.63 1 1 1 300 21.90 1 1 0
302 33.90 1 0 1 302 14.80 1 0 0 315 1.73 1 0 1 315 6.20 1 0 0
324 46.90 0 1 1 324 22.00 1 1 0 328 31.13 0 0 1 328 31.13 0 0 0
335 30.20 1 0 1 335 22.00 1 0 0 342 70.90 0 0 1 342 70.90 0 0 0
349 25.80 1 1 1 349 13.87 1 1 0 357 5.73 1 1 1 357 48.30 1 1 0
368 53.43 0 0 1 368 53.43 0 0 0 385 1.90 1 0 1 385 51.10 0 0 0
396 9.90 1 1 1 396 9.90 1 1 0 405 34.20 0 0 1 405 34.20 0 0 0
409 46.73 0 1 1 409 2.67 1 1 0 419 18.73 0 1 1 419 13.83 1 1 0
429 32.03 0 1 1 429 4.27 1 1 0 433 69.87 0 1 1 433 13.90 1 1 0
445 66.80 0 0 1 445 66.80 0 0 0 454 64.73 0 0 1 454 64.73 0 0 0
468 1.70 1 0 1 468 1.70 1 0 0 480 1.77 1 0 1 480 43.03 1 0 0
485 29.03 0 0 1 485 29.03 0 0 0 491 56.57 0 1 1 491 56.57 0 1 0
503 8.30 1 1 1 503 8.30 1 1 0 515 21.57 0 1 1 515 18.43 1 1 0
522 31.57 0 0 1 522 31.57 0 0 0 538 31.63 0 1 1 538 31.63 1 1 0
547 39.77 0 1 1 547 39.77 0 1 0 550 18.70 1 0 1 550 6.53 1 0 0
554 18.90 0 0 1 554 18.90 0 0 0 557 56.80 0 0 1 557 22.23 1 0 0
561 55.60 0 0 1 561 14.00 1 0 0 568 42.17 1 0 1 568 42.17 1 0 0
572 10.70 0 0 1 572 5.33 1 0 0 576 66.33 0 0 1 576 59.80 1 0 0
581 52.33 0 1 1 581 5.83 1 1 0 606 58.17 0 0 1 606 2.17 1 0 0
610 14.30 1 0 1 610 48.43 1 0 0 615 25.83 0 0 1 615 25.83 0 0 0
618 45.40 0 0 1 618 45.40 0 0 0 624 47.60 0 0 1 624 47.60 0 0 0
631 13.33 1 0 1 631 9.60 1 0 0 636 42.10 0 0 1 636 42.10 0 0 0
645 39.93 0 0 1 645 39.93 0 0 0 653 14.27 1 0 1 653 7.60 1 0 0
662 34.57 1 0 1 662 1.80 1 0 0 664 65.80 0 0 1 664 4.30 1 0 0
683 4.10 1 1 1 683 12.20 1 1 0 687 60.93 0 0 1 687 60.93 0 0 0
701 57.20 0 0 1 701 57.20 0 0 0 706 38.07 0 1 1 706 12.73 1 1 0
717 54.10 0 1 1 717 54.10 1 1 0 722 59.27 0 0 1 722 9.40 1 0 0
731 21.57 1 0 1 731 9.90 1 0 0 740 54.10 0 0 1 740 54.10 0 0 0
749 50.47 0 1 1 749 50.47 0 1 0 757 46.17 0 0 1 757 46.17 0 0 0
760 46.30 0 0 1 760 46.30 0 0 0 766 38.83 0 1 1 766 38.83 0 1 0
769 44.60 0 0 1 769 44.60 0 0 0 772 43.07 0 0 1 772 43.07 0 0 0
778 26.23 1 1 1 778 40.03 0 1 0 780 41.60 0 0 1 780 18.03 1 0 0
793 38.07 0 1 1 793 38.07 0 1 0 800 65.23 0 1 1 800 65.23 0 1 0
804 7.07 1 1 1 804 66.77 0 1 0 810 13.77 1 0 1 810 13.77 1 0 0
815 9.63 0 1 1 815 9.63 1 1 0 832 46.23 0 0 1 832 46.23 0 0 0
834 45.73 0 0 1 834 1.50 1 0 0 838 33.63 1 1 1 838 33.63 1 1 0
857 40.17 0 0 1 857 40.17 0 0 0 866 63.33 1 1 1 866 27.60 1 1 0
887 38.47 1 1 1 887 1.63 1 1 0 903 55.23 0 1 1 903 55.23 0 1 0
910 52.77 0 1 1 910 25.30 1 1 0 920 57.17 0 0 1 920 46.20 1 0 0
925 9.87 0 1 1 925 1.70 1 1 0 931 57.90 0 0 1 931 57.90 0 0 0
936 5.90 0 0 1 936 5.90 0 0 0 945 32.20 0 0 1 945 32.20 0 0 0
949 10.33 1 0 1 949 0.83 1 0 0 952 6.13 1 0 1 952 50.90 0 0 0
962 43.67 0 0 1 962 25.93 1 0 0 964 38.30 0 0 1 964 38.30 0 0 0
971 38.77 0 1 1 971 19.40 1 1 0 978 38.07 0 0 1 978 21.97 1 0 0
983 38.30 0 0 1 983 38.30 0 0 0 987 26.20 1 0 1 987 70.03 0 0 0
1002 62.57 0 0 1 1002 18.03 1 0 0 1017 13.83 1 1 1 1017 1.57 1 1 0
1029 46.50 0 1 1 1029 13.37 1 1 0 1034 11.07 1 0 1 1034 1.97 1 0 0
1037 42.47 0 1 1 1037 22.20 1 1 0 1042 38.73 0 1 1 1042 38.73 0 1 0
1069 51.13 0 1 1 1069 51.13 0 1 0 1074 6.10 1 0 1 1074 46.50 0 0 0
1098 2.10 1 0 1 1098 11.30 1 0 0 1102 17.73 1 0 1 1102 42.30 0 0 0
1112 26.47 0 0 1 1112 26.47 0 0 0 1117 10.77 0 0 1 1117 10.77 0 0 0
1126 55.33 0 1 1 1126 55.33 0 1 0 1135 58.67 0 0 1 1135 58.67 0 0 0
1145 12.93 1 1 1 1145 4.97 1 1 0 1148 54.20 0 1 1 1148 26.47 1 1 0
1167 49.57 0 0 1 1167 49.57 0 0 0 1184 24.43 1 1 1 1184 9.87 1 1 0
1191 50.23 0 1 1 1191 50.23 0 1 0 1205 13.97 1 0 1 1205 30.40 1 0 0
1213 43.33 0 0 1 1213 43.33 1 0 0 1228 42.23 0 1 1 1228 42.23 0 1 0
1247 74.93 0 0 1 1247 74.93 0 0 0 1250 66.93 0 1 1 1250 66.93 0 1 0
1253 73.43 0 0 1 1253 73.43 0 0 0 1267 67.47 0 1 1 1267 38.57 1 1 0
1281 3.67 0 1 1 1281 3.67 1 1 0 1287 48.87 1 0 1 1287 67.03 0 0 0
1293 65.60 0 0 1 1293 65.60 0 0 0 1296 15.83 0 0 1 1296 15.83 1 0 0
1309 20.07 0 1 1 1309 8.83 1 1 0 1312 67.43 0 0 1 1312 67.43 0 0 0
1317 1.47 0 0 1 1317 1.47 0 0 0 1321 62.93 0 0 1 1321 22.13 1 0 0
1333 6.30 1 0 1 1333 56.97 0 0 0 1347 59.70 0 0 1 1347 18.93 1 0 0
1361 13.80 1 0 1 1361 19.00 1 0 0 1366 55.13 0 1 1 1366 55.13 0 1 0
1373 13.57 1 0 1 1373 5.43 1 0 0 1397 42.20 0 1 1 1397 42.20 0 1 0
1410 38.27 0 1 1 1410 38.27 0 1 0 1413 7.10 0 0 1 1413 7.10 1 0 0
1425 63.63 0 1 1 1425 26.17 1 1 0 1447 59.00 0 0 1 1447 24.73 1 0 0
1461 54.37 0 1 1 1461 54.37 0 1 0 1469 54.60 0 1 1 1469 10.97 1 1 0
1480 63.87 0 1 1 1480 21.10 1 1 0 1487 62.37 0 1 1 1487 43.70 1 1 0
1491 62.80 0 1 1 1491 62.80 0 1 0 1499 63.33 0 1 1 1499 14.37 1 1 0
1503 58.53 0 1 1 1503 58.53 0 1 0 1513 58.07 0 1 1 1513 58.07 0 1 0
1524 58.50 0 1 1 1524 58.50 0 1 0 1533 1.50 1 1 1 1533 14.37 0 1 0
1537 54.73 0 0 1 1537 38.40 1 0 0 1552 50.63 0 0 1 1552 2.83 1 0 0
1554 51.10 0 1 1 1554 51.10 0 1 0 1562 49.93 0 1 1 1562 6.57 1 1 0
1572 46.27 0 1 1 1572 46.27 1 1 0 1581 10.60 0 1 1 1581 10.60 0 1 0
1585 42.77 0 1 1 1585 42.77 0 1 0 1596 34.37 1 0 1 1596 42.27 0 0 0
1600 42.07 0 0 1 1600 42.07 0 0 0 1603 38.77 0 0 1 1603 38.77 0 0 0
1619 74.97 0 1 1 1619 61.83 1 1 0 1627 6.57 1 0 1 1627 66.97 0 0 0
1636 38.87 1 0 1 1636 68.30 0 0 0 1640 42.43 1 0 1 1640 46.63 1 0 0
1643 67.07 0 0 1 1643 67.07 0 0 0 1649 2.70 1 0 1 1649 2.70 0 0 0
1666 63.80 0 0 1 1666 63.80 0 0 0 1672 32.63 0 0 1 1672 32.63 0 0 0
1683 62.00 0 1 1 1683 62.00 0 1 0 1688 13.10 1 0 1 1688 54.80 0 0 0
1705 8.00 0 0 1 1705 8.00 0 0 0 1717 51.60 0 1 1 1717 42.33 1 1 0
1727 49.97 0 1 1 1727 2.90 1 1 0 1746 45.90 0 0 1 1746 1.43 1 0 0
1749 41.93 0 1 1 1749 41.93 0 1 0
;
data partition;
input int_1-int_9;
datalines;
0.1 6.545 13.95 26.47 38.8 45.88 54.35 62 80
;
%let n = 8;
data _a;
set blind;
if _n_ eq 1 then set partition;
array int[*] int_:;
array Y[&n];
array dN[&n];
do k = 1 to (dim(int)-1);
Y[k] = (time - int[k] + 0.001 >= 0);
dN[k] = Y[k] * ( int[k+1] - time - 0.001 >= 0) * status;
end;
output;
drop int_: k;
run;
proc print data=_a(obs=10);
run;
data _b;
set _a;
array y[*] y:;
array dn[*] dn:;
do i = 1 to (dim(y));
y_val = y[i];
dn_val = dn[i];
int_index = i;
output;
end;
keep y_: dn_: diabetictype treatment int_index id;
run;
data _b;
set _b;
rename y_val=Y dn_val=dN;
run;
proc print data=_b(obs=10);
run;
data inputdata;
set _b;
if Y > 0;
run;
proc mcmc data=inputdata nmc=10000 outpost=postout seed=12351
maxtune=5;
ods select PostSumInt ESS;
parms beta1-beta3 0;
prior beta: ~ normal(0, var = 1e6);
random lambda ~ gamma(0.01, iscale = 0.01) subject=int_index;
bZ = beta1*treatment + beta2*diabetictype + beta3*treatment*diabetictype;
idt = exp(bz) * lambda;
model dN ~ poisson(idt);
run;
ods graphics on;
%cater(data=postout, var=lambda_:);
ods graphics off;
ods select none;
proc mcmc data=inputdata nmc=10000 outpost=postout seed=12351
stats=summary diag=none;
parms beta1-beta3 0 s2;
prior beta: ~ normal(0, var = 1e6);
prior s2 ~ igamma(0.01, scale=0.01);
random lambda ~ gamma(0.01, iscale = 0.01) subject=int_index;
random u ~ normal(0, var=s2) subject=id;
bZ = beta1*treatment + beta2*diabetictype + beta3*treatment*diabetictype + u;
idt = exp(bZ) * lambda;
model dN ~ poisson(idt);
run;