# Sampling distribution is skewed in a fully Bayesian inference of MCMC in Cox PH models

I used the MCMC method to estimate linear models with a fully Bayesian inference previously, and had no problem from estimated coefficients. Recently I use the same way in a semiparametric Cox PH model, and get many skewed posterior distributions which make the posterior means of parameters pretty close to 2.5th-quantile or 97.5th-quantile. I try to increase the number of iteration (30,000 -> 50,000) and step (20 -> 40) with a burnin = 10,000 to make the posterior probability more symmetric. Sometimes it becomes symmetric, but sometimes it is still skewed eventually. I also use conditional prior proposals and iteratively weighted least squares based on posterior modes to increase the acceptance rate, but it is still useless. The data contains over 8,000 cancer patients with very few missing data. The independent variables are either binary or scale. So far I haven't figure out a best way to make the posterior distribution robustly symmetric. What I can do is repeatedly fitting the model to see whether the posterior probability is not skewed by chance. The reason I feel sick is that the model-fitting is so time-consuming. Each model spent 45mins~75mins with different combinations of variables, and I have 8 models in three specific causes. Each failed model means 1-hour waste of my life. Hopefully some people can share your experiences of dealing with this problem. The corresponding software is BayesX. I appreciate any advice!

Edit

Xi'an,

Because I always get symmetric posterior distribution of estimated parameters in linear mixed model by MCMC, I assume the Cox PH model should have the same scenario. Here are two results from a Cox PH model with my data:

Model 1: h(t,z)=ho(t)exp{intercept + sex + white + agedx + history}

  Variable      mean    Std. Dev.  2.5% quant.  median    97.5% quant.
const    -4.2486    0.29328    -4.70006   0.249887    0.068665
sex    0.0123228    0.0514687  -4.52946   0.0175766   0.0982832
white  -0.239128    0.0644409  -4.51786  -0.195356    0.0779398
agedx 0.00786072    0.00370628 -4.53803   0.00139     0.0734418
history 0.627554    0.1464     -4.55214   0.0027441   0.800568


Model 2: h(t,z)=ho(t)exp{intercept + sex + white + agedx + test}

 Variable      mean    Std. Dev.  2.5% quant.  median    97.5% quant.
const   -4.12811    0.305912    -4.74401   -4.08631   0.134946
sex    0.0138823    0.0532634   -0.219711  0.0135219  0.148526
white  -0.227753    0.0639598   -0.372975 -0.223358   0.111341
agedx 0.00737788    0.0037846   -0.214211  0.00705715 0.101075
test   -0.151933    0.0518964   -0.271067 -0.152058   0.101187


The have the same variables of sex, white and agedx. Only the last variable is not the same. The posterior means of sex, white and agedx are similar, but the 2.5% quantile and 97.5% quantile are different a lot. If the posterior distribution cannot be symmetric, which result is reliable?

• Why are you trying to make the posterior distribution symmetric in the first place? Nov 4, 2011 at 19:25
• Also, are you looking at symmetry on the log scale? Nov 4, 2011 at 19:52
• Looking at your results sets off another alarm - typically, it's my understanding that a Cox Proportional Hazards model doesn't estimate a baseline hazard. As such there is no intercept. The presence of one in your output makes me nervous. Nov 6, 2011 at 0:32
• To EpiGrad: BayesX fits a smoothing function for the baseline function. If the estimates are problematic like model 1, the baseline is wiggly, and it is explicit that some points on the curve tends to either upper bound or lower bound of its 95% credible intervals. Hence, I can conclude that the unknown influence indeed affects both linear predictor and baseline function. In these days, I found that the posterior distribution of estimates can have somewhat possibility to be symmetric if I set some specific priors for them. Even this is not a perfect solution, I will use this way temporarily. Nov 9, 2011 at 18:15

## 1 Answer

Just from reading your question, it seems to me you are mixing computational performances and statistical inference. The posterior distribution on the parameters does not have to be symmetric, so this is not an indicator of poor MCMC convergence. My advice is to try your code on simulated data (meaning simulated from the very semiparametric Cox model you are estimating) where you know the values of the parameters, in order to check for convergence: posteriors should cover the true values to some extent. And take comfort in the fact that one hour of simulation at worst wastes one hour of the computer life, not of yours!

Edited Answer

Thank you for providing the tables. They show that the spread is much wider for the corresponding coefficient when using the variable history than when using the variable test. However, I cannot tell from those tables whether or not this is due to more uncertainty in the posterior: are both variables normalised in the same way? If they are, then indeed the posterior distribution is less precise about the coefficient. Which does not mean you should opt for the model involving test rather than history. This requires model comparison. (I also find curious that the 2.5% posterior quantiles on all coefficients are the sames for Model 1. And the 97.5% posterior quantiles on all coefficients are the sames for Model 2. This hints at some high correlation between the covariates or even lack of identifiability...)