Introduction to brea

Hazard Rates and Cox Proportional Hazards Models

With continuous time, the most common approach for analyzing event times is to model the hazard rate of event occurrence. Let \(T\) denote a continuous event time variable, taking some value in the interval \((0,\infty)\). The continuous time hazard rate \(\lambda(t)\) is defined to be the “instantaneous risk” of event occurrence at time \(t\) given the event of interest has not yet occurred prior to \(t\): \[ \lambda(t)=\lim_{dt\rightarrow 0} \frac{P(t\leq T \leq t+dt)}{dt\cdot P(T\geq t)} \] The continuous time Cox proportional hazards model then models these hazard rates with: \[ \lambda(t)=\lambda_0(t)\text{exp}(\beta_1X_1+\cdots\beta_MX_M)=\text{exp}(f_0(t)+\beta_1X_1+\cdots\beta_MX_M) \] where \(\lambda_0(t)\) is an arbitrary function of time \(t\) called a baseline hazard (and \(f_0(t)=\text{log}(\lambda_0(t))\) is this same baseline hazard on a log scale) and \(X_1,\ldots,X_M\) are covariate values (Cox 1972). Taking logarithms of both sides of this equation yields an equivalent form of our continuous time Cox model: \[ \text{log}(\lambda(t)) = \eta(t) = f_0(t)+\beta_1X_1+\cdots\beta_MX_M \] where \(\eta(t)\) is the linear predictor incorporating both the time effects and other covariate effects.

As with continuous time, with discrete time we will also model the hazard rate of event occurrence. Let \(T\) now denote the discrete time of event occurrence, and we assume the possible discrete time points are labeled with the positive integers. So for example, if the event of interest occurs at the third discrete time point, then \(T=3\). If there is only a single event type (i.e., no competing risks), we define the discrete time hazard rate \(\lambda(t)\) to be the probability of event occurrence at time \(t\) given the event has not already occurred at an earlier time point: \[\lambda(t)=P(T=t|T\geq t)\] Unlike with the continuous time hazard rate (which may take any value between 0 and \(\infty\)), the discrete time hazard rate is a probability and is thus bounded between 0 and 1. Hence, we relate \(\lambda(t)\) to a linear predictor \(\eta(t)\) encapsulating covariate effects (including the effect of time \(t\)) using the logit link function (just as in logistic regression): \[ \text{log}\left(\frac{\lambda(t)}{1-\lambda(t)}\right) = \eta(t) = f_0(t)+\beta_1X_1+\cdots\beta_MX_M \] Here, for the discrete time model, proportionality is now technically in the odds of the hazard. However, for simplicity we will still refer to this as the discrete time Cox proportional hazards model.

Expanding Data to Person-Time Format

The discrete Cox model is recognizable as being identical to logistic regression. However, the difference here is that we must include one logistic regression observation for each discrete time point \(t\) each sample member was at risk for event occurrence in order to model the aggregate risk experienced by each sample member. To illustrate, suppose we have the following data set:

where ID is a subject id number, \(T\) is study time (the time of either event occurrence or right censoring), \(\delta\) is an event indicator (1 for event occurrence, 0 for right censoring), and age and gender are covariates. For instance, the first study subject did not experience the event of interest at times \(t=1\) or \(t=2\) but did experience the event of interest at time \(t=3\).

To fit the discrete Cox model to this data, we must first expand the data set to person-time format, with one row of data for each discrete time point each sample member was at risk for event occurrence. The above data set would become:

where \(t\) is the discrete time for that specific person-time observation, and \(Y\) is an outcome indicator (\(Y=0\) for no event occurring at that person-time point, and \(Y=1\) for event occurrence). Note that when \(t=T\) (i.e., it is the last discrete time point of observation for a sample member), \(Y=\delta\). For example, study subject 02 was observed for two discrete time points, and because her \(\delta\) value is 0, her \(Y\) value at her second person-time point of observation is also 0. For all three study subjects depicted here, we assumed their ages remained constant across all discrete time points of observation. However, if we knew otherwise, we could easily calculate potentially distinct age values for each person-time observation.

Formatting Data for the brea Package

Once a data set has been expanded into person-time format, it may be fit via conventional logistic regression software (e.g., the glm() function with the family=binomial option; see the example later in this tutorial). The brea_mcmc() function from the brea package is based on logistic regression, but has one additional requirement: All predictor variables (including time \(t\)) must be specified in a categorical manner, either by giving positive integer values \(k\) (where \(k=1,\ldots,K\) for some \(K\)) or by using an R factor data type. The logistic regression model implementing the discrete Cox model then becomes:

\[ \text{log}\left(\frac{\lambda(t)}{1-\lambda(t)}\right) = \eta(t) = \beta_0 + \beta_1[X_1(t)] + \beta_2[X_2(t)] + \cdots + \beta_M[X_m(t)] \] Here each predictor value \(X_m(t)\) takes a positive integer value \(k\) between \(1\) and \(K_m\), and so the effects of the \(m^\text{th}\) predictor are modeled on the logit scale by a collection of \(K_m\) parameters \(\beta_m[1]\), \(\beta_m[2]\),…,\(\beta_m[K_m]\).

Note that while the specification of covariate values must be categorical, quantitative covariates are still fully supported. To include a quantitative covariate, the user must first discretize the covariate (e.g., using the cut() function). This effectively means the potentially arbitrary smooth relationship between the covariate and the hazard will be approximated by a step function. To increase the amount of smoothness, a finer discretization (with a larger number \(K_m\) of categories) may be used. Gaussian Markov random field (GMRF) prior distributions are then used for the parameter values \(\beta_m[1]\), \(\beta_m[2]\),…,\(\beta_m[K_m]\) to ensure regularity/smoothness of the functional relationship even if the discretization is fine. See (King and Weiss 2021) and the practical example later in this tutorial for more details.

Prior and Posterior Distributions

The probability distribution we use to capture knowledge we have about model parameters prior to seeing the sample data is called the prior distribution. For example, if had reason to believe that a model parameter \(\mu\) was equal to 5, then we could use a prior distribution for \(\mu\) with a mean of 5. Furthermore, if we were roughly 95% sure (prior to seeing the sample data) that \(\mu\) is between 1 and 9, then we could use a \(N(5,2^2)\) (i.e., a normal distribution with a mean of 5 and a standard deviation of 2) prior distribution for \(\mu\), since for this distribution roughly 95% of its probability is between 1 and 9:

pnorm(9,5,2)-pnorm(1,5,2)

## [1] 0.9544997

qnorm(c(0.025,0.975),5,2)

## [1] 1.080072 8.919928

With our modeling approach, we will generally use noninformative priors that specify little, if any, prior information about the model parameters (other than, for example, the fact that adjacent step heights in a step function approximation of a smooth curve are close together). We will illustrate how to do this in the practical example, and users of the brea package will generally not have to be concerned with developing prior distribution specifications beyond “default” noninformative priors.

The posterior distribution combines the information about the model parameters from the prior distribution with the information from the sample data to produce a single probability distribution that encapsulates our most up-to-date state of knowledge about the parameters. The posterior distribution is the primary “output” or “result” or “inference” in a Bayesian analysis.

Inference via MCMC Sampling from the Posterior Distribution

Mathematically, the posterior distribution is completely determined by the combination of the prior distribution, statistical model, and sample data. In all but the simplest cases, this probability distribution doesn’t have some simple form that we can determine using algebra or calculus. Instead, the only practical way to learn about this distribution is to take a random sample from the posterior distribution. Note that here we are talking about using computer simulation algorithms to draw a random sample from the posterior probability distribution of the model parameters, and this sample is completely different from the sample data that we are analyzing.

The most common way of selecting a random sample from the posterior distribution is via a numerical algorithm called Markov chain Monte Carlo, or MCMC for short. This is the algorithm used by the brea_mcmc() function we will use later. One important difference between a random sample drawn using an MCMC algorithm and a simple random sample (SRS) that the user may be familiar with is that unlike with an SRS sample, the values from an MCMC sample are not independent of each other. Specifically, MCMC generates values for the sample as a sequence, and successive elements of the sequence tend to be more similar to each other than independently drawn values (i.e., the sequence of sampled values exhibits positive autocorrelation).

As a result, for a given sample size \(S\), there is usually much less information in a sample drawn using MCMC than an independent SRS sample. So while we may draw a sample of size \(S\) using MCMC, the effective sample size (the size of an independent SRS sample that would contain the same amount of information as the MCMC sample) is often 10 or even 100 times smaller than \(S\). In the practical example later, we will show how to calculate the effective sample size using the effectiveSize() function from the coda package; for publication-quality inferences, we recommend choosing the MCMC sample size \(S\) so that the effective sample size is at least 1,000 (which may necessitate choosing \(S=10,000\) or even \(S=100,000\)).

Once we have obtained a sample from the posterior distribution, it is very easy to obtain the types of inferential results the user is familiar with from fitting frequentist models. Specifically, we may obtain a point estimate of a specific model parameter by calculating the mean or median of the MCMC sample for that parameter. Similarly, we may obtain a confidence interval (CI, which Bayesians often call credible intervals) for a model parameter by calculating quantiles of the sampled values of that parameter. For example, we could find the endpoints of a 95% CI by calculating the \(2.5^\text{th}\) percentile and \(97.5^\text{th}\) percentile of the posterior MCMC sample for the desired parameter. We illustrate this in the practical example below.

Polynomial Time Effect

With linear modeling of the effect of time \(t\), our model becomes: \[ \text{log}\left(\frac{\lambda(t)}{1-\lambda(t)}\right) = \beta_0 + \beta_1X_1(t) + \beta_2X_2 \] where \(X_1(t)=t\) is discrete time and \(X_2\) is treatment group (0 for 6-MP, 1 for placebo). We may fit this logistic regression model to our person-time data set using the glm() function with the family=binomial option:

linear_fit <- glm(y ~ t + treat,family=binomial,data=expanded_data)
summary(linear_fit)

## 
## Call:
## glm(formula = y ~ t + treat, family = binomial, data = expanded_data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.01627    0.50971  -7.880 3.29e-15 ***
## t            0.02768    0.02745   1.008    0.313    
## treat        1.77337    0.44334   4.000 6.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 231.84  on 540  degrees of freedom
## Residual deviance: 213.32  on 538  degrees of freedom
## AIC: 219.32
## 
## Number of Fisher Scoring iterations: 6

We may then obtain a point estimate of the hazard rate ratio (the exponentiated coefficient of the treatment variable \(X_2\)) using:

beta2_hat <- coef(linear_fit)["treat"]
exp(beta2_hat)

##    treat 
## 5.890668

Note that technically this is an odds ratio estimate, not a hazard ratio estimate, since the discrete Cox model sets the log of the odds of the hazard equal to the linear predictor \(\eta(t)\). Odds ratios will generally tend to be more extreme (further from 1) than corresponding risk/hazard ratios, though when the probabilities are small (as they are for most discrete survival models), the differences are modest. We will informally still interpret odds ratio estimates from discrete Cox models as though they were hazard rate ratios.

We may also try quadratic and cubic formulations for the time effect:

quadratic_fit <- glm(y ~ poly(t,2) + treat,family=binomial,data=expanded_data)
exp(coef(quadratic_fit)["treat"])

##    treat 
## 5.550794

cubic_fit <- glm(y ~ poly(t,3) + treat,family=binomial,data=expanded_data)
exp(coef(cubic_fit)["treat"])

##    treat 
## 5.530024

Step Function Time Effect

Next, we consider grouping discrete time t into non-overlapping intervals and treating the variable as categorical. This is equivalent to approximating the baseline hazard \(f_0(t)\) with a step function. When doing frequentist inference with this approach, we must have at least one observed event in each interval of \(t\) values, or else the corresponding parameter estimate for that interval will be \(-\infty\) (since the hazard appears to be zero in that interval). Our dataset only contains 30 observed events, and the last time of event occurrence is \(t=23\). With such limited data, we follow the suggestion of (Cox 1972) of dividing the variable range into intervals of length 10 (with the exception of the last interval, which will have length 15 in order to include the last time of observation):

expanded_data$t_cat <- cut(expanded_data$t,c(0,10,20,35))

By doing a cross tabulation of the discrete time variable t in the expanded person-time data set with the new grouped version t_cat, we can verify that all discrete time observations have been grouped into the correct intervals:

with(expanded_data,table(t,t_cat))

##     t_cat
## t    (0,10] (10,20] (20,35]
##   1      42       0       0
##   2      40       0       0
##   3      38       0       0
##   4      37       0       0
##   5      35       0       0
##   6      33       0       0
##   7      29       0       0
##   8      28       0       0
##   9      24       0       0
##   10     23       0       0
##   11      0      21       0
##   12      0      18       0
##   13      0      16       0
##   14      0      15       0
##   15      0      15       0
##   16      0      14       0
##   17      0      13       0
##   18      0      11       0
##   19      0      11       0
##   20      0      10       0
##   21      0       0       9
##   22      0       0       9
##   23      0       0       7
##   24      0       0       5
##   25      0       0       5
##   26      0       0       4
##   27      0       0       4
##   28      0       0       4
##   29      0       0       4
##   30      0       0       4
##   31      0       0       4
##   32      0       0       4
##   33      0       0       2
##   34      0       0       2
##   35      0       0       1

Indeed, we observe for example that all observations with t values between 1 and 10 have been correctly assigned to the interval (0,10]. Our discrete time Cox model with a step function baseline hazard is now as follows: \[ \text{log}\left(\frac{\lambda(t)}{1-\lambda(t)}\right) = \beta_0 + \beta_1X_1(t) + \beta_2X_2(t) + \beta_3X_3 \] where \(X_1(t)\) is a dummy/indicator variable for the second interval \((10,20]\) (i.e., \(X_1(t)=1\) if \(t\in (10,20]\) and \(X_1(t)=0\) otherwise), \(X_2(t)\) is a dummy/indicator variable for the third interval \((20,30]\) (so the first interval \((0,10]\) is the reference group), and \(X_3\) is the treatment group variable. We may now fit this model:

step_fit <- glm(y ~ t_cat + treat,family=binomial,data=expanded_data)
summary(step_fit)

## 
## Call:
## glm(formula = y ~ t_cat + treat, family = binomial, data = expanded_data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.9850     0.4311  -9.244  < 2e-16 ***
## t_cat(10,20]   0.3326     0.4526   0.735    0.462    
## t_cat(20,35]   0.9364     0.6306   1.485    0.138    
## treat          1.8372     0.4459   4.120 3.79e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 231.84  on 540  degrees of freedom
## Residual deviance: 212.16  on 537  degrees of freedom
## AIC: 220.16
## 
## Number of Fisher Scoring iterations: 6

Our estimated hazard rate ratio \(\text{exp}(\hat \beta_3)\) for the treatment effect is:

exp(coef(step_fit)["treat"])

##    treat 
## 6.279227

We may also obtain a 95% CI for the hazard rate ratio as follows:

confint(step_fit)

## Waiting for profiling to be done...

##                   2.5 %    97.5 %
## (Intercept)  -4.9127924 -3.212701
## t_cat(10,20] -0.6038059  1.193752
## t_cat(20,35] -0.4157477  2.121095
## treat         1.0011180  2.768136

exp(confint(step_fit)["treat",])

## Waiting for profiling to be done...

##     2.5 %    97.5 % 
##  2.721323 15.928911

Hence, using this model we estimate that the hazard of remission ending is 6.28 times higher in the placebo group compared to the 6-MP group. In addition, we can be 95% sure the hazard is between 2.72 and 15.93 times higher among placebo patients when compared to patients receiving the active treatment.

Note that the effect estimate here is substantially larger in magnitude than for any of the other formulations we have considered thus far. This is due to our step function approximation of the baseline hazard (the time effect) being too crude (only having three steps). With our Bayesian approach in the next section, we will see how to fix this.

Setting Up an MCMC Run

The brea_mcmc() function has the following function signature:

brea_mcmc(x, y, priors = NULL, S = 1000, B = 100, n = NULL, K = NULL, store_re = FALSE)

The only required arguments the user must specify are x and y, which specify, respectively, the predictor values and event occurrence status for each person-time observation of the expanded data set. More specifically, y is a vector of length N or a matrix of height N with a single column (or with multiple columns if there are competing risks), where N is the number of person-time observations. Each entry of y is either 0 (if the event of interest did not occur at that person-time observation) or 1 (if the event of interest did occur). The y variable we created earlier meets these requirements without further modification.

The x variable used to specify covariates at the person-time level may either be a matrix or dataframe with N rows. All predictors must use a categorical formulation (though often the categories represent steps in a step function approximation of an arbitrary smooth function of a numerical predictor). The categories of each categorical predictor are represented either by positive whole number values or by R factor variables (factors are only possible if x is a data frame; if x is a matrix, then x must be numeric).

We will now construct a numeric matrix x that meets the above requirements. First, we will use the same grouping of discrete time t into three intervals as used in the previous section. We will also change the coding of the treatment indicator to take values 1 and 2 instead of 0 and 1:

x_brea <- matrix(0,N,2)
x_brea[,1] <- cut(expanded_data$t,c(0,10,20,35),labels=FALSE) # grouped time t
x_brea[,2] <- expanded_data$treat + 1 # treatment group

Typically, users will also specify the priors argument, which determines both the covariate type (categorical fixed effects, quantitative covariates with step function approximations of the effects, or random effects) and any parameters for its corresponding prior distribution. However, for this first brea model, we will be treating both predictors as categorical fixed effects and be using corresponding noninformative priors. Since this is the default specification (when NULL is supplied for the priors argument), we do not have to provide this argument.

We may then load the brea package and obtain a sample from the posterior distribution of the model parameters (“fitting” the model in the Bayesian sense) as follows:

library(brea)
set.seed(1234)
fit <- brea_mcmc(x_brea,y)

Note that since the MCMC algorithm uses random number generation, we set the seed of the random number generator to ensure our results are reproducible.

Extracting and Examining Results

The structure of the R object returned by brea_mcm() is below:

str(fit)

## List of 4
##  $ b_0_s: num [1, 1:1000] -2.77 -3.1 -2.98 -2.95 -2.48 ...
##  $ b_m_s:List of 2
##   ..$ : num [1, 1:3, 1:1000] -0.2452 -0.059 0.3043 0.0227 0.2089 ...
##   ..$ : num [1, 1:2, 1:1000] -0.795 0.795 -0.735 0.735 -0.879 ...
##  $ s_m_s:List of 2
##   ..$ : NULL
##   ..$ : NULL
##  $ b_m_a:List of 2
##   ..$ : int [1:3] 450 446 483
##   ..$ : int [1:2] 457 470

This object is a list with two components of interest to us (as well as a few that do not concern us now): b_0_s contains the posterior sampled values of the intercept parameter \(\beta_0\), and b_m_s contains posterior sampled values of the parameters \(\beta_m[k]\), where \(\beta_m[k]\) represents the effect of the \(k^\text{th}\) covariate value (category) of the \(m^\text{th}\) predictor on the hazard of event occurrence. Specifically, the \(s^\text{th}\) sampled value (post burn-in) of \(\beta_m[k]\) is stored in b_m_s[[m]][1,k,s] (the 1 index in [1,k,s] denotes that we are examining the first competing risk, which must be specified even in cases like this where there is only one competing risk). In particular, for our model in question, m=2 represents the treatment variable (since treatment is the second column of the x_brea matrix we supplied), k=1 represents the 6-MP treatment, and k=2 represents the placebo control group. Hence, MCMC samples of the treatment effect (on the linear predictor scale) may be calculated as follows:

b_treatment <- fit$b_m_s[[2]][1,1,]
b_control <- fit$b_m_s[[2]][1,2,]
d <- b_control - b_treatment # sampled values of treatment effect on logit scale

We may then obtain point estimates and a standard error estimate of the treatment effect on the linear predictor (logit) scale, and compare it to the results from fitting a frequentist version of this model using glm() earlier:

mean(d) # posterior mean point estimate

## [1] 1.873343

median(d) # posterior median point estimate

## [1] 1.841043

sd(d) # posterior standard deviation (standard error)

## [1] 0.439673

summary(step_fit) # identical frequentist model fit with glm() earlier

## 
## Call:
## glm(formula = y ~ t_cat + treat, family = binomial, data = expanded_data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.9850     0.4311  -9.244  < 2e-16 ***
## t_cat(10,20]   0.3326     0.4526   0.735    0.462    
## t_cat(20,35]   0.9364     0.6306   1.485    0.138    
## treat          1.8372     0.4459   4.120 3.79e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 231.84  on 540  degrees of freedom
## Residual deviance: 212.16  on 537  degrees of freedom
## AIC: 220.16
## 
## Number of Fisher Scoring iterations: 6

As we can see, “fitting” the Bayesian version of this discrete Cox model gives almost exactly the same results as fitting the frequentist version using glm().

Assessing MCMC Performance and Adjusting MCMC Run Length S

As described earlier, with our Bayesian approach, we obtain inferences by taking a sample from the posterior distribution of the relevant model parameters using an MCMC algorithm. How accurately this posterior sample describes the posterior distribution depends on two things: how much serial correlation (autocorrelation) is present in the MCMC sample, and the size of the MCMC sample \(S\) (which is the length of the MCMC chain). We may visualize the sequence of sampled values using a trace plot, which plots the MCMC iteration number along the horizontal axis and the sampled parameter value along the vertical axis:

par(cex=0.66,mgp=c(1.75,0.5,0),mai=c(0.4,0.4,0.1,0.1))
plot(d,type="l",xlab="MCMC Interation Number s",
     ylab="Sampled Value of Treatment Effect Parameter")

From the plot, the serial correlation does not seem too severe. As mentioned earlier though, it’s useful to calculate an effective sample size estimate, which is the size of an independent sample that contains the same information as contained in our MCMC sample.

library(coda)
effectiveSize(d)

##     var1 
## 207.7105

The brea_mcmc() function uses a default posterior sample size of \(S=1,000\), and this sample produced an effective sample size of 207.7105468, meaning our efficiency was around 20%. As mentioned earlier, we typically should seek to have effective sample sizes over 1,000, and since the MCMC run is already very fast, we can increase the chain length to \(S=10,000\) (and also increase the initial portion of the chain discarded as burn-in to \(B=1,000\)):

set.seed(1234)
fit_10k <- brea_mcmc(x_brea,y,S=10000,B=1000)
d <- fit_10k$b_m_s[[2]][1,2,] - fit_10k$b_m_s[[2]][1,1,]
effectiveSize(d)

##     var1 
## 1695.931

mean(d)

## [1] 1.847461

median(d)

## [1] 1.834604

sd(d)

## [1] 0.4576245

exp(median(d))

## [1] 6.262654

Our effective sample size estimate is now comfortably over 1,000, though our estimates are not substantially changed (and are still quite close to the corresponding frequentist estimates, as we would expect).

Better Step Function Approximation using the GMRF Prior Type

We hypothesized above that our step function approximation of the time effect which only used three steps was too crude to adequately model the time effect. However, with the frequentist approach, we were restrained from including a more fine-grained step function by the limited amount of data available. With our Bayesian approach however, we can circumvent this restriction by using a prior distribution that mathematically “ties together” adjacent steps, which allows us to estimate step heights even if there is little to no data directly influencing each step height. One such prior distribution is called a Gaussian Markov random field (GMRF), and the brea package includes this as a possible prior type.

First, we will re-group the discrete time variable \(t\) into intervals of length 3 instead of intervals of length 10:

x_brea[,1] <- cut(expanded_data$t,seq(0,36,3),labels=FALSE)

Next, we must explicitly set the priors argument of the brea_mcmc() function; priors is a list with one element for each covariate (column of x). Each element, in turn, is a list where the first element specifies the covariate type (gmrf for quantitative covariates via GMRF, cat for categorical fixed effects, or re for random effects). The gmrf option requires two additional numerical parameters; specification of these is technical, and we recommend values of 3 and 0.01 as a “default” noninformative prior specification that should work well in most cases. Similarly, the cat option requires one additional numerical parameter, and a value of 4 results in a fairly noninformative prior distribution. We may now rerun the MCMC as follows:

set.seed(1234)
priors <- list(list("gmrf",3,0.01),list("cat",4))
fit_gmrf <- brea_mcmc(x_brea,y,priors,S=10000,B=1000)
d_gmrf <- fit_gmrf$b_m_s[[2]][1,2,] - fit_gmrf$b_m_s[[2]][1,1,]
effectiveSize(d_gmrf)

##     var1 
## 2058.032

The effective sample for this MCMC run is once again well over 1,000, so we may now proceed to giving practical inferences.

Inferences

The posterior mean, median, and standard deviation (standard error) of the treatment effect parameter on the linear predictor scale are as follows:

mean(d_gmrf)

## [1] 1.674498

median(d_gmrf)

## [1] 1.669348

sd(d_gmrf)

## [1] 0.4155349

Exponentiating, we find the following hazard ratio estimate:

median(exp(d_gmrf))

## [1] 5.308705

We estimate the hazard of cancer remission cessation is 5.31 times lower for patients receiving 6-MP compared to patients receiving placebo. We may also calculate a 95% CI for the hazard rate ratio by calculating the \(2.5^\text{th}\) and \(97.5^\text{th}\) percentile of the posterior sample of exponentiated treatment effect parameters:

quantile(exp(d_gmrf),c(0.025,0.975))

##      2.5%     97.5% 
##  2.412046 12.557243

Hence, based on these results, we are 95% confident that 6-MP treatment reduces the hazard of remission ending by a factor of between 2.41 and 12.56.

Comparison with Other Approaches

We have obtained inferential results from a total of nine different model fits to the leukemia remission data from (Cox 1972). First, we fit a classic frequentist continuous time Cox proportional hazards model with three different options for the handling of tied observations (the Efron, Breslow, and exact methods). Second, we fit a frequentist discrete time Cox model with four different ways of modeling the effect of discrete time \(t\): linear, quadratic, and cubic polynomials, and a step function with three steps. Finally, we obtained inferences from a Bayesian discrete time Cox model using the brea package with two different time effect formulations: a step function with 3 steps and a step function with 12 steps (the latter case using a GMRF prior distribution for the step heights).

For all nine of these model fits, we obtained point estimates and standard errors of the treatment effect parameter on the linear predictor scale. We also obtained corresponding hazard rate ratio estimates comparing the rate of remission cessation in the placebo control group to the rate in the 6-MP treatment group. These results are summarized in the table below:

While the results of all nine model fits were qualitatively similar, there were clear patterns regarding which results were most similar to others. For example, the frequentist and Bayesian discrete time Cox models the 3-step formulation for the time effect yielded almost identical results, which makes sense considering these are essentially the same model but fit in drastically different ways.

In addition, the continuous time Cox model fit with the exact method of handling ties performed quite similarly to the brea model with a flexible (12-step) formulation for the baseline hazard function. Again, this makes sense, as the models are quite similar, with the main difference being that we expect the log odds ratio estimate for the brea model (1.67) to be slightly further from the null value of 0 than a corresponding log hazard ratio from the continuous CPH model (1.63), which is what was observed.