Specifying a Multinomial Logistic Regression with Unpooled Estimates

Modeling
, ,
1

Apologies if there are simple mistakes in this code, this is probably my second Stan model ever and I still have a lot to learn i.t.o data types/structures and other Stan model specifications.

I am trying to fit a Multinomial Logistic Regression model, where I have a global intercept for each of the score equations and then an individual specific offset term for each of the score equations. I then use the softmax function to pass the result to the multinomial distribution function.

Attempt1:

data {
    int<lower = 0> N; // number of individuals
    int<lower = 2> K; // number of possible scores
    int UID_length;
    int UID[N];
    matrix[N, K] counts;
}
parameters {
    vector[K-1] a; // intercepts
    matrix[K-1, UID_length] b; // association of income with choice
    vector[K-1] a_global;
    
}
model {
    matrix[N, K] p;
    matrix[N, K] s;
    a_global ~ normal( 0 , 1 );

    b[,1] ~ normal( 0 , 0.5 );
    b[,2] ~ normal( 0 , 0.5 );
    b[,3] ~ normal( 0 , 0.5 );
    b[,4] ~ normal( 0 , 0.5 );
    for (i in 1:N){
        s[i, 1] = a_global[1] + b[UID[i], 1];
        s[i, 2] = a_global[2] + b[UID[i], 2];
        s[i, 3] = a_global[3] + b[UID[i], 3];
        s[i, 4] = a_global[4] + b[UID[i], 4];
        s[i, 5] = 0; // pivot
        p[i,] = to_row_vector(softmax(s[i,]'));
        counts[i,] ~ multinomial(p[i, ]);
        }
}
}

An example of the data I have to fit to this model, looks like:

{'counts': array([[ 0,  1,  0,  2,  0],
        [ 0,  0,  1,  1,  0],
        [ 0,  0,  0,  1,  0],
        [ 0,  0,  1,  0,  0],
        [ 0,  6,  5,  7,  0],
        [ 0,  2,  1,  3,  0],
        [ 0,  2,  5,  5,  0],
        [ 0,  2,  4,  3,  0],
        [ 0,  1,  1,  1,  0],
        [ 1,  0,  0,  0,  0]]),
 'UID': array([ 0,  1,  2,  3,  3,  5,  6,  7,  8,  8, 10]),
 'UID_length': 10}

Apart from being unsure about how you specify the indices, I get an error for the last multinomial statement when I run stan.build function:

Building: Semantic error:
   -------------------------------------------------
    30:          s[i, 5] = 0; // pivot
    31:          p[i,] = to_row_vector(softmax(s[i,]'));
    32:          counts[i] ~ multinomial(p[i, ]);
                 ^
    33:          }
    34:  }
   -------------------------------------------------

Ill-typed arguments to '~' statement. No distribution 'multinomial' was found with the correct signature.

Any help or guidance i.t.o. this problem will be very much appreciated, or any reference to relevant Stan documentation.

2 
data {
    int<lower = 0> N; // number of individuals
    int<lower = 2> K; // number of possible scores
    int UID_length;
    int UID[N];
    int counts[N, K];
}
parameters {
    vector[K-1] a; // intercepts
    matrix[K-1, UID_length] b; // association of income with choice
    vector[K-1] a_global;
    
}
model {
    matrix[N, K] p;
    matrix[N, K] s;
    a_global ~ normal( 0 , 1 );

    b[,1] ~ normal( 0 , 0.5 );
    b[,2] ~ normal( 0 , 0.5 );
    b[,3] ~ normal( 0 , 0.5 );
    b[,4] ~ normal( 0 , 0.5 );
    for (i in 1:N){
        s[i, 1] = a_global[1] + b[UID[i], 1];
        s[i, 2] = a_global[2] + b[UID[i], 2];
        s[i, 3] = a_global[3] + b[UID[i], 3];
        s[i, 4] = a_global[4] + b[UID[i], 4];
        s[i, 5] = 0; // pivot
        p[i,] = to_row_vector(softmax(s[i,]'));
        counts[i,] ~ multinomial(p[i, ]');
        }
}

multinomial is defined for array[] int y and vector only.

3

Thanks @andre.pfeuffer that is the solution, I’m glad it was such a small issue and change that I had to make to get this working. I probably misunderstood the documentation, it is sometimes a bit confusing to keep track of data types, sizes of data structure, and when something can be vectorized and when not.

4

An updated, cleaner working version (incorporating your change) for future reference:

Attempt 2:

data {
    int<lower = 0> N;
    int<lower = 2> K;
    int counts[N, K];
    int UID_length;
    int UID[UID_length];
}
parameters{
    vector[UID_length] b; 
    vector[K-1] a_global;
    
}

model{
    vector[K] p;
    vector[K] s;
    a_global ~ normal( 0 , 1 );
    b ~ normal( 0 , 0.5 );

    for (i in 1:N){
        s[1] = a_global[1] + b[UID[i]];
        s[2] = a_global[2] + b[UID[i]];
        s[3] = a_global[3] + b[UID[i]];
        s[4] = a_global[4] + b[UID[i]];
        s[5] = 0; // pivot
        p = softmax(s);
        counts[i,] ~ multinomial(p);
        }
}