The PHREG Procedure

The proportional hazards spline model (Royston and Parmar 2002) with fixed covariate vector has a cumulative hazard function at time t,

where is the vector of regression coefficients and is a cubic spline function as described for the SPLINE option in the BAYES statement. The number of knots is equal to one plus the degrees of freedom (which you specify in the DF= option). The knots are determined by the sequence of the distinct event times of the data. PROC PHREG places the terminal knots at the minimum and maximum of the sequence and selects the m internal knots as the kth percentiles of the sequence, where . For example, if you specify DF=4, the three internal knots are the 25th, 50th, and 75th percentiles of the distinct event times.

Let be the observed data. The likelihood function of the proportional hazards spline model is given by

where and . Note that

Let be the observed data. Let be a partition of the time axis.

Hazards in Original Scale

The hazard function for subject i is

where

The baseline cumulative hazard function is

where

The log likelihood is given by

where .

Note that for , the full conditional for is log-concave only when , but the full conditionals for the ’s are always log-concave.

For a given , gives

Substituting these values into gives the profile log likelihood for

where . Since the constant c does not depend on , it can be discarded from in the optimization.

The MLE of is obtained by maximizing

with respect to , and the MLE of is given by

For , let

The partial derivatives of are

The asymptotic covariance matrix for is obtained as the inverse of the information matrix given by

See Example 6.5.1 in Lawless (2003) for details.

Hazards in Log Scale

By letting

you can build a prior correlation among the ’s by using a correlated prior , where .

The log likelihood is given by

Then the MLE of is given by

Note that the full conditionals for ’s and ’s are always log-concave.

The asymptotic covariance matrix for is obtained as the inverse of the information matrix formed by

For a Cox model, the model parameters are the regression coefficients. For a piecewise exponential model, the model parameters consist of the regression coefficients and the hazards or log-hazards. The priors for the hazards and the priors for the regression coefficients are assumed to be independent, while you can have a joint multivariate normal prior for the log-hazards and the regression coefficients. For a proportional hazards spline model, the model parameters consist of the regression coefficients and the cubic spline parameters. You can have a joint multivariate normal prior for the cubic spline parameters and the regression coefficients. Otherwise, the prior for the cubic spline parameters and the prior for the regression coefficients are assumed to be independent.

Cubic Spline Parameters

Uniform Prior

The joint prior density is given by

Normal Prior

Assume has a multivariate normal prior with mean vector and covariance matrix . The joint prior density is given by

Hazard Parameters

Let be the constant baseline hazards.

Improper Prior

The joint prior density is given by

This prior is improper (nonintegrable), but the posterior distribution is proper as long as there is at least one event time in each of the constant hazard intervals.

Uniform Prior

The joint prior density is given by

This prior is improper (nonintegrable), but the posteriors are proper as long as there is at least one event time in each of the constant hazard intervals.

Gamma Prior

The gamma distribution has a PDF

where a is the shape parameter and is the scale parameter. The mean is and the variance is .

Independent Gamma Prior

Suppose for , has an independent prior. The joint prior density is given by

AR1 Prior

are correlated as follows:

The joint prior density is given by

Log-Hazard Parameters

Write .

Uniform Prior

The joint prior density is given by

Note that the uniform prior for the log-hazards is the same as the improper prior for the hazards.

Normal Prior

Assume has a multivariate normal prior with mean vector and covariance matrix . The joint prior density is given by

Regression Coefficients

Let be the vector of regression coefficients.

Uniform Prior

The joint prior density is given by

This prior is improper, but the posterior distributions for are proper.

Normal Prior

Assume has a multivariate normal prior with mean vector and covariance matrix . The joint prior density is given by

Joint Multivariate Normal Prior for Log-Hazards and Regression Coefficients

Assume has a multivariate normal prior with mean vector and covariance matrix . The joint prior density is given by

Joint Multivariate Normal Prior for Cubic Spline Parameters and Regression Coefficients

Assume has a multivariate normal prior with mean vector and covariance matrix . The joint prior density is given by

Zellner’s g-Prior

Assume has a multivariate normal prior with mean vector and covariance matrix , where is the design matrix and g is either a constant or it follows a gamma prior with density where a and b are the SHAPE= and ISCALE= parameters. Let k be the rank of . The joint prior density with g being a constant c is given by

The joint prior density with g having a gamma prior is given by

Dispersion Parameter for Frailty Model

Improper Prior

The density is

Inverse Gamma Prior

The inverse gamma distribution has a density

where a and b are the SHAPE= and SCALE= parameters, respectively.

Gamma Prior

The gamma distribution has a density

where a and b are the SHAPE= and ISCALE= parameters, respectively.

Denote the observed data as D.

Cox Model

Proportional Hazards Spline Model

where is the likelihood function with cubic spline parameters and regression coefficients as parameters.

Frailty Model

Based on the framework of Sargent (1998),

where the joint density of the random effects is given by

Piecewise Exponential Model

Hazard Parameters

where is the likelihood function with hazards and regression coefficients as parameters.

Log-Hazard Parameters

where is the likelihood function with log-hazards and regression coefficients as parameters.

For the Gibbs sampler, PROC PHREG uses the ARMS (adaptive rejection Metropolis sampling) algorithm of Gilks, Best, and Tan (1995) to sample from the full conditionals. This is the default sampling scheme. Alternatively, you can requests the random walk Metropolis (RWM) algorithm to sample an entire parameter vector from the posterior distribution. For a general discussion of these algorithms, see section Markov Chain Monte Carlo Method in Chapter 8, Introduction to Bayesian Analysis Procedures.

You can output these posterior samples into a SAS data set by using the OUTPOST= option in the BAYES statement, or you can use the following SAS statement to output the posterior samples into the SAS data set Post:

 ods output PosteriorSample=Post;

The output data set also includes the variables LogLike and LogPost, which represent the log of the likelihood and the log of the posterior log density, respectively.

Let be the parameter vector. For the Cox model, the ’s are the regression coefficients ’s, and for the piecewise constant baseline hazard model, the ’s consist of the baseline hazards ’s (or log baseline hazards ’s) and the regression coefficients ’s. Let be the likelihood function, where D is the observed data. Note that for the Cox model, the likelihood contains the infinite-dimensional baseline hazard function, and the gamma process is perhaps the most commonly used prior process (Ibrahim, Chen, and Sinha 2001). However, Sinha, Ibrahim, and Chen (2003) justify using the partial likelihood as the likelihood function for the Bayesian analysis. Let be the prior distribution. The posterior is proportional to the joint distribution .

Gibbs Sampler

The full conditional distribution of is proportional to the joint distribution; that is,

For example, the one-dimensional conditional distribution of , given , is computed as

Suppose you have a set of arbitrary starting values . Using the ARMS algorithm, an iteration of the Gibbs sampler consists of the following:

• draw from

• draw from

• 

• draw from

After one iteration, you have . After n iterations, you have . Cumulatively, a chain of n samples is obtained.

When the BAYES statement is specified, PROC PHREG generates one Markov chain that contains the approximate posterior samples of the model parameters. Additional chains are produced when the Gelman-Rubin diagnostics are requested. Starting values (initial values) can be specified in the INITIAL= data set in the BAYES statement. If the INITIAL= option is not specified, PROC PHREG picks its own initial values for the chains based on the maximum likelihood estimates of and the prior information of .

Denote as the integral value of x.

Constant Baseline Hazard Parameters ’s

For the first chain that the summary statistics and diagnostics are based on, the initial values are

For subsequent chains, the starting values are picked in two different ways according to the total number of chains specified. If the total number of chains specified is less than or equal to 10, initial values of the rth chain () are given by

with the plus sign for odd r and minus sign for even r. If the total number of chains is greater than 10, initial values are picked at random over a wide range of values. Let be a uniform random number between 0 and 1; the initial value for is given by

Regression Coefficients and Log-Hazard Parameters ’s

The ’s are the regression coefficients ’s, and in the piecewise exponential model, include the log-hazard parameters ’s. For the first chain that the summary statistics and regression diagnostics are based on, the initial values are

If the number of chains requested is less than or equal to 10, initial values for the rth chain () are given by

with the plus sign for odd r and minus sign for even r. When there are more than 10 chains, the initial value for the is picked at random over the range ; that is,

where is a uniform random number between 0 and 1.

Denote the observed data by D. Let be the vector of parameters of length k. Let be the likelihood. The deviance information criterion (DIC) proposed in Spiegelhalter et al. (2002) is a Bayesian model assessment tool. Let Dev. Let and be the corresponding posterior means of and , respectively. The deviance information criterion is computed as

Also computed is

where pD is interpreted as the effective number of parameters.

Note that defined here does not have the standardizing term as in the section Deviance Information Criterion (DIC) in Chapter 8, Introduction to Bayesian Analysis Procedures. Nevertheless, the DIC calculated here is still useful for variable selection.

Let be the parameter vector. For the Cox model, the ’s are the regression coefficients ’s; for the proportional hazards spline model, the ’s consist of the cubic spline parameters ’s and the regression coefficients ’s; for the piecewise constant baseline hazard model, the ’s consist of the baseline hazards ’s (or log baseline hazards ’s) and the regression coefficients ’s.

Let be the chain that represents the posterior distribution for .

Consider a quantity of interest that can be expressed as a function of the parameter vector . You can construct the posterior distribution of by evaluating the function for each in . The posterior chain for is Summary statistics such as mean, standard deviation, percentiles, and credible intervals are used to describe the posterior distribution of .

Survival Distribution

Let be a covariate vector of interest.

Cox Model

Let be the observed data. Define

Consider the rth draw of . The baseline cumulative hazard function at time t is given by

For the given covariate vector , the cumulative hazard function at time t is

and the survival function at time t is

Proportional Hazards Spline Model

Consider the rth draw in , where consists of and . The baseline cumulative hazard function at time t is

where is a cubic spline function as described for the SPLINE option in the BAYES statement. For the given covariate vector , the cumulative hazard function at time t is

and the survival function at time t is

Piecewise Exponential Model

Let be a partition of the time axis. Consider the rth draw in , where consists of and . The baseline cumulative hazard function at time t is

where

For the given covariate vector , the cumulative hazard function at time t is

and the survival function at time t is