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

Question

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?

Answer 1

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...)