Let $\bf{y}$ denote a vector of binary data, such as whether a group of individuals suffer from a particular disease, and let $\bf{X}$ denote a matrix of potential predictors, including an intercept term. A common approach to this kind of problem is logistic regression, which aims to determine the effect of each predictor in $\bf{X}$ on the log odds of developing the disease.

Within a Bayesian framework, we could write out logistic regression like so

$$ \boldsymbol{\beta}\sim\mathrm{Normal}\left(0,\sigma^2\right)\\ \boldsymbol{\pi}=\mathrm{logistic}\left(\bf{X}\boldsymbol{\beta}\right)\\ \bf{y}\stackrel{\textrm{obs}}{\sim}\mathrm{Bernoulli}\left(\boldsymbol{\pi}\right)\textrm{,} $$ where $\sigma^2$ is a suitable prior variance on the regression coefficients, $\boldsymbol{\beta}$ is the vector of regression coefficients, and $\boldsymbol{\pi}$ is the vector of log odds. Such a model is quite easy to code in modern probabilistic programming packages such as Stan or PyMC3.

However in applied research, in particular when modeling diseases, a liability threshold model is sometimes used instead. I would like to construct a Bayesian version of such a model. As far as I understand, a linear model is placed on a latent variable which, when some threshold is exceeded, causes the disease. The latent variable is usually assumed to be normal with unit variance and the threshold is usually 0.

Here is my attempt to write out a liability threshold model:

$$ \boldsymbol{\beta}\sim\mathrm{Normal}\left(0,\sigma^2\right)\\ \boldsymbol{\psi}\sim\mathrm{Normal}\left(\mathbf{X}\boldsymbol{\beta},1\right)\\ \bf{y}\begin{cases} 0 \textrm{ if }\boldsymbol{\psi}\le0 \\ 1 \textrm{ if }\boldsymbol{\psi}>0 \end{cases}$$

My questions are:

  1. What form does the likelihood of $\bf{y}$ take? It is not normal, because $\mathbf{y} \in [0, 1]$, and it is not Bernoulli, because the values of $\boldsymbol{\psi}$ guarantee an outcome of 0 or 1, rather than giving a probability.
  2. How does one actually code this model in Stan/PyMC?

If the threshold is zero and you integrate $\psi$ out of the posterior, the second model becomes a probit model. Gibbs samplers actually estimate the $\psi$ parameters, but that is not a good idea (if you can avoid it) for Hamiltonian Monte Carlo samplers such as Stan. A Stan program with a probit log-likelihood would be

data {
  int<lower = 0> N;
  int<lower = 0> K;
  matrix[N, K] X; // includes a column of ones
  int<lower = 0, upper = 1> y[N];
parameters {
  vector[K] beta;
  real<lower = 0> sigma; // not squared
model {
  target += normal_lpdf(beta | 0, sigma);
  target += bernoulli_lpmf(y | Phi(X * beta));
  // prior on sigma
  • $\begingroup$ Can you say why this is not a good idea with HMC? $\endgroup$
    – sammosummo
    Oct 17 '19 at 10:24
  • $\begingroup$ You have to estimate N more parameters, which are constrained to be positive or negative depending on the corresponding element of y. $\endgroup$ Oct 17 '19 at 13:52
  • $\begingroup$ Is the issue simply that it is more expensive or will HMC give you biased or poor sampling? $\endgroup$
    – sammosummo
    Jan 23 '20 at 1:18
  • $\begingroup$ Certainly more expensive but dynamic HMC may give poor sampling, but if so, there should be warnings from Stan. $\endgroup$ Jan 23 '20 at 3:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.