Sum of neg_binomial_2, what about the parameters?

[Ste]

Hello Community,

there are some publication that address this question for other parametrizations.

For a sum of neg_binomial_2 is a another neg_binomial_2 with 

- mu_tot = sum(mu)

- phi_tot = ?

What is the resulting resulting phi? Is what function of phi[n]?

Thanks a lot! 

[Marco Inacio]

I tryed the lazy way, and was not able to find a relation between the two when we change both \phi and \mu:

phi1 = sample(c((1000:1)/100, (1:100)/100), 200, TRUE)
phi2 = sample(c((1000:1)/100, (1:100)/100), 200, TRUE)

phi1 = c(phi1, runif(1000, 0.1, 100))
phi2 = c(phi2, runif(1000, 0.1, 100))

var_ = mapply(function(phi1i, phi2i) { var(c(rnbinom(1e4, phi1i, mu=1) + rnbinom(1e4, phi2i, mu=3))) }, phi1, phi2)

phi_new <- 4**2/(var_-4)

> summary(lm(phi_new ~ phi1 + phi2 + I(phi1^2) + I(phi2^2) + phi1*phi2 - 1))

Call:
lm(formula = phi_new ~ phi1 + phi2 + I(phi1^2) + I(phi2^2) +
    phi1 * phi2 - 1)

Residuals:
   Min     1Q Median     3Q    Max
-556.3  -29.3    2.4   11.8 8128.4

Coefficients:
           Estimate Std. Error t value Pr(>|t|)   
phi1       0.454870   0.916078   0.497 0.619604   
phi2      -0.429264   0.938322  -0.457 0.647409   
I(phi1^2) -0.008243   0.010014  -0.823 0.410618   
I(phi2^2)  0.012547   0.010302   1.218 0.223500   
phi1:phi2  0.030667   0.008202   3.739 0.000193 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 272.8 on 1195 degrees of freedom
Multiple R-squared:  0.1567,    Adjusted R-squared:  0.1532
F-statistic: 44.42 on 5 and 1195 DF,  p-value: < 2.2e-16

This math exchange answer https://math.stackexchange.com/questions/1054048/negative-binomial-distribution-sum-of-two-random-variables tells us that for the special case where \phi (remember that \phi \equiv r) does change but p doesn't, then we have another neg_binomial for the sum of such iid random variables.

Unfortunately, this is not the case in general when you change \phi, but keeps \mu. There might be other special cases, but if we follow the one in such answer, then p=\mu / (\mu+\phi) must be keep constant for the sum to be another negative binomial.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, 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.

[stefano]

Thanks Marco,

thanks for your answer. So a mix of neg_binomial_2 can only be in the form 

y[n] ~ neg_binomial_2( mu_1 + mu_2, phi );

where phi is a unit number, right?

Bw. 

To unsubscribe from this group and stop receiving emails from it, 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.

--
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/MqLQH3jeUCg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Marco Inacio]

Not really, according to that example, it's p that must be constant, where p is the probability of sucess of the usual negative binomial parameterization of undergraduate courses, here I'll call it neg_binomial_undergrad (it doesn't exist on Stan, of course).

Then if you have something like:

y[1] ~ neg_binomial_undergrad(p, phi1)

y[2] ~ neg_binomial_undergrad(p, phi2)

Then:

(y[1] + y[2]) ~ neg_binomial_undergrad(p, phi1 + phi2)

On Stan neg_binomial_2, you'll have to use the fact that:

mu = p * phi/ (1 - p)

Then, you have:

y[1] ~ neg_binomial_2(p * phi1/ (1 - p), phi1)

y[2] ~ neg_binomial_2(p * phi2/ (1 - p), phi2)

(y[1] + y[2]) ~ neg_binomial_2(p * (phi1 + phi2)/ (1 - p), phi1 + phi2)

Not that useful, but as I said there might be other special cases of sum of independent negative binomial random variables that gives another negative binomial that I'm not aware of.

A quick test to check the results on R:

p = .4
phi1 = 11.5
phi2 = 4.6
first_ = rnbinom(1e6, phi1, mu=p * phi1/ (1 - p)) + rnbinom(1e4, phi2, mu=p * phi2/ (1 - p))
second_ = rnbinom(1e6, phi1 + phi2, mu=p * (phi1 + phi2)/ (1 - p))
mean(first_)
mean(second_)
var(first_)
var(second_)
mean(first_**3)
mean(second_**3)
quantile(first_)
quantile(second_)

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Marco Inacio]

Well, it's good to mention, if don't already know, that there's an import special case that might be useful, it's when they independent and identically distributed:

y[1] ~ neg_binomial_undergrad(p, phi)

y[2] ~ neg_binomial_undergrad(p, phi)

(y[1] + y[2]) ~ neg_binomial_undergrad(p, 2*phi)

On Stan we would have:

y[1] ~ neg_binomial_2(mu, phi)

y[2] ~ neg_binomial_2(mu, phi)

(y[1] + y[2]) ~ neg_binomial_2(2*mu, 2*phi)

mu = 4.4
phi = 11.5
first_ = rnbinom(1e6, phi, mu=mu) + rnbinom(1e6, phi, mu=mu)
second_ = rnbinom(1e6, 2*phi, mu=2*mu)

mean(first_)
mean(second_)
var(first_)
var(second_)
mean(first_**3)
mean(second_**3)
quantile(first_)
quantile(second_)

> mean(first_)
[1] 8.803444
> mean(second_)
[1] 8.80471
> var(first_)
[1] 12.15042
> var(second_)
[1] 12.19663
> mean(first_**3)
[1] 1024.479
> mean(second_**3)
[1] 1026.228
> quantile(first_)
  0%  25%  50%  75% 100%
   0    6    9   11   32
> quantile(second_)
  0%  25%  50%  75% 100%
   0    6    9   11   31

[stefano]

Thanks a lot.

I guess the best compromise is to assume that all neg_binomial_2 in the sum have same phi in the way to have

(y[1] + y[2]) ~ neg_binomial_2(mu_1 + mu_2, 2*phi);

Thanks your answers were really helpful.

Stefano

Stefano Mangiola

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

Hovens lab - Papenfuss lab

Royal Melbourne Hospital - WEHI

+61 (0)450135954

[Marco Inacio]

No, this will be true only in the limit of phi going to infinity, where both random variable will be Poisson distributed and therefore their sum will also be Poisson.

See for instance when phi=0.1

set.seed(12)
phi = 0.1
mu1 = 1
mu2 = 10
first_ = rnbinom(1e6, phi, mu=mu1) + rnbinom(1e6, phi, mu=mu2)
second_ = rnbinom(1e6, 2*phi, mu=mu1+mu2)

mean(first_)
mean(second_)
var(first_)
var(second_)
mean(first_**3)
mean(second_**3)
quantile(first_)
quantile(second_)


> mean(first_)

[1] 11.03025
> mean(second_)
[1] 10.98656
> var(first_)
[1] 1027.162
> var(second_)
[1] 613.6651
> mean(first_**3)
[1] 240507.6
> mean(second_**3)
[1] 89080.94

> quantile(first_)
  0%  25%  50%  75% 100%

   0    0    1    6  902

> quantile(second_)
  0%  25%  50%  75% 100%

   0    0    1   10  653



Now look how similar they are when phi=100

set.seed(12)
phi = 100
mu1 = 1
mu2 = 10
first_ = rnbinom(1e6, phi, mu=mu1) + rnbinom(1e6, phi, mu=mu2)
second_ = rnbinom(1e6, 2*phi, mu=mu1+mu2)

mean(first_)
mean(second_)
var(first_)
var(second_)
mean(first_**3)
mean(second_**3)
quantile(first_)
quantile(second_)



> mean(first_)

[1] 10.99724
> mean(second_)
[1] 10.99785
> var(first_)
[1] 12.02925
> var(second_)
[1] 11.62556
> mean(first_**3)
[1] 1741.239
> mean(second_**3)
[1] 1726.74

> quantile(first_)
  0%  25%  50%  75% 100%

   0    9   11   13   32

> quantile(second_)
  0%  25%  50%  75% 100%

   0    9   11   13   33

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[stefano]

Thanks Marco,

1) but then what is the best approach to model sums of neg_binomials

2) I am still trying to understand how good is this approximation

Let's say that in median I have 800 as count value, and empirically phy = 2, for those parameters

set.seed(12)
phi = 2
mu1 = 50
mu2 = 800

first_ = rnbinom(1e6, phi, mu=mu1) + rnbinom(1e6, phi, mu=mu2)
second_ = rnbinom(1e6, 2*phi, mu=mu1+mu2)
mean(first_)
mean(second_)
var(first_)
var(second_)

quantile(first_)
quantile(second_)

> mean(first_)
[1] 850.4116
> mean(second_)
[1] 849.0126
> var(first_)
[1] 323373.7
> var(second_)
[1] 180821

P.S. Why do you analyse mean(first_**3) and mean(second_**3)?

With low means

set.seed(12)
phi = 2
mu1 = 5
mu2 = 20

first_ = rnbinom(1e6, phi, mu=mu1) + rnbinom(1e6, phi, mu=mu2)
second_ = rnbinom(1e6, 2*phi, mu=mu1+mu2)
mean(first_)
mean(second_)

var(first_)
var(second_)

> mean(first_)
[1] 25.01095
> mean(second_)
[1] 24.98901
> var(first_)
[1] 237.6719
> var(second_)
[1] 180.9081

Now what do you think about this is an imperfect assumption but could still be acceptable, or is wrong and that's it. In this case, how can I model sums of neg_binomials?

Thanks!

[Marco Inacio]

On 07/04/2017 00:54, stefano wrote:

Thanks Marco,

1) but then what is the best approach to model sums of neg_binomials

It will depend on your assumptions and on what you have.

Do you observe the variables individually? If so, just model them individually on Stan.

Or do you only observe their sum and know that they are individually negative binomial distributed? If that's the case, maybe you could model the individual components as parameters in Stan, constrained so that their sum is your observed sum and attempt to place all the information you have on the relation among them as prior distributions.

2) I am still trying to understand how good is this approximation

Yeah, obviously, the lower the value of phi, the worse should be the approximation is since we are getting further apart from the Poisson case and also the further apart mu1 and mu2 are from each other, the worse it should since you are getting further apart from the special case of i.i.d random variables. But how worse it gets it possibly pretty hard question.

Let's say that in median I have 800 as count value, and empirically phy = 2, for those parameters

set.seed(12)
phi = 2
mu1 = 50
mu2 = 800
first_ = rnbinom(1e6, phi, mu=mu1) + rnbinom(1e6, phi, mu=mu2)
second_ = rnbinom(1e6, 2*phi, mu=mu1+mu2)
mean(first_)
mean(second_)
var(first_)
var(second_)
quantile(first_)
quantile(second_)

> mean(first_)
[1] 850.4116
> mean(second_)
[1] 849.0126
> var(first_)
[1] 323373.7
> var(second_)
[1] 180821

P.S. Why do you analyse mean(first_**3) and mean(second_**3)?

The third moment and the quantiles are just an extra to make "sure" the random variable have the same or different distribution.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[S]

On Friday, April 7, 2017 at 2:43:33 PM UTC+10, Marco Inacio wrote:

On 07/04/2017 00:54, stefano wrote:

Thanks Marco,

1) but then what is the best approach to model sums of neg_binomials

It will depend on your assumptions and on what you have.

Do you observe the variables individually? If so, just model them individually on Stan.

Or do you only observe their sum and know that they are individually negative binomial distributed? If that's the case, maybe you could model the individual components as parameters in Stan, constrained so that their sum is your observed sum and attempt to place all the information you have on the relation among them as prior distributions.

Indeed this is the case. Your solution can be a good one. The aspect I want to observe is if I assume they all have the same dispersion (even though they wouldn't necessarily) how the accuracy of the model drop. This rough approximation could make things a lot easier and maybe is good enough.

 

Bw. 

To unsubscribe from this group and stop receiving emails from it, 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.

--
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/MqLQH3jeUCg/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.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, 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.

--
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/MqLQH3jeUCg/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.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, 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.

--
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/MqLQH3jeUCg/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.