Hierarchical multinomial model

Hi all,

I’d appreciate help on the code below.

I’m trying to fit a hierarchical multinomial logit model.

I tried to adapt to my needs this example http://discourse.mc-stan.org/t/speeding-up-a-hierarchical-multinomial-logit-model/1538/4 but something went wrong.

I have 800+ observations of respondents, 3 response choices, 9 independent variables. The data is nested within 7 regions.

Something must be wrong with the syntax I provide. The error I get is:

[Error in new_CppObject_xp(fields$.module, fields$.pointer, ...) : 
  no valid constructor available for the argument list
failed to create the sampler; sampling not done]

If you see obvious errors in the code below, please point out.

 model<-"data {
  int<lower=3> C; // Number of alternatives (choices) in each scenario
int<lower=1> K; // Number of respondent level covariates
int<lower=1> R; // Number of respondents
int<lower=1> S; // Number of regions
int<lower=1,upper=C> rel[R, S]; // religion
int<lower=1,upper=C> eth[R, S]; // ethnicity
int<lower=1,upper=C> mix[R, S]; // mixed
matrix[C, K] X[R, S]; // matrix of attributes for each obs

parameters {
vector[K] Beta[R];
vector[K - 1] Theta_raw;
cholesky_factor_corr[K] L_Omega;
vector<lower=0, upper=pi()/2>[K] L_sigma_unif;

transformed parameters {
vector<lower=0>[K] L_sigma;
matrix[K, K] L_Sigma;
vector[C] XB[R, S];
vector[K] Theta;

for (k in 1:K) {
L_sigma[k] = 2.5 * tan(L_sigma_unif[k]);

L_Sigma = diag_pre_multiply(L_sigma, L_Omega);

Theta[1] = 0;
for (k in 1:(K-1)) {
Theta[k + 1] = Theta_raw[k];

for (r in 1:R) {
for (s in 1:S) {
XB[r,s] = X[r,s] * Beta[r];

model {
Theta_raw ~ normal(0, 10);
L_Omega ~ lkj_corr_cholesky(4);

Beta ~ multi_normal_cholesky(Theta, L_Sigma);
for (r in 1:R) {
for (s in 1:S) {
rel[r,s] ~ categorical_logit(XB[r,s]);
eth[r,s] ~ categorical_logit(-XB[r,s]);
mix[r,s] ~ categorical_logit(XB[r,s]);

What’s your data look like?

This is similar to a question from a few days ago: Help! trying deprecated constructor; please alert package maintainer Error in new_CppObject_xp(fields$.module, fields$.pointer, ...) . In that case the data being passed into the model didn’t match the dimensions declared in the data block.

The model itself compiles fine for me.

That’s good, thanks:)

Then, something must wrong in the run syntax, I guess. Did it specify the Ys correctly?

write(model, file = "multinom_hier.stan")

#number of regression params
#the regression slopes

R<-867 #sample size
C<-3 #number of choices
S<-21 #number of alternatives of each respondent (7 regions X 3 response choices)

#the response, Y1,Y2,Y3

#fit the model
m_mult<-stan(file = "multinom_hier.stan",data = list(R=R,C=C,S=S,K=K,X=X,rel=rel,

rel, eth, and mix all have size 2. They’re declared in the model as

int<lower=1,upper=C> rel[R, S]; // religion
int<lower=1,upper=C> eth[R, S]; // ethnicity
int<lower=1,upper=C> mix[R, S]; // mixed

Which would be 867x21 arrays of integers, so that’s one problem.

Could you check dim(X) as well?

I see, thanks. Still throws the same error though.

If the Ys are dichotomous (0,1) indicators, should they be declared as


dim(X) is 10 that’s what it’s supposed to be, right?

Hmm, that definitely doesn’t sound right.

Have you had a look at the Github repo for that model you linked? There is some code in there which might help you set up your data: https://github.com/ksvanhorn/ART-Forum-2017-Stan-Tutorial/blob/master/3_hmnl/hmnl.R

At least you could look and see how that ksvanhorn did it and copy-paste to get started. There’s a section on generating synthetic data to test the model in the presentation: https://cdn.rawgit.com/ksvanhorn/ART-Forum-2017-Stan-Tutorial/9bffc2b9/3_hmnl/hmnl.html#1

Thanks for the heads up on those links. Hopefully, they will help. I should certainly try the model on simulated data.

You might take a look at the code from this paper: https://link.springer.com/article/10.1007/s00265-017-2363-8

This doesn’t look like it’ll produce a lower-triangular matrix.

This can be vectorized as

for (r in 1:R) {
  rel[r] ~ categorical_logit(XB[r, ])

I’d think you could maybe get one more layer of vectorization out of this one, but I’m not sure as the result is already a vector (it’s much faster when things can be done as matrix multiplication rather than repeated vector multiplication).