Multilevel Regression too computationally expensive, probably implemented incorrectly

[Pedro Carvalho]

I first asked this question on Cross Validated, where it probably has better visualisation.

I have, for N_countries countries and for each day during a certain period, a list of purchases from my website and how much each purchase was. I think that money spent as a function of days since the start of this period follows a sigmoid.

I want to model this according to a multilevel model, in which some of the parametres are drawn from distributions whose hyperparametres are drawn from another distribution (my population level is the country). The way I did it is, however, incredibly computationally expensive:

   

data {
      int<lower = 1> N_people_day;
      int<lower = 1> N_countries;
      int<lower = 1, upper = N_countries> country[N_people_day];
      int<lower = 0> dsi[N_people_day];
      vector<lower = 0>[N_people_day] net_revenue;
    }
   
    parameters {
      // parametres
      vector<lower = 0>[N_countries] scale;                       // the sigmoid's scale parametre
      real<lower = 0> rate;                                       // rate parametre
      real<lower = 0> x_location;                                 // location parametre
      vector<lower = 0>[N_countries] y_location;                  // location parametre
      real<lower = 0> sigma;                                      // dispersion parametre
     
      // hyperparametres
      real<lower = 0> y_location_hyper_shape;
      real<lower = 0> y_location_hyper_rate;
      real<lower = 0> scale_hyper_shape;
      real<lower = 0> scale_hyper_rate;
    }
   
    transformed parameters {
      vector[N_people_day] mu;
     
      for(i in 1:N_people_day){
        mu[i] = scale[country[i]] / (1 + exp(-(rate * dsi[i] + x_location))) + y_location[country[i]];
      }
    }
   
    model {
      net_revenue ~ normal(mu, sigma);
      y_location ~ gamma(y_location_hyper_shape, y_location_hyper_rate);
      x_location ~ gamma(10, 0.5);
      scale ~ gamma(scale_hyper_shape, scale_hyper_rate);
    }

I'm almost certain there's a better way to do it and that the biggest problem is that one for in the transformed parameters block. How should I do it instead, though?

[Andrew Gelman]

Hi, just a couple quick comments.  There's something kinda weird about your model because the data are positive but you're fitting a normal model with additive errors.  Can't net_revenue be negative sometimes?  Why must it be positive.  Also in your expression for mu, you can simplify the code by vectorzing (just omit the subscripts [i] and take it out of a loop), also use the inv_logit function if that's what you want.  But I don't think this inv_logit thing makes much sense either.  Also it doesn't look right that y_location is additive with a gamma distribution.  Lots of things about this model look odd.  I recommend you re-think it from first principles, as what you have below looks like a mash-up of distributions with no underlying model.

A

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

[Pedro Carvalho]

Hello!

net_revenue cannot be negative - a purchase can only give me money, it cannot take money from me - and the model assumption is "a purchase on a given day will give me a value of money that follows a logistic function of the day plus random noise."

Trying to omit the subscripts [i] gives me an error: it says there are no matches for real * int[]. If I convert all my int[] to vectors I get no matches for vector / vector, plus "Container index must be integer" at scale[country].

The inv_logit thing is mostly because, empirically, it looks like my data behaves like that: the average is mostly constant for a while, then it grows, then it becomes mostly constant at a new value afterwards.

The model doesn't really have a lot going for it in the way of first principles - what I want to do is just a regression, because it seemed like what made the most sense but I'm not super confident in this.

Pedro

[Andrew Gelman]

Hi, I still don't understand how net revenue can't be negative (usually when something is "net", it's gross minus some correction), but if it really can't be negative then I think you'd want a multiplicative, not random error.  Regaarding the int-to-real conversion:  that I'm not so sure about.  I know Stan has such functions but I always get confused, as I'm used to R, where everything's real and there are no integers.

A

[Pedro Carvalho]

In this case the "net" is gross times 0.7. How would I model the multiplicative, exactly?

And my biggest problem is exactly with the vectorisation part: that for takes forever to run, but even the example linear regression model on git uses a for rather than vectorisation so I don't know if it's possible at all.

[Bob Carpenter]

No, Stan provides no real to integer conversion on purpose.
Its algorithms are aimed at continuously differentiable densities.

> On Dec 16, 2016, at 5:31 PM, Pedro Carvalho <pedro...@gmail.com> wrote:
>
> In this case the "net" is gross times 0.7. How would I model the multiplicative, exactly?

lognormal.

That's an unusual net function. Usually it's income minus costs
and the costs aren't constant proportions of the gross.

> And my biggest problem is exactly with the vectorisation part: that for takes forever to run, but even the example linear regression model on git uses a for rather than vectorisation so I don't know if it's possible at all.

Yes, you can vectorize. We just added a bunch of language features
like general vectorization in 2.13 so we'll sweep through and update
all the models soon, but we probably won't be done until well into next
year.

- Bob

[Pedro Carvalho]

How do you vectorise? Just removing the subscript and taking it out of the for loop gave me the error I mentioned (about being unable to divide vector by vector etc).

Pedro

[Andrew Gelman]

To do pointwise multiplicaiton or division of vectors you do .* or ./

[Pedro Carvalho]

Trying to do that returned me no matches for real .* vector and suggested vector .* vector, row vector .* row vector, or matrix .* matrix at rate .* dsi.

[Daniel Lee]

We follow math conventions. Real * vector is just that: a * b. If you need vector multiplied by vector, that must be elementwise multiplication. 

The functions, signatures, return types, and doc are all in the manual. If you're using rstan, you can also use the lookup() function. 

Daniel

[Pedro Carvalho]

Right, but when I try doing real * vector it gives me the vector / vector error (which cannot really be corrected by pointwise division because the vectors' sizes are different), and there's also the problem where scale[country] doesn't work because country is a vector (even though what I really mean is scale[country[i]] for every i, as defined in the for loop).

[Bob Carpenter]

> On Dec 17, 2016, at 9:21 AM, Pedro Carvalho <pedro...@gmail.com> wrote:
>
> Right, but when I try doing real * vector it gives me the vector / vector error

If you have (real * vector), the result is vector. Same with (vector * real).
You can do (row_vector * vector) and get a real, or (vector * row_vector)
and get a matrix. But you can't do (vector * vector) because the dimensions
don't conform to matrix multiplication.

> (which cannot really be corrected by pointwise division because the vectors' sizes are different),

Do you have a simple example of the types you want
to define and operations you want to perform? We can help
with specifics of error messages.

> and there's also the problem where scale[country] doesn't work because country is a vector (even though what I really mean is scale[country[i]] for every i, as defined in the for loop).

country is going to be an array if it's an index. and if
it's an array of K elements, then scale[country] will also
have K elements, (scale[country[1]], ..., scale[country[K]])

- Bob

[Pedro Carvalho]

The line in question is:

for(i in 1:N_people_day){
    mu[i] = scale[country[i]] / (1 + exp(-(rate * dsi[i] + x_location))) + y_location[country[i]];
}

mu, country, and dsi are vectors of size N_people_day, scale and y_location are vectors of size N_country (which is the upper bound on the vector country). When I try changing that to:

mu = scale[country] / (1 + exp(-(rate * dsi + x_location))) + y_location[country];

It gives me an error saying that I "Container index must be integer" and also that I can't divide vector by vector... but it now occurs to me I should perhaps try ./ rather than just /. I'll try that.

[Bob Carpenter]

> On Dec 20, 2016, at 1:29 PM, Pedro Carvalho <pedro...@gmail.com> wrote:
>
> The line in question is:
>
> for(i in 1:N_people_day){
> mu[i] = scale[country[i]] / (1 + exp(-(rate * dsi[i] + x_location))) + y_location[country[i]];
> }
>
> mu, country, and dsi are vectors of size N_people_day,

You can't index with a vector---it needs to be an array of integers.
Unlike R, but like C++, Stan distinguishes integer and real values.

- Bob

[Pedro Carvalho]

Yeah, using the pointwise division there worked, and I got this to do what I wanted, thanks!