I am doing a Bayesian regression. I have groups of data $(y_1 ~X_1), (y_2~X_2),...$, where each $y$ and $X$ is a vector. The subscript is regarded as group number.

The completely unpooled regression would be $(y_1, y_2,...)~(X_1, X_2,...)$. The completely pooled regression would be to do regression on each $X_i, y_i$ separately.

I am interested in a partial pooling scheme. The thing here is that I believe similar groups, defined by similarity of group number (subscript), will have similar regression coefficients.

Say, $\beta_1$ is a coefficient vector for $y_1 \sim \beta_1^TX_1$, $\beta_2$ is a coefficient for $y_2 \sim \beta_2^TX_2$ and so on. $\beta_{1j}$ is the $jth$ element of vector $\beta_{1j}$.

What I do is to build a covariance matrix using something like an RBF kernel so that $$ K= \begin{bmatrix} k(1,1) & k(1,2) & \dots & k(1,n) \\ k(2,1) & k(2,2) & \dots & k(2,n) \\ \vdots \\ k(n,1) & k(n,2) & \dots & k(n,n) \end{bmatrix} $$, where $k(.,.)$ is the kernel function.

Then I can say that $$ \begin{bmatrix} \beta_{1i} \\ \beta_{2i} \\ \vdots \\ \beta_{ni} \end{bmatrix} \sim N(0, K\sigma^2) $$, where $N(.,.)$ is a multivariate normal distribution. This is the prior for regression coefficiets.

The problem is that I have about 100 groups, but when I try to fit the model using MCMC (pymc3), this is extremely slow.

I was thinking of using a random walk to express the slow change of regression coefficients as we go from one group to the other. The problem don't like this approach is that the variance of beta coefficients of groups with larger group number is much larger.

I wonder if someone has some ideas to do this kind of thing in practice.

  • 1
    $\begingroup$ How did you initialize your problem? Perhaps using a MAP solution (if possible) can then make convergence faster. $\endgroup$ – Vladislavs Dovgalecs Apr 27 '16 at 6:37
  • $\begingroup$ Also, in PyMC3 the NUTS algorithm is available. If your variables are real-valued, may be the convergence can be achieved faster than with MCMC. $\endgroup$ – Vladislavs Dovgalecs Apr 27 '16 at 6:40
  • $\begingroup$ I tried that, but the optimizer cannot even find a MAP solution. There is some numerical instability. If I reduce the number of variables to under 10, then it works. I suspect I have too many variables to optimize and do MCMC. Not sure if 1000 variables is too much for a hierarchical regression - each regression has 10 variables and I have 100 groups. So a total of 1,000 variables. I am using NUTS. $\endgroup$ – Tom Bennett Apr 27 '16 at 6:41
  • $\begingroup$ About those numerical instabilities - do you have data for each individual regression? Is observed data at least 10 points long? May be the MAP problem is ill-posed. $\endgroup$ – Vladislavs Dovgalecs Apr 27 '16 at 6:46
  • $\begingroup$ What if you use the posterior for the solution with 10 variables as a crude prior for your model with 1,000 variables? $\endgroup$ – Vladislavs Dovgalecs Apr 27 '16 at 6:53

Your Answer

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

Browse other questions tagged or ask your own question.