Hi, I am following Jim Savage rank ordered logit approach generally. However, instead of his beta_individual definition, I want to generate beta_individual (a matrix) from a multinomial distribution. Then I get the error because stan can only generate a vector or row_vector rather than a matrix under the distribution. Any suggestions? (The sigma equation is not correct while I just list an example to show sigma is a ixp matrix.)
Stan code for the ranked choice model
// saved as ranked_rcl.stan
functions {
real rank_logit_lpmf(int[] rank_order, vector delta) {
// We reorder the raw utilities so that the first rank is first, second rank second...
vector[rows(delta)] tmp = delta[rank_order];
real out;
// ... and sequentially take the log of the first element of the softmax applied to the remaining
// unranked elements.
for(i in 1:(rows(tmp) - 1)) {
if(i == 1) {
out = tmp[1] - log_sum_exp(tmp);
} else {
out += tmp[i] - log_sum_exp(tmp[i:]);
}
}
// And return the log likelihood of observing that ranking
return(out);
}
}
data {
int N; // number of rows
int T; // number of inidvidual-choice sets/task combinations
int I; // number of Individuals
int P; // number of covariates that vary by choice
int P2; // number of covariates that vary by individual
int K; // number of choices
int rank_order[N]; // The vector describing the index (within each task) of the first, second, third, ... choices.
// In R, this is order(-utility) within each task
matrix[N, P] X; // choice attributes
matrix[I, P2] X2; // individual attributes
int task[T]; // index for tasks
int task_individual[T]; // index for individual
int start[T]; // the starting observation for each task
int end[T]; // the ending observation for each task
}
parameters {
vector[P] beta; // hypermeans of the part-worths
matrix[P, P2] Gamma; // coefficient matrix on individual attributes
vector<lower = 0>[P] tau; // diagonal of the part-worth covariance matrix
matrix[I, P] z; // individual random effects (unscaled)
matrix[I, P] beta_individual ;
}
transformed parameters {
// here we use the reparameterization discussed on slide 30
matrix[I, P] sigma=X2*Gamma';
}
model {
// priors on the parameters
tau ~ normal(0, .5);
beta ~ normal(0, 1);
to_vector(z) ~ normal(0, 1);
to_vector(Gamma) ~ normal(0, 1);
**beta_individual ~ multi_normal(beta,sigma);**
// log probabilities of each choice in the dataset
for(t in 1:T) {
vector[K] utilities; // tmp vector holding the utilities for the task/individual combination
// add utility from product attributes with individual part-worths/marginal utilities
utilities = X[start[t]:end[t]]*beta_individual[task_individual[t]]';
rank_order[start[t]:end[t]] ~ rank_logit(utilities);
}
}