Hierarchical multinomial model

Modeling
1

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 {
//priors
Theta_raw ~ normal(0, 10);
L_Omega ~ lkj_corr_cholesky(4);

//likelihood
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]);
}
}
}"
2

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.

3

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")

X<-model.matrix(~X1+X2+X3+X4+X5+X6+X7+X8+
                  X9,mydata)
				  
#number of regression params
K<-dim(X)[2] 
#the regression slopes
Beta<-c(0,0,0,0,0,0,0,0,0,0)
#pars
Theta<-0.1
L_Sigma<-0.1
L_Omega<-0.1

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
eth<-c(0,1)
rel<-c(0,1)
mix<-c(0,1)

#fit the model
m_mult<-stan(file = "multinom_hier.stan",data = list(R=R,C=C,S=S,K=K,X=X,rel=rel,
                                                     eth=eth,mix=mix),
             pars=c("L_Sigma","Theta","L_Omega","Beta"))
4

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?

5

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

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

int<lower=0,upper=1>?

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

6

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

7

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

8

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

9

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