Quantile Regression

[Brendan]

Dear List,

Quantile regression has come up once or twice before, but I can't seem to find a reproducible example of anyone who's got it working.

The only MCMC code I can seem to find online is this example from a JAGS model, based on a data augmentation procedure described by Reed & Yu (2009 )[1].

I've converted it to Stan code (.Rmd file attached, and on GitHub, here, html output here), and compared the results to the JAGS output, and results from the bayesQR and quantreg R packages.

Unfortunately I'm getting whacky results out of Stan. I can't seem to get the model to converge, and the coefficients look pretty far off. My questions:

1. Can anyone spot a problem in my Stan code? Some of the warning messages are unexpected. For example, I'm getting warnings that the exponential distribution function is being fed negative values, which shouldn't be possible (tau is constrained to be positive, as far as I've understood my own code).

2. Could it be that this particular algorithm is only suitable for Gibbs sampling? (From skimming the paper, it seems like this may be the case, however, I'm no expert!)

3. If this quantile regression algorithm is peculiar to GIbbs sampling, what procedure, if any, might be appropriate for Stan?

I'm keen to use Stan (as opposed to falling back to JAGS/Gibbs), as the final model I'd like to fit is QR on-top-of a large hierarchical model.

A minimal reproducible example of quantile regression would be an excellent addition to the Stan list/manual!

Best,

Brendan

Reed, C., & Yu, K. (2009). A Partially collapsed Gibbs sampler for Bayesian quantile regression. 

URL: http://bura.brunel.ac.uk/bitstream/2438/3593/1/fulltext.pdf

[Brendan]

In case Rmarkdown is a bit fiddly, here are the model specs:

JAGS:

  model{
    for(i in 1:length(y)){
      mu[i] <- alpha + beta * x[i]
      w[i]  ~ dexp(tau)
      me[i] <- (1 - 2 * p) / (p * (1 - p)) * w[i] + mu[i]
      pe[i] <- (p * (1 - p) * tau) / (2 * w[i])
      y[i]  ~ dnorm(me[i], pe[i])
    }

    # Regression Priors
    alpha ~ dnorm(0, 1E-6)
    beta  ~ dnorm(0, 1E-6)

    lsigma ~ dunif(-5, 15)
    sigma  <- exp(lsigma / 2)
    tau    <- pow(sigma, -2)
  }

Stan:

data {
  int<lower=0> N;
  real<lower=0, upper=1> p;
  vector[N] x;
  vector[N] y;
}
parameters {
  real alpha;
  real beta;
  real<lower=0.001, upper=10> sigma;
  vector[N] w;
}
transformed parameters {
  real<lower=0> tau;
  vector[N] mu;
  vector[N] me;
  vector[N] pe;

  tau <- pow(sigma, -2);
  mu  <- alpha + beta * x;
  me  <- (1 - 2 * p) / (p * (1 - p)) * w + mu;
  pe  <- (p * (1 - p) * tau) ./ (2 * w);
}
model {
  # Priors
  alpha ~ normal(0, 100);
  beta  ~ normal(0, 100);
  sigma ~ cauchy(0, 2.5);

  # Data Augmentation
  w     ~ exponential(tau);

  # The model
  y     ~ normal(me, pe);
}

[Jonah]

Hi Brendan, 

Looking at the stan code there are several things that would lead to some problems. 

1) Constraints need to be declared for everything that has constraints in parameters and transformed parameters. For example, if w has an exponential distribution then it needs to be positive so you need 

vector<lower=0>[N] w;

and same for pe in transformed parameters, since it's also constrained to be positive. 

2) We don't recommend hard constraints like lower=0.001, upper=10 for scale parameters. Better to declare sigma just to have lower=0 and then use a soft constraint imposed by the prior. 

3) Also regarding pe, it looks like you're treating it the same in the Jags code as in the Stan code, but for Stan the normal distribution is parameterized in terms of the standard deviation whereas (if I remember correctly) Jags uses the precision.

4) We also don't recommend priors like normal(0,100) unless that big of a scale actually makes sense. Here you'd be better off using a much smaller scale for the normal priors on alpha and beta, e.g. normal(0,5) or even normal(0,1). 

With those changes let us know if Stan gives better answers! Hope that helps,

Jonah

[Brendan]

Thanks for the quick response Jonah, and sorry about making such goofy mistakes! It's been a while since I've worked with Stan.

Your suggestions worked a treat. 5k iterations now bring the model to convergence in less than 30 seconds.

Thanks again!

Brendan

P.S. For anyone searching for a reference QR model, the altered model spec is below, and a self-contained dummy-data example is attached.

data {
  int<lower=0> N;
  real<lower=0, upper=1> p;
  vector[N] x;
  vector[N] y;
}
parameters {
  real alpha;
  real beta;

  real<lower=0> sigma;
  vector<lower=0>[N] w;

}
transformed parameters {
  real<lower=0> tau;
  vector[N] mu;
  vector[N] me;
  vector[N] pe;

  vector[N] pe2;

  tau <- pow(sigma, -2);
  mu  <- alpha + beta * x;
  me  <- (1 - 2 * p) / (p * (1 - p)) * w + mu;
  pe  <- (p * (1 - p) * tau) ./ (2 * w);

 
  for(n in 1:N)
    pe2[n] <- inv_sqrt(pe[n]);
}
model {
  # Priors
  alpha ~ normal(0, 2);
  beta  ~ normal(0, 2);

  sigma ~ cauchy(0, 2.5);

  # Data Augmentation
  w     ~ exponential(tau);
 
  # The model

  y     ~ normal(me, pe2);
}

[Jonah Sol Gabry]

Hi Brendan, 

Glad it's working. Just out of curiosity, did you try it for fewer iterations, find that it didn't converge, and then increase to 5000? Or did you start at 5000? It's certainly possible that it needs that many, but my hunch is that you could actually run it for a lot less than that. 

Jonah

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/iw5u9wYK2hI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.