Problem with multi-logit regression parameterization

Modeling
1

Hi, I’m trying to fit a particular parameterization of a multi-logit model. In particular, I have some survey data and am trying to recreate an analysis similar to what was done here. I can get the model to run (and fit perfectly well) if I use the reference manual’s ordered logistic function (page 138) or ordered probit (page 139), but as best I can tell neither of those give me the mu parameter that I really want.

So I’ve tried to recreate the model using normal_cdf, which parameterizes in terms of mu and sigma, but I start running into all sorts of output errors - depending on what I do exactly, I get rejected initial values, invalid ordered vector values, and/or not a valid simplex errors. I’ve put in a simple version of the code below; can someone give me a hand? If you want to run the model, the simple data set I’m starting with is one survey question with six responses: n_cat = 5, n_resp = 6, y = 3,3,3,3,2,5 .

Thanks!

data {
  int<lower=2> n_cat;
  int<lower=1> n_resp;
  int<lower=1,upper=n_cat> y[n_resp];
}
parameters {
  real beta;
  real<lower=0> sigma;
  ordered[n_cat-1] thresh;
}
model {
  vector[n_cat] theta;
  vector[2] junk;
  junk[1]=0.0;
  beta ~ normal((n_cat+0.5)/2,n_cat);
  sigma ~ cauchy(0,20);
  for (k in 2:(n_cat-2))
    thresh[k] ~ normal(k+0.5,2);
  for (n in 1:n_resp) {
    theta[1] = 1-normal_cdf(thresh[1],beta,sigma);
    for (k in 2:(n_cat-1)) {
      junk[2] = normal_cdf(thresh[k-1],beta,sigma) - normal_cdf(thresh[k],beta,sigma);
      theta[k] = max(junk);
    }
    theta[n_cat] = normal_cdf(thresh[n_cat-1],beta,sigma);
    y[n] ~ categorical(theta);
  }
}

[edit: escaped code and auto-indented]

2

I’m not sure what that models’ trying to do with max(junk).

Are you sure that the theta[n] form simplexes? If so, then the last entry can just subtract the previous entries.

The ordering can fail if you wind up collapsing two values wind up collapsing on each other or drifting off to infinity.

3

In Kruschke’s code, the middle (2 to n_cat-1) thetas have a ‘max’ term to ensure that they’re at least 0. It was written for JAGS, and I don’t know enough to say what exactly it’s for.

After some additional reading and testing I got a model that works, including on a bigger data set. The key appears to be specifying some initial values. I used

initf ← function() list(beta=mean(y), sigma=2, thresh= 1:(n_cat-1))

Blockquote
data {
int<lower=2> n_cat;
int<lower=1> n_resp;
int<lower=1,upper=n_cat> y[n_resp];
}
parameters {
real beta;
real<lower=0> sigma;
ordered[n_cat-1] thresh;
}
model {
vector[n_cat] theta;
beta ~ normal((n_cat+1)/2,2);
sigma ~ cauchy(0,20);
for (k in 1:(n_cat-1))
thresh[k] ~ normal(k+.5,2);
for (n in 1:n_resp) {
theta[1] = normal_cdf(thresh[1],beta,sigma);
theta[n_cat] = 1-normal_cdf(thresh[n_cat-1],beta,sigma);
for (k in 2:(n_cat-1))
theta[k] = normal_cdf(thresh[k],beta,sigma) - normal_cdf(thresh[k-1],beta,sigma);
y[n] ~ categorical(theta);
}
}

4

Maybe to avoid with rounding issues accidentally making something negative.

Great. That’s usually because of numerical issues in the regression terms. Rescaling everything so you expect the posteriors for all the parameters to be roughly unit scale (normal(0, 1) or something like uniform(-2, 2) like we use for init).

1 Like