Time Independent Model

The regression model considered in this example uses the following formula:

The hard part of coding this in PROC MCMC is the construction of the risk set . contains all observations that have survival time greater than or equal to . First suppose that there are no ties in time. Sorting the data set by the variable time into descending order gives you that is in the right order. Observation ’s risk set consists of all data points such that in the data set. You can cumulatively increment S in the SAS statements.

With potential ties in time, at observation , you need to know whether any subsequent observations, and so on, have the same survival time as . Suppose that the th, the th, and the th observations all have the same survival time; all three of them need to be included in the risk set calculation. This means that to calculate the likelihood for some observations, you need to access both the previous and subsequent observations in the data set. There are two ways to do this. One is to use the LAG function; the other is to use the option JOINTMODEL.

The LAG function returns values from a queue (see SAS Language Reference: Dictionary). So for the th observation, you can use LAG1 to access variables from the previous row in the data set. You want to compare the lag1 value of time with the current time value. Depending on whether the two time values are equal, you can add correction terms in the calculation for the risk set S.

The following statements sort the data set by time into descending order, with the largest survival time on top:

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;

The following statements run PROC MCMC and produce Output 52.8.1:

proc mcmc data=myelomaM outpost=outi nmc=50000 ntu=3000 seed=1;
   ods select PostSummaries PostIntervals;
   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;

The symbol bZ is the regression term, which is independent of the time variable. The symbol ind indexes observation numbers in the data set. The symbol S keeps track of the risk set term for every observation. The symbol l calculates the log likelihood for each observation. Note that the value of l for observation ind is not necessarily the correct log likelihood value for that observation, especially in cases where the observation ind is in the tied times. Correction terms are added to subsequent values of l when the time variable becomes different in order to make up the difference. The total sum of l calculated over the entire data set is correct. The symbol v keeps track of the sum of vstatus, as censored data do not enter the likelihood and need to be taken out.

You use the function LAG1 to detect if two adjacent time values are different. If they are, you know that the current observation is in a different risk set than the last one. You then need to add a correction term to the log likelihood value of l. The IF-ELSE statements break the observations into three parts: the first observation, the last observation and everything in the middle.

An alternative to using the LAG function is to use the PROC option JOINTMODEL. With this option, the log-likelihood function you specify applies not to a single observation but to the entire data set. See Modeling Joint Likelihood for details on how to properly use this option. The basic idea is that you store all necessary data set variables in arrays and use only the arrays to construct the log likelihood of the entire data set. This approach works here because for every observation , you can use index to access different values of arrays to construct the risk set S. To use the JOINTMODEL option, you need to do some additional data manipulation. You want to create a stop variable for each observation, which indicates the observation number that should be included in S for that observation. For example, if observations 4, 5, 6 all have the same survival time, the stop value for all of them is 6.

The following statements generate a new data set myelomaM that contains the stop variable:

data myelomaM;
   merge myelomaM freqs(drop=percent);
   by descending time;
   retain stop;
   if first.time then do;
      stop = _n_ + count - 1;
   end;
run;

The following SAS program fits the same Cox model by using the JOINTMODEL option:

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;

No output tables were produced because this PROC MCMC run produces identical posterior samples as does the previous example.

Because the JOINTMODEL option is specified here, you do not need to specify myelomaM as the input data set. An empty data set a is used to speed up the procedure run.

Multiple ARRAY statements allocate array symbols that are used to store the parameters (beta), the response and the covariates (data, timeA, vstatusA, and stopA), and the work space (bZ and S). The data, timeA, vstatusA, and stopA arrays are declared with the /NOSYMBOLS option. This option enables PROC MCMC to dynamically resize these arrays to match the dimensions of the input data set. See the section READ_ARRAY Function. The bZ and S arrays store the regression term and the risk set term for every observation.

The BEGINCNST and ENDCNST statements enclose programming statements that read the data set variables into these arrays. The rest of the programming statements construct the log likelihood for the entire data set.

The CALL MULT function calculates the regression term in the model and stores the result in the array bZ. In the first DO loop, you sum the risk set term S as if there are no ties in time. This underevaluates some of the S elements. For observations that have a tied time, you make the necessary correction to the corresponding S values. The correction takes place in the second DO loop. Any observation that has a tied time also has a stopA[i] that is different from i. You add the right terms to S and sum up the joint log likelihood jl. The MODEL statement specifies that the log likelihood takes on the value of jl.

To see that you get identical results from these two approaches, use PROC COMPARE to compare the posterior samples from two runs:

proc compare data=outi compare=outa;
   ods select comparesummary;
   var beta1-beta9;
run;

The output is not shown here.

Generally, the JOINTMODEL option can be slightly faster than using the default setup. The savings come from avoiding the overhead cost of accessing the data set repeatedly at every iteration. However, the speed gain is not guaranteed because it largely depends on the efficiency of your programs.

PROC PHREG fits the same model, and you get very similar results to PROC MCMC. The following statements run PROC PHREG and produce Output 52.8.2:

proc phreg data=Myeloma;
   ods select PostSummaries PostIntervals;
   ods output posteriorsample = phout;
   model Time*VStatus(0)=LogBUN HGB Platelet Age LogWBC
                         Frac LogPBM Protein Scalc;
   bayes seed=1 nmc=10000;
run;

Output 52.8.3 shows kernel density plots that compare the posterior marginal distributions of all the beta parameters from the PROC MCMC run and the PROC PHREG run. The following statements generate the comparison:

proc kde data=outi;
   ods exclude all;
   univar beta1 beta2 beta3 beta4 beta5 beta6 beta7 beta8 beta9
      / out=m1 (drop=count);
run;
ods exclude none;
 
%reshape(m1, mcmc, suffix1=, suffix2=md);
 
data phout;
   set phout(drop = LogPost Iteration);
   beta1 = LogBUN;  beta2 = HGB;      beta3 = Platelet;
   beta4 = Age;     beta5 = LogWBC;   beta6 = Frac;
   beta7 = LogPBM;  beta8 = Protein;  beta9 = SCalc;
   drop LogBUN HGB Platelet Age LogWBC Frac LogPBM Protein SCalc;
run;
 
proc kde data=phout;
   ods exclude all;
   univar beta1 beta2 beta3 beta4 beta5 beta6 beta7 beta8 beta9
      / out=m2 (drop=count);
run;
ods exclude none;
 
%reshape(m2, phreg, suffix1=p, suffix2=pd);
 
data all;
   merge mcmc phreg;
run;
 
proc template;
   define statgraph threebythree;
      %macro plot;
         begingraph;
            layout lattice / rows=3 columns=3;
               %do i = 1 %to 9;
                  layout overlay /yaxisopts=(label=" ");
                     seriesplot y=beta&i.md x=beta&i
                        / connectorder=xaxis
                          lineattrs=(pattern=mediumdash color=blue)
                        legendlabel = "MCMC" name="MCMC";
                     seriesplot y=beta&i.pd x=beta&i.p
                        / connectorder=xaxis lineattrs=(color=red)
                        legendlabel = "PHREG" name="PHREG";
                  endlayout;
               %end;
               Sidebar / align = bottom;
                  discretelegend "MCMC" "PHREG";
               endsidebar;
            endlayout;
         endgraph;
      %mend; %plot;
   end;
run;
 
proc sgrender data=all template=threebythree;
   title "Kernel Density Comparison";
run;

The macro %RESHAPE is defined in the example Logistic Regression Random-Effects Model. The posterior densities are almost identical to one another.

Output 52.8.3 Comparing Estimates from PROC MCMC and PROC PHREG

Time Dependent Model

To model as a function of the survival time, you can relate time to covariates by using this formula:

For illustrational purposes, only three explanatory variables, LOGBUN, HBG, and PLATELET, are used in this example.

Since depends on , every term in the summation of is a product of the current time and all observations that are in the risk set. You can use the JOINTMODEL option, as in the last example, or you can modify the input data set such that every row contains not only the current observation but also all observations that are in the corresponding risk set. When you construct the log likelihood for each observation, you have all the relevant data at your disposal.

The following statements illustrate how you can create a new data set with different risk sets at different rows:

title 'Cox Model with Time Dependent Covariates';
ods select none;
proc freq data=myeloma;
   tables time / out=freqs;
run;
ods select all;
 
proc sort data = freqs;
   by descending time;
run;
 
data myelomaM;
   set myeloma;
   ind = _N_;
run;
 
data myelomaM;
   merge myelomaM freqs(drop=percent); by descending time;
   retain stop;
   if first.time then do;
      stop = _n_ + count - 1;
   end;
run;
 
%macro array(list);
   %global mcmcarray;
   %let mcmcarray = ;
   %do i = 1 %to 32000;
      %let v = %scan(&list, &i, %str( ));
      %if %nrbquote(&v) ne %then %do;
         array _&v[&n];
         %let mcmcarray = &mcmcarray array _&v[&n] _&v.1 - _&v.&n%str(;);
            do i = 1 to stop;
               set myelomaM(keep=&v) point=i;
               _&v[i] = &v;
            end;
      %end;
      %else %let i = 32001;
   %end;
%mend;
 
data z;
   set myelomaM;
   %array(logbun hgb platelet);
   drop vstatus logbun hgb platelet count stop;
run;
 
data myelomaM;
   merge myelomaM z; by descending time;
run;

The data set myelomaM contains 65 observations and 209 variables. For each observation, you see added variables stop, _logbun1 through _logbun65, _hgb1 through _hgb65, and _platelet1 through _platelet65. The variable stop indicates the number of observations that are in the risk set of the current observation. The rest are transposed values of model covariates of the entire data set. The data set contains a number of missing values. This is due to the fact that only the relevant observations are kept, such as _logbun1 to _logbunstop. The rest of the cells are filled in with missing values. For example, the first observation has a unique survival time of 92 and stop is 1, making it a risk set of itself. You see nonmissing values only in _logbun1, _hgb1, and _platelet1.

The following statements fit the Cox model by using PROC MCMC:

proc mcmc data=myelomaM outpost=outi nmc=50000 ntu=3000 seed=17
          missing=ac;
   ods select PostSummaries PostIntervals;
   array beta[6];
   &mcmcarray
   parms (beta:) 0;
   prior beta: ~ normal(0, prec=1e-6);
 
   b = (beta1 + beta2 * time) * logbun +
       (beta3 + beta4 * time) * hgb +
       (beta5 + beta6 * time) * platelet;
   S = 0;
   do i = 1 to stop;
      S = S + exp( (beta1 + beta2 * time) * _logbun[i] +
                   (beta3 + beta4 * time) * _hgb[i] +
                   (beta5 + beta6 * time) * _platelet[i]);
   end;
   loglike = vstatus * (b - log(S));
 
   model general(loglike);
run;

Note that the option MISSING= is set to AC. This is due to missing cells in the input data set. You must use this option so that PROC MCMC retains observations that contain missing values.

The macro variable &mcmcarray is defined in the earlier part in this example. You can use a %put statement to print its value:

%put &mcmcarray;

This statement prints the following:

   array _logbun[65] _logbun1 - _logbun65; array _hgb[65] _hgb1 - _hgb65; array
   _platelet[65] _platelet1 - _platelet65;

The macro uses the ARRAY statement to allocate three arrays, each of which links their corresponding data set variables. This makes it easier to reference these data set variables in the program. The PARMS statement puts all the parameters in the same block. The PRIOR statement gives them normal priors with large variance. The symbol b is the regression term, and S is cumulatively added from 1 to stop for each observation in the DO loop. The symbol loglike completes the construction of log likelihood for each observation and the MODEL statement completes the model specification.

Posterior summary and interval statistics are shown in Output 52.8.4.

You can also use the option JOINTMODEL to get the same inference and avoid transposing the data for every observation:

proc mcmc data=myelomaM outpost=outa nmc=50000 ntu=3000 seed=17 jointmodel;
   ods select none;
   array beta[6];      array timeA[&n];      array vstatusA[&n];
   array logbunA[&n];  array hgbA[&n];       array plateletA[&n];
   array stopA[&n];    array bZ[&n];         array S[&n];
 
   begincnst;
      timeA[ind]=time;          vstatusA[ind]=vstatus;
      logbunA[ind]=logbun;      hgbA[ind]=hgb;
      plateletA[ind]=platelet;  stopA[ind]=stop;
   endcnst;
 
   parms (beta:) 0;
   prior beta: ~ normal(0, prec=1e-6);
 
   jl = 0;
   do i = 1 to &n;
      v1 = beta1 + beta2 * timeA[i];
      v2 = beta3 + beta4 * timeA[i];
      v3 = beta5 + beta6 * timeA[i];
      bZ[i] = v1 * logbunA[i] + v2 * hgbA[i] + v3 * plateletA[i];
 
      /* sum over risk set without considering ties in time. */
      S[i] = exp(bZ[i]);
      if (i > 1) then do;
         do j = 1 to (i-1);
            b1 = v1 * logbunA[j] + v2 * hgbA[j] + v3 * plateletA[j];
            S[i] = S[i] + exp(b1);
         end;
      end;
   end;
 
   /* make correction to the risk set due to ties in time. */
   do i = 1 to &n;
      if(stopA[i] > i) then do;
         v1 = beta1 + beta2 * timeA[i];
         v2 = beta3 + beta4 * timeA[i];
         v3 = beta5 + beta6 * timeA[i];
         do j = (i+1) to stopA[i];
            b1 = v1 * logbunA[j] + v2 * hgbA[j] + v3 * plateletA[j];
            S[i] = S[i] + exp(b1);
         end;
      end;
      jl = jl + vstatusA[i] * (bZ[i] - log(S[i]));
   end;
   model general(jl);
run;

The multiple ARRAY statements allocate array symbols that are used to store the parameters (beta), the response (timeA), the covariates (vstatusA, logbunA, hgbA, plateletA, and stopA), and work space (bZ and S). The bZ and S arrays store the regression term and the risk set term for every observation. Programming statements in the BEGINCNST and ENDCNST statements input the response and covariates from the data set to the arrays.

Using the same technique shown in the example Time Independent Model, the next DO loop calculates the regression term and corresponding S for every observation, pretending that there are no ties in time. This means that the risk set for observation involves only observation to . The correction terms are added to the corresponding S[i] in the second DO loop, conditional on whether the stop variable is greater than the observation count itself. The symbol jl cumulatively adds the log likelihood for the entire data set, and the MODEL statement specifies the joint log-likelihood function.

The following statements run PROC COMPARE and show that the output data set outa contains identical posterior samples as outi:

proc compare data=outi compare=outa;
   ods select comparesummary;
   var beta1-beta6;
run;

The results are not shown here.

The following statements use PROC PHREG to fit the same time dependent Cox model:

proc phreg data=Myeloma;
   ods select PostSummaries PostIntervals;
   ods output posteriorsample = phout;
   model Time*VStatus(0)=LogBUN z2 hgb z3 platelet z4;
   z2 = Time*logbun;
   z3 = Time*hgb;
   z4 = Time*platelet;
   bayes seed=1 nmc=10000;
run;

Coding is simpler than PROC MCMC. See Output 52.8.5 for posterior summary and interval statistics:

Output 52.8.6 shows a kernel density comparison plot that compares posterior marginal distributions of all the beta parameters from the PROC MCMC run and the PROC PHREG run. The following statements generate Output 52.8.6:

proc kde data=outi;
   ods exclude all;
   univar beta1 beta2 beta3 beta4 beta5 beta6 / out=m1 (drop=count);
run; 
ods exclude none;

%reshape(m1, mcmc, suffix1=, suffix2=md);
 
data phout;
   set phout(drop = LogPost Iteration);
   beta1 = LogBUN;  beta2 = z2;       beta3 = HGB;
   beta4 = z3;      beta5 = Platelet; beta6 = z4;
   drop LogBUN HGB Platelet z2-z4;
run;
 
proc kde data=phout;
   ods exclude all;
   univar beta1 beta2 beta3 beta4 beta5 beta6 / out=m2 (drop=count);
run;
ods exclude none;
 
%reshape(m2, phreg, suffix1=p, suffix2=pd);
 
data all;
   merge mcmc phreg;
run;
 
proc template;
   define statgraph twobythree;
      %macro plot;
         begingraph;
            layout lattice / rows=2 columns=3;
                %do i = 1 %to 6;
                  layout overlay /yaxisopts=(label=" ");
                     seriesplot y=beta&i.md x=beta&i
                        / connectorder=xaxis
                          lineattrs=(pattern=mediumdash color=blue)
                        legendlabel = "MCMC" name="MCMC";
                     seriesplot y=beta&i.pd x=beta&i.p
                        / connectorder=xaxis lineattrs=(color=red)
                        legendlabel = "PHREG" name="PHREG";
                  endlayout;
                %end;
                Sidebar / align = bottom;
                   discretelegend "MCMC" "PHREG";
                endsidebar;
            endlayout;
         endgraph;
      %mend; %plot;
   end;
run;
 
proc sgrender data=all template=twobythree;
   title "Kernel Density Comparison";
run;

The macro %RESHAPE is defined in the example Logistic Regression Random-Effects Model.

Output 52.8.6 Comparing Estimates from PROC MCMC and PROC PHREG