Documentation Example 14 for PROC MCMC
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: MCMCEX14 */
/* TITLE: Documentation Example 14 for PROC MCMC */
/* Time Independent Cox Model */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: MCMC */
/* DATA: */
/* */
/* SUPPORT: Fang Chen */
/* REF: PROC MCMC, EXAMPLE 14 */
/* MISC: */
/****************************************************************/
data Myeloma;
input Time Vstatus LogBUN HGB Platelet Age LogWBC Frac
LogPBM Protein SCalc;
label Time='survival time'
VStatus='0=alive 1=dead';
datalines;
1.25 1 2.2175 9.4 1 67 3.6628 1 1.9542 12 10
1.25 1 1.9395 12.0 1 38 3.9868 1 1.9542 20 18
2.00 1 1.5185 9.8 1 81 3.8751 1 2.0000 2 15
2.00 1 1.7482 11.3 0 75 3.8062 1 1.2553 0 12
2.00 1 1.3010 5.1 0 57 3.7243 1 2.0000 3 9
3.00 1 1.5441 6.7 1 46 4.4757 0 1.9345 12 10
5.00 1 2.2355 10.1 1 50 4.9542 1 1.6628 4 9
5.00 1 1.6812 6.5 1 74 3.7324 0 1.7324 5 9
6.00 1 1.3617 9.0 1 77 3.5441 0 1.4624 1 8
6.00 1 2.1139 10.2 0 70 3.5441 1 1.3617 1 8
6.00 1 1.1139 9.7 1 60 3.5185 1 1.3979 0 10
6.00 1 1.4150 10.4 1 67 3.9294 1 1.6902 0 8
7.00 1 1.9777 9.5 1 48 3.3617 1 1.5682 5 10
7.00 1 1.0414 5.1 0 61 3.7324 1 2.0000 1 10
7.00 1 1.1761 11.4 1 53 3.7243 1 1.5185 1 13
9.00 1 1.7243 8.2 1 55 3.7993 1 1.7404 0 12
11.00 1 1.1139 14.0 1 61 3.8808 1 1.2788 0 10
11.00 1 1.2304 12.0 1 43 3.7709 1 1.1761 1 9
11.00 1 1.3010 13.2 1 65 3.7993 1 1.8195 1 10
11.00 1 1.5682 7.5 1 70 3.8865 0 1.6721 0 12
11.00 1 1.0792 9.6 1 51 3.5051 1 1.9031 0 9
13.00 1 0.7782 5.5 0 60 3.5798 1 1.3979 2 10
14.00 1 1.3979 14.6 1 66 3.7243 1 1.2553 2 10
15.00 1 1.6021 10.6 1 70 3.6902 1 1.4314 0 11
16.00 1 1.3424 9.0 1 48 3.9345 1 2.0000 0 10
16.00 1 1.3222 8.8 1 62 3.6990 1 0.6990 17 10
17.00 1 1.2304 10.0 1 53 3.8808 1 1.4472 4 9
17.00 1 1.5911 11.2 1 68 3.4314 0 1.6128 1 10
18.00 1 1.4472 7.5 1 65 3.5682 0 0.9031 7 8
19.00 1 1.0792 14.4 1 51 3.9191 1 2.0000 6 15
19.00 1 1.2553 7.5 0 60 3.7924 1 1.9294 5 9
24.00 1 1.3010 14.6 1 56 4.0899 1 0.4771 0 9
25.00 1 1.0000 12.4 1 67 3.8195 1 1.6435 0 10
26.00 1 1.2304 11.2 1 49 3.6021 1 2.0000 27 11
32.00 1 1.3222 10.6 1 46 3.6990 1 1.6335 1 9
35.00 1 1.1139 7.0 0 48 3.6532 1 1.1761 4 10
37.00 1 1.6021 11.0 1 63 3.9542 0 1.2041 7 9
41.00 1 1.0000 10.2 1 69 3.4771 1 1.4771 6 10
41.00 1 1.1461 5.0 1 70 3.5185 1 1.3424 0 9
51.00 1 1.5682 7.7 0 74 3.4150 1 1.0414 4 13
52.00 1 1.0000 10.1 1 60 3.8573 1 1.6532 4 10
54.00 1 1.2553 9.0 1 49 3.7243 1 1.6990 2 10
58.00 1 1.2041 12.1 1 42 3.6990 1 1.5798 22 10
66.00 1 1.4472 6.6 1 59 3.7853 1 1.8195 0 9
67.00 1 1.3222 12.8 1 52 3.6435 1 1.0414 1 10
88.00 1 1.1761 10.6 1 47 3.5563 0 1.7559 21 9
89.00 1 1.3222 14.0 1 63 3.6532 1 1.6232 1 9
92.00 1 1.4314 11.0 1 58 4.0755 1 1.4150 4 11
4.00 0 1.9542 10.2 1 59 4.0453 0 0.7782 12 10
4.00 0 1.9243 10.0 1 49 3.9590 0 1.6232 0 13
7.00 0 1.1139 12.4 1 48 3.7993 1 1.8573 0 10
7.00 0 1.5315 10.2 1 81 3.5911 0 1.8808 0 11
8.00 0 1.0792 9.9 1 57 3.8325 1 1.6532 0 8
12.00 0 1.1461 11.6 1 46 3.6435 0 1.1461 0 7
11.00 0 1.6128 14.0 1 60 3.7324 1 1.8451 3 9
12.00 0 1.3979 8.8 1 66 3.8388 1 1.3617 0 9
13.00 0 1.6628 4.9 0 71 3.6435 0 1.7924 0 9
16.00 0 1.1461 13.0 1 55 3.8573 0 0.9031 0 9
19.00 0 1.3222 13.0 1 59 3.7709 1 2.0000 1 10
19.00 0 1.3222 10.8 1 69 3.8808 1 1.5185 0 10
28.00 0 1.2304 7.3 1 82 3.7482 1 1.6721 0 9
41.00 0 1.7559 12.8 1 72 3.7243 1 1.4472 1 9
53.00 0 1.1139 12.0 1 66 3.6128 1 2.0000 1 11
57.00 0 1.2553 12.5 1 66 3.9685 0 1.9542 0 11
77.00 0 1.0792 14.0 1 60 3.6812 0 0.9542 0 12
;
proc sort data = Myeloma;
by descending time;
run;
data _null_;
set Myeloma nobs=_n;
call symputx('N', _n);
stop;
run;
title 'Cox Model with Time Independent Covariates';
proc freq data=myeloma;
ods select none;
tables time / out=freqs;
run;
proc sort data = freqs;
by descending time;
run;
data myelomaM;
set myeloma;
ind = _N_;
run;
ods select all;
proc mcmc data=myelomaM outpost=outi nmc=50000 ntu=3000 seed=1;
ods select PostSumInt;
array beta[9];
parms beta: 0;
prior beta: ~ normal(0, var=1e6);
bZ = beta1 * LogBUN + beta2 * HGB + beta3 * Platelet
+ beta4 * Age + beta5 * LogWBC + beta6 * Frac +
beta7 * LogPBM + beta8 * Protein + beta9 * SCalc;
if ind = 1 then do; /* first observation */
S = exp(bZ);
l = vstatus * bZ;
v = vstatus;
end;
else if (1 < ind < &N) then do;
if (lag1(time) ne time) then do;
l = vstatus * bZ;
l = l - v * log(S); /* correct the loglike value */
v = vstatus; /* reset v count value */
S = S + exp(bZ);
end;
else do; /* still a tie */
l = vstatus * bZ;
S = S + exp(bZ);
v = v + vstatus; /* add # of nonsensored values */
end;
end;
else do; /* last observation */
if (lag1(time) ne time) then do;
l = - v * log(S); /* correct the loglike value */
S = S + exp(bZ);
l = l + vstatus * (bZ - log(S));
end;
else do;
S = S + exp(bZ);
l = vstatus * bZ - (v + vstatus) * log(S);
end;
end;
model general(l);
run;
data myelomaM;
merge myelomaM freqs(drop=percent);
by descending time;
retain stop;
if first.time then do;
stop = _n_ + count - 1;
end;
run;
data a;
run;
proc mcmc data=a outpost=outa nmc=50000 ntu=3000 seed=1 jointmodel;
ods select none;
array beta[9];
array data[1] / nosymbols;
array timeA[1] / nosymbols;
array vstatusA[1] / nosymbols;
array stopA[1] / nosymbols;
array bZ[&n];
array S[&n];
begincnst;
rc = read_array("myelomam", data, "logbun", "hgb", "platelet",
"age", "logwbc", "frac", "logpbm", "protein", "scalc");
rc = read_array("myelomam", timeA, "time");
rc = read_array("myelomam", vstatusA, "vstatus");
rc = read_array("myelomam", stopA, "stop");
endcnst;
parms (beta:) 0;
prior beta: ~ normal(0, var=1e6);
jl = 0;
/* calculate each bZ and cumulatively adding S as if there are no ties.*/
call mult(data, beta, bZ);
S[1] = exp(bZ[1]);
do i = 2 to &n;
S[i] = S[i-1] + exp(bZ[i]);
end;
do i = 1 to &n;
/* correct the S[i] term, when needed. */
if(stopA[i] > i) then do;
do j = (i+1) to stopA[i];
S[i] = S[i] + exp(bZ[j]);
end;
end;
jl = jl + vstatusA[i] * (bZ[i] - log(S[i]));
end;
model general(jl);
run;
ods select all;
proc compare data=outi compare=outa;
ods select comparesummary;
var beta1-beta9;
run;
proc phreg data=Myeloma;
ods select PostSumInt;
model Time*VStatus(0)=LogBUN HGB Platelet Age LogWBC
Frac LogPBM Protein Scalc;
bayes seed=1 nmc=10000 outpost=phout;
run;