Hierarchical negative binomial model

[Meghna Krishnadas]

Hi,

I am new to Stan and am trying to construct a hierarchical model as per the following requirements:

response Y ~ neg_binomial_2( mu.y , phi.y )  ## the dots are only for explanation purposes here; the parameters are not named with dots in the Stan code.

where, mu.y = exp(a) * X^b  ## This is a biologically meaningful and established relationship for the variables being measured in this system.

and, X is itself drawn from a negative binomial distribution    X ~ neg_binomial_2( lamda , phi.x) with lamda ~ lognormal( lamda.bar , sigma)

lamda.bar = mean(X.obs.plot), where X.obs is a vector of the observed values of the predictor at that location/plot, measured over multiple years. The expected number of X for a plot is lamda and these were drawn from a lognormal distribution with mean lamda.bar and variance sigma sq.

FYI, to give you an idea of why I need to undertake this joint modeling of the probability distribution of Y: 

In this data, the response Y was measured in plots adjacent to where the predictor X.obs was measured. This spatial mismatch implies that the values of Y do not correspond to the true value of X at that location, but have to be estimated from X.obs. Additionally, in many cases, the value of X.obs is 0 but there is a non-zero value for Y. This again occurs due to the spatial mismatch mentioned above. Therefore, the actual values of X used to obtain mu.y for the negative binomial distribution for estimating the response Y, have to be drawn from the values of X.obs. The response Y has to modeled from a joint probability distribution of the expected values for X and Y.

I've pasted my code and error messages below. There might be something very basic in the coding that I am unable to understand and hence going wrong with. Any help will be much appreciated! 

Thank you for your time.

Best,

Meghna

perspp.nb.hier<-"

data{

  int<lower=0> N; //no.of data points

  vector[N] seed; //predictor: positive real number, range= 0 - 50000 (~60% values are zeros)

  int<lower=0> sdl[N]; //response: integer, range = 0 - 212 (~50% values are zeros)

}

parameters{

  real sigma;

  real<lower=0> lamda;

  real<lower=0> phisd;

  real alpha;

  real beta;

  real<lower=0> phisdl;

}

model{

  int seed2[N];

  sigma ~ uniform(0,10);

  phisd ~ uniform(0,10);

  alpha ~ uniform(0,100);

  beta ~ uniform(0,10);

  phisdl ~ uniform(0,10);

  for(i in 1:N){

    lamda[i] ~ lognormal(seed[i], sigma);

    seed2[i] ~ neg_binomial_2(lamda[i], phisd);

    }

    sdl ~ neg_binomial_2_log(exp(alpha)*seed2^beta, phisdl);

  }

}"

DIAGNOSTIC(S) FROM PARSER:

Warning (non-fatal):

Left-hand side of sampling statement (~) may contain a non-linear 

transform of a parameter or local variable.

If so, you need to call increment_log_prob() with the log absolute 

determinant of the Jacobian of the transform.

Left-hand-side of sampling statement:

    seed2[i] ~ neg_binomial_2(...)

SAMPLING FOR MODEL 'eb5d4d92b044bf68050a15a91450a8dc' NOW (CHAIN 1).

 [1] "Rejecting initial value:"                                       

                                                                      

 

 [2] "  Error evaluating the log probability at the initial value."   

                                                                      

 

 [3] "Informational Message: The current Metropolis proposal is about 

to be rejected because of the following issue:"                       

 

 [4] "Exception thrown at line 30: stan::math::neg_binomial_2_log: 

Location parameter is 0, but must be > 0!"                            

    

 [5] "If this warning occurs sporadically, such as for highly 

constrained variable types like covariance matrices, then the sampler 

is fine,"

 [6] "but if this warning occurs often then your model may be either 

severely ill-conditioned or misspecified."                            

  

 [7] "Rejecting initial value:"                                       

                                                                      

 

 [8] "  Error evaluating the log probability at the initial value."   

                                                                      

 

 [9] "Initialization between (-2, 2) failed after 100 attempts. "     

                                                                      

 

[10] " Try specifying initial values, reducing ranges of constrained 

values, or reparameterizing the model."                               

  

[1] "error occurred during calling the sampler; sampling not done"

[Bob Carpenter]

May I ask what you read before you tried to write this
model? The reason I ask is that the two problems in your
code are quite common and we'd like to be able to more
clearly formulate the issues.

Issue 1.

The variables in Stan need to be declared with constraints
so that any value satisfying the constraints has finite
log density. So when you have a parameter declared
as

> real sigma;

and a sampling statement

> sigma ~ uniform(0,10);

you produce a problem where sigma < 0 or sigma > 10 satsify
the declared constraints (there aren't any).

If you insist on this prior, you need to declare

real<lower=0, upper=10> sigma;

and similarly for the other parameters.

Issue 2.

You are using the local variable seed2 before it is defined:

int seed2[N];
...

seed2[i] ~ neg_binomial_2(lamda[i], phisd);

Stan doesn't work like BUGS --- this just increments the posterior
log density. In BUGS, it defines an edge in a directed graphical
model.

I think the only solution is to marginalize seed2 out of your
model. I'm not sure how a compound negative binomial works. You're
giving your model a double-dose of overdispersion.

Issue 3.

You can't raise an array to a power, so seed2^beta is going to
be ill-formed as well. The legal signatures are defined in the
index to the manual.

- Bob

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

[Meghna Krishnadas]

Hi Bob,

Thanks for your response and for outlining the major issues with my model. I apologize for my delay in reverting to you - I am undertaking fieldwork in a remote location in India and have very limited connectivity of any sort. 

​To 

answer

​​

your questions:

1. I used the Stan manual and online examples for constructing negative binomial models

​.​

2. The idea to jointly model the probability distributions comes from a publication that implemented this method using INLA in R.

​I have a couple of more questions in ​response to your suggestions:

​1. Can I define seed2 as a vector instead and use the ​generated quantities option to determine the values of seed2 from the distribution that describes it? I couldn't determine a way to draw random values from a prescribed distribution in Stan, as can be done in R.

2. Or, could I put the main response within the loop that computes seed2, writing it as:

​

 for(i in 1:N){
>     lamda[i] ~ lognormal(seed[i], sigma);
>     seed2[i] ~ neg_binomial_2(lamda[i], phisd);

>     sdl[i] ~ neg_binomial_2_log(exp(alpha)*seed2[i]^beta, phisdl);

       }

I had previously written it as this, but that was throwing errors too. 

Thanks,

Meghna

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------

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

[Bob Carpenter]

I don't quite understand (2) --- it's much easier for us if
we have the entire model to look at. You can probably vectorize
that code to be much more efficient if lambda, seed2, and sdl are
1D containers (vectors or 1D arrays).

In this case, for Stan, seed2 and sdl would need to be data.
Stan can't sample discrete parameters.

- Bob