prior over NB regression predicted points - jacobian

[Ste]

Hello Community,

I have a negative binomial regression with log link, and I would like to put a prior on the predicted values

"e" is the predicted regression values

 

transformed parameters {

  vector<lower=0>[S] e[G,P];//the linear predictor

  for(g in 1:G) for(p in 1:P) e[g,p] = exp(X*gamma[g,p]); //using the log link 

}

"gamma" is the regression parameters and does NOT have prior on it

I would like to put

for(p in 1:P) for(g in 1:G) e[g, p] ~ gamma(alpha[p,g], beta[p,g]);

The question is do I have to add the Jacobian adjustment?

target += gamma * X

Thanks a lot.

[Ben B]

Yeah, I think you need a Jacobian here. It's Chapters 20 (section "Jacobians not Necessary" and 20.3) and 33 in the manual that cover this.

I'm assuming X is a constant [S, K] matrix, and each gamma[g, p] is a length K vector.

If K > 1, I'm not super sold on the adjustment being just a bunch of X * gammas (not clear from the model if it is or isn't).

I've just been thinking through this stuff recently myself, so take my word with a grain of salt, but, if I'm thinking clearly:

If things were scalars, we want to sample on gamma, and we have specified a distribution on exp(X * gamma), then f^-1 = exp(X * gamma), df^-1/dgamma = X * exp(X * gamma), log(abs(df^-1/dgamma)) = log(abs(X)) + X * gamma, and you could drop the log(abs(X)) because it is constant.

If you're dealing with matrices/vectors, you have a multivariate transform for each exp(X * gamma[g, p]). It's not obvious to me your Jacobian would in general be triangular or diagonal (because of the matrix product X * gamma[g, p]) so check that (Chaper 33, "Multivariate Changes of Variables").

Does that make sense at all?

Ben

[andre.p...@googlemail.com]

Ben B,

for those who have access to Mathematica:

(Notes: I replaced gamma with g, to avoid confusion with the Gamma function, Re[g x] = real of g * x)

PowerExpand[Log[Abs[D[Exp[x*g], g]]]]

returns: Log[Abs[x]] + Re[g x]

Thus you calculation is right.

Andre

[Ben B]

Thanks for the spot check Andre!

Ben

[stefano]

Thanks guys, I have read the manual already but still I can't grasp the concept for vectors.

let me see if I got it.

if X and g were scalar (one real number each) then my jacobian (target += g* X) is right.

In my case the realistic scenario is, for example

X = [1,  0

        1,  0

        1,  0

        1,  0

        1,  1

        1,  1

        1,  1

        1,  1]

OR

X = [1,  1.5

        1,  2.0

        1,  1.0

        1,  0.5

        1,  5.4

        1,  5.5

        1,  6.6

        1,  2.1]

and 

g = [2.3,  5.0] (intercept and slope respectively)   (there will be many regressions, many g's, in the model but i guess it's not important here)

(in this case target += g* X would result into target += g0 + g1* X1)

Is this correct in my case? Anyway target += g* X provide a model that does not converge.

This is the model without prior for comparison

I also attach the full model with data

====================================================================

model_code="

data {

int<lower=1> N;                                           // number of regressions

int<lower=1> S;                                           // number of points for regression

  int y[N,S];                                                            // observed

  int K;                                                                    // k = 2

  matrix[S,K] X;                                                     // the model matrix

}

parameters {

vector<lower=0>[N] sigma;                       // overdispersion

real<lower=0> sigma_mu;                        // prior overdispersion

real<lower=0> sigma_sd;                         // prior overdispersion

vector[K] g[N];                                              // the regression parameters

vector<lower=1>[N] alpha;                       // prior predicted points

vector<lower=0>[N] beta;                         // prior predicted points

real<lower=1> alpha_alpha;                   // Hyper prior predicted points

real<lower=0> alpha_beta;                     // Hyper prior predicted points

real<lower=1> beta_alpha;                     // Hyper prior predicted points

real<lower=0> beta_beta;                       // Hyper prior predicted points

}

transformed parameters {

  vector<lower=0>[S] e[N];                                //the linear predictor

  for(n in 1:N) e[n] = exp(X*g[n]);                     //using the log link 

}

model {

//-----------------------------------------------------//

//--------------REGRESSION--------------------//

//-----------------------------------------------------//

for(n in 1:N) for(s in 1:S) y[n,s] ~ neg_binomial_2( e[n, s], sigma[n]);

sigma ~ normal(sigma_mu, sigma_sd);

//-----------------------------------------------------//

//--------------JACOBIAN-------------------------//

//-----------------------------------------------------//

for(n in 1:N) e[n] ~ gamma(alpha[n], beta[n]);

for(n in 1:N) target += X*g[n];

// //-----------------------------------------------------//

// //-----------------------------------------------------//

alpha ~ gamma(alpha_alpha,alpha_beta);

beta ~ gamma(beta_alpha,beta_beta);

}

"

fit = rstan:::sampling(

model_code=model_code, 

data= data,

iter=                             500

)

====================================================================

Thanks for the support!

2017-04-24 3:51 GMT+10:00 Ben B <bbba...@gmail.com>:

Thanks for the spot check Andre!

Ben

--
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/HaWznIobD2w/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[andre.p...@googlemail.com]

Ste,

your target makes sense. If you code it in matrix form. It's easy to add. You don't specify priors on some of your hyperparameters.

This is just a 2 cents statement: I coded the lindley-poisson, some bivariate poisson combound from Steinitz, the poisson shanker,

etc. and I marginalized them out, but I never see in those papers, some transformation. So I question, the NB is already a compound gamma 

with poisson, so why use the gamma again? And where does the idea come from? What is the math behind?

Andre

[stefano]

Hello Andre,

thanks for your message. I will try to add hyper-hyper parameters on alpha_alpha, alpha_beta, beta_alpha, and beta_beta. 

- Do you confirm that it should definitely work? 

- What do you mean with " If you code it in matrix form. It's easy to add." could you explicit the code that I could use?

========== Reason ===========================================================================

The reason of hyper-prior on "e" is more logical than mathematical. But since it worked beautifully for categorical grouping instead of regression it should work for regression as well.

I explain: instead of a regression with 

X = [1,  0

        1,  0

        1,  0

        1,  0

        1,  1

        1,  1

        1,  1

        1,  1]

I can structure the model separating X[,2] = 0 and X[,2] = 1 (in this simple case regression and grouping is exactly the same)

1. I calculate the mu for both groups

y[g,s] ~ neg_binomial_2( e[COVARIATE_1_o_2], sigma);

2. I put the two e[1] and e[2] under a Gamma prior 

for(i in 1:2) ex_mean[i] ~ gamma(alpha, beta);

3. I force alphas and betas (many for many set of points, many "regressions") to share priors

for(p in 1:P) alpha[p] ~ gamma(alpha_alpha,alpha_beta);

for(p in 1:P) beta[p] ~ gamma(beta_alpha,beta_beta);

This results is a clear clustering of mode and sd of the gamma distributions, in which the background (mode/sd) is clearly visible and those regressions that are interestingly altered (more variable) are clearly clustered differentially

That's why I want to put Gamma prior not on two mu's but on the regression predicted values. It is conceptually the same

===========================================================================================

 

Someone can give me more information why I get missing mixing of chains in my case?

Thanks

Stefano Mangiola

Data analyst (MSc, MPhil, PhD candidate)
@ Melbourne University

Hovens lab - Papenfuss lab

Royal Melbourne Hospital - WEHI

+61 (0)450135954

[Ben B]

This looks like a glm of some sort. It looks like you're using R, so have you see Rstanarm? Check this out: https://cran.r-project.org/web/packages/rstanarm/vignettes/count.html

There's a section on fitting count data with negative binomials there as well (https://cran.r-project.org/web/packages/rstanarm/vignettes/count.html#poisson-and-negative-binomial-regression-example). Maybe this package will let you do your regression more easily?

But back to your model, I think the distributions on 'e' you're specifying are in a weird place. I think I agree with Andre that you've got an extra gamma distribution in there that might not need to be there. Check out: https://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma.E2.80.93Poisson_mixture . Can the gamma in this case just get integrated out? I think doing that transform specifies both parameters of the negative Poisson as well (so there's no free sigma).

Now to the Jacobians, I think the Jacobians not making sense is a result of the model.

For example, just taking three rows of your X:

X = [1, 2.0

1, 1.0

1, 0.5]

If y = Xg

y0 = g0 + 2 * g1

y1 = g0 + g1

y2 = g0 + 0.5 * g1

This means

y0 = 3 * y1 - 2 * y2

So if you specify distributions on y1 and y2, you already have a distribution on y0!

So applying the priors where you are I think is creating issues cause there are bunch of different conflicting soft constraints. I don't deeply understand the Math behind when you can easily apply the Jacobian adjustment or not, but it's not always.

Ben

[Bob Carpenter]

> On Apr 23, 2017, at 5:41 AM, Ste <mangiol...@gmail.com> wrote:
>
>
> Hello Community,
>
> I have a negative binomial regression with log link, and I would like to put a prior on the predicted values

Not a good idea as it breaks the generative nature of the model.

> "e" is the predicted regression values
>
> transformed parameters {
> vector<lower=0>[S] e[G,P];//the linear predictor
> for(g in 1:G) for(p in 1:P) e[g,p] = exp(X*gamma[g,p]); //using the log link
> }
>
> "gamma" is the regression parameters and does NOT have prior on it
>
> I would like to put
>
> for(p in 1:P) for(g in 1:G) e[g, p] ~ gamma(alpha[p,g], beta[p,g]);

No matter how you try to put a distribution on exp(X * gamma[g, p]), it
won't be a proper Bayesian prior because it involves data. Thus it breaks
the interpretation of a Bayesian model as generating parameters from the
prior and then generating data from the parameters.

If you absolutely must include constraints (usually they're a bad
idea for both statistical and computational reasons), you're better
off putting them on the gamma parameter so that the model at least
remains generative (that is, you can generate the data from the
parameters without rejection sampling!).

Or you should include them in some kind of link function.

>
> The question is do I have to add the Jacobian adjustment?
>
> target += gamma * X

If you want to put a distribution on e, then you indeed have to
apply the Jacobian. If X is a vector, you'll need to do this in
loop with X * gamma[g, p].

- Bob

[stefano]

Thanks Guys,

I think Bob, that indeed while the regression is my generative process, the gamma observation is generated and not generating. This discussion is getting interesting, let me show you what it does work, with tested data, outperforming current methods. The framework is not regression but grouping. Then maybe we can generalise it to continuous.

Let's suppose I have many worlds (W). In each of those, old and young people [P] (two discrete groups) have a "distinct" OR "same" number of hair ~ NB(). This could be modelled in this way:

Let young and old hair be generated by two different NB 

hair_young[p,w] ~ NB(mu_young[w], phi[w]);

hair_young[p,w] ~ NB(mu_old[w], phi[w]);

Let the two estimates mu be generated by a common Gamma

mu_young[w] ~ gamma(alpha[w], beta[w]);

mu_old[w] ~ gamma(alpha[w], beta[w]);

Now this does the trick, let those gammas for all words be generated by the same process

alpha[w] ~ gamma(a, b);

beta[w] ~ gamma(c, d);

This allows to tell apart the worlds where young and old have different number of hair. Yes because the mere difference mu_young[w]  - mu_old[w] and maybe its uncertainty (like in a frequentist approach) does NOT separate the two cases (i.e., same/  different).

But alphas and betas DO! That was for me the revelation. Because 

-- gamma generated by a age invariant process (regression y = a) in which there is no change follow a precise background alpha/beta ratio,

-- gammas generated by an age variant process  (regression y = a1 + a2 * x) have a really distinct alpha/beta ratio, having bigger variance.

Figure: mode/sd of gammas; red being age dependent

This is how I solved the hypothesis testing problem in the Bayes context, for two categories. Now this underlie a continuous generalization, but I could not crack it. What about if age is continuous and not discrete?

Thanks to all!

[stefano]

Eureka!

I found what was my misconception was in my working model. It is working but is not doing what I was thinking :)

mu_young[w] ~ gamma(alpha[w], beta[w]);

mu_old[w] ~ gamma(alpha[w], beta[w]);

is indeed not sampling the transformed (implicitly) values

mu_young = g0

mu_old = exp(g0 + g1 * x) = exp(g0) * exp(g1 * x)

but is indirectly sampling

exp(g0) ~ gamma(a1, b1);

exp(g0) * exp(g1 * x) ~ gamma(a2, b2);

therefore

exp(g1 * x) ~ gamma(a2, b2) / gamma(a1, b1);

therefore

exp(g1 * x) ~ beta_prime() OR ~ f_distribution()

references

https://stats.stackexchange.com/questions/207595/distribution-of-the-ratio-of-two-gamma-random-variables

https://stats.stackexchange.com/questions/143311/distribution-of-the-quotient-of-two-gamma-random-variables-with-different-rate-p

I will soon test this hypothesis and letting you know ;) This would have been a great deductive/reverse-engineering/phylosophical victory!

[andre.p...@googlemail.com]

I did one time the gamma count distribution in Stan:

real gammacount(int y , real la, real a) {

    if(y>0)

      return gamma_p(y*a, a*la) - gamma_p((y+1)*a, a*la);

    return gamma_q((y+1)*a, a*la);

  }

target += log(y, exp(mu[i]), lambda[1]);

It uses the difference of the gamma CDF to get the PMF. 

Not tested: If you sample mu.nb ~ gamma(?!?!, a*la) you get the PDF

with la = exp(mu). Note la > 0, a > 0, a ~ gamma(b,c). 

?!?! left as homework. (Sorry I don't have time to dig into this now,

but I'm sure its not complicated.)

Andre

[Ben B]

I don't know the full story here, but I think this is some sort of hierarchical regression problem.

You have collected hair count data for a number of different people of different ages (regression part, difference in two ages vs. continuum of ages shouldn't be that big). You have done this experiment in many different worlds, but you expect the relations to be similar (hierarchical part).

I'd start by just building that regression with Rstanarm and going from there (https://cran.r-project.org/web/packages/rstanarm/vignettes/count.html, https://cran.r-project.org/web/packages/rstanarm/vignettes/glmer.html).

It might be simpler to use the posteriors of the regression coefficients to analyzer your data than trying to figure out some kind of hypothesis test.

Good luck!

[stefano]

Yes, I think all comes down to this.

Can a hierarchical regression be better than a non-hierarchical, in my scenario?

News soon.. you guys were already all really helpful.

--