Hi all,

I am estimating a multinomial regression model with a 4 category outcome variables. The predictors are

```
male: (1=male)
ASBG04: # of books in home 1) 0-10, 2) 11-25, 3) 26-100, 4) 101-200 5) More than 200
ASBGSMR: Students motivated to read (4pt, strongly agree to strongly disagree)
```

The model set up is as follows

```
N <- nrow(Canadareg1)
f <- as.formula("ASBR07A ~ ASBG04 + ASBGSMR + male")
M <- model.matrix(f,Canadareg1)
data.list <- list(N=nrow(M),
K=length(unique(Canadareg1[,1])),
D=ncol(M),
x=M,
ASBR07A=as.numeric(Canadareg1[,1]))
modelString <- "
data {
int <lower=2> K; // This is 4, the number of outcomes categories
int <lower=0> N;
int <lower=1> D; // This is the number of columns in the design matrix: 4
int <lower = 1, upper = K> ASBR07A[N];
matrix[N, D] x; // This will be N by 2 matrix of data
}
parameters {
matrix[D, K] beta; // This is a 4 x 4 matrix of betas
}
transformed parameters {
matrix[N, K] x_beta = x * beta; // N x 4
}
model {
to_vector(beta) ~ normal(0, 5);
for (i in 1:N)
ASBR07A[i] ~ categorical_logit(x_beta[i]');
}
generated quantities {
int<lower=1,upper=K> ASBR07A_rep[N];
for (i in 1:N){
ASBR07A_rep[i] = categorical_logit_rng(x_beta[i]');
}
}
"
```

Everything runs find, and the output looks is as follows.

```
## Inference for Stan model: 87a49ecccc72c4a5df74a69624386dba.
## 1 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=500.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## beta[1,1] -2.33 0.25 2.26 -6.21 -3.99 -2.42 -0.80 2.06 83 1.00
## beta[1,2] 1.00 0.24 2.23 -3.01 -0.61 0.84 2.50 5.34 86 1.00
## beta[1,3] 0.90 0.25 2.26 -3.20 -0.79 0.85 2.44 5.04 84 1.00
## beta[1,4] -0.17 0.24 2.26 -4.18 -1.80 -0.17 1.28 4.33 87 1.00
## beta[2,1] 0.75 0.44 2.69 -4.85 -1.01 1.08 2.68 5.17 37 1.00
## beta[2,2] 0.50 0.44 2.68 -4.98 -1.26 0.87 2.45 5.00 37 1.00
## beta[2,3] 0.49 0.44 2.69 -5.14 -1.28 0.85 2.47 4.98 37 1.00
## beta[2,4] 0.17 0.45 2.69 -5.40 -1.61 0.57 2.15 4.65 36 1.00
## beta[3,1] 0.11 0.44 2.53 -5.02 -1.68 0.09 1.85 5.44 33 1.02
## beta[3,2] -0.17 0.44 2.53 -5.29 -1.96 -0.20 1.57 5.19 33 1.02
## beta[3,3] -0.19 0.44 2.53 -5.28 -1.98 -0.22 1.56 5.13 33 1.02
## beta[3,4] 0.03 0.44 2.53 -5.13 -1.74 0.02 1.76 5.38 33 1.02
## beta[4,1] 1.19 0.39 2.54 -3.81 -0.47 1.17 3.23 5.79 42 1.01
## beta[4,2] 0.65 0.39 2.53 -4.20 -1.12 0.60 2.67 5.22 42 1.01
## beta[4,3] 0.35 0.39 2.53 -4.49 -1.34 0.30 2.38 4.90 42 1.01
## beta[4,4] 0.41 0.39 2.54 -4.48 -1.29 0.33 2.40 5.07 42 1.00
##
## Samples were drawn using NUTS(diag_e) at Mon Nov 16 21:06:41 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
```

As you can see, I am getting the results for ASBGSMR for each of the four categories, but not the others. I tried a different ordering of the predictors but receive the same results. Iâ€™m not sure what is going on here.

Thanks