Multi-level ordered logistic regression: where to put the varying part of the model?

Modeling
,
1

Hello there,

Firstly, congrats on the new forum. Looks very streamlined and stylish.
Now my query:

I’m fitting a multi-level ordered logistic regression model, which seems to converge well. But when I started checking on the old forum, and noticed that the varying intercept part is not on the lineal model part, but rather on the cutpoints part (hope that’s understandable). The model I fitted though has the varying intercepts on the lineal model part, something like this:

where the last bit, b1i, is the varying incercept. This kind parameterization (and picture) I took from a frequentist book though (Applied Longitudinal Analysis, by Fitzmaurice). Is this an OK way to add varying intercepts to a ordered logistic regression model?
Please see below the full code for my model:

data{
int<lower=1> N;
int<lower=1> N_sample;
int<lower=1> N_cage;
int histo[N];
real AGD[N];
real Des[N];
real cyst[N];
real PRV[N];
real pox[N];
real Ten[N];
int sample[N];
int cage[N];
real Temp[N];
real time_w[N];
matrix[N_sample,N_sample] Dmat;
}
parameters{
ordered[3] cutpoints;
real bt;
real bwater_t;
real btemp_timew;
real bagd;
real bdes;
real bcyst;
real bprv;
real bpox;
real bten;
real bagd_des;
real bagd_cyst;
real bagd_prv;
real bagd_pox;
real bagd_ten;
vector[N_sample] a_sample;
vector[N_cage] a_cage_raw;
real<lower=0> sigma_cage;
real<lower=0> etasq;
real<lower=0> rhosq;
real<lower=0> sigma;
}
transformed parameters{
vector[N_cage] a_cage;
a_cage = a_cage_raw*sigma_cage;
}

model{
matrix[N_sample,N_sample] SIGMA_Dmat;
vector[N] phi;
sigma ~ cauchy( 0 , 1 );
sigma_cage ~ cauchy( 0 , 1 );
rhosq ~ cauchy( 0 , 1 );
etasq ~ cauchy( 0 , 1 );
for ( i in 1:(N_sample-1) )
for ( j in (i+1):N_sample ) {
SIGMA_Dmat[i,j] = etasqexp(-rhosqpow(Dmat[i,j],2));
SIGMA_Dmat[j,i] = SIGMA_Dmat[i,j];
}
for ( k in 1:N_sample )
SIGMA_Dmat[k,k] = etasq + sigma;
a_sample ~ multi_normal( rep_vector(0,N_sample) , SIGMA_Dmat );
a_cage_raw ~ normal(0, 1);
cutpoints ~ normal( 0 , 5 );
bagd_ten ~ normal( 0 , 1 );
bagd_pox ~ normal( 0 , 1 );
bagd_prv ~ normal( 0 , 1 );
bagd_cyst ~ normal( 0 , 1 );
bagd_des ~ normal( 0 , 1 );
bten ~ normal( 0 , 1 );
bpox ~ normal( 0 , 1 );
bprv ~ normal( 0 , 1 );
bcyst ~ normal( 0 , 1 );
bdes ~ normal( 0 , 1 );
bagd ~ normal( 0 , 1 );
btemp_timew ~ normal( 0 , 1 );
bwater_t ~ normal( 0 , 1 );
bt ~ normal( 0 , 1 );
for ( i in 1:N ) {
phi[i] = bt * Temp[i] + bwater_t * time_w[i] + btemp_timew * Temp[i] * time_w[i] + bagd * AGD[i]

  • bdes * Des[i] + bcyst * cyst[i] + bprv * PRV[i] + bpox * pox[i] + bten * Ten[i]
  • bagd_des * AGD[i] * Des[i] + bagd_cyst * AGD[i] * cyst[i]
  • bagd_prv * AGD[i] * PRV[i] + bagd_pox * AGD[i] * pox[i] + bagd_ten * AGD[i] * Ten[i]
  • a_sample[sample[i]] + a_cage[cage[i]];
    }
    for ( i in 1:N )
    histo[i] ~ ordered_logistic( phi[i] , cutpoints );
    }
    generated quantities{
    matrix[N_sample,N_sample] SIGMA_Dmat;
    vector[N] phi;
    real log_lik[N];

for ( i in 1:(N_sample-1) )
for ( j in (i+1):N_sample ) {
SIGMA_Dmat[i,j] = etasqexp(-rhosqpow(Dmat[i,j],2));
SIGMA_Dmat[j,i] = SIGMA_Dmat[i,j];
}
for ( k in 1:N_sample )
SIGMA_Dmat[k,k] = etasq + sigma;

for ( i in 1:N ) {
phi[i] = bt * Temp[i] + bwater_t * time_w[i] + btemp_timew * Temp[i] * time_w[i] + bagd * AGD[i]

  • bdes * Des[i] + bcyst * cyst[i] + bprv * PRV[i] + bpox * pox[i] + bten * Ten[i]
  • bagd_des * AGD[i] * Des[i] + bagd_cyst * AGD[i] * cyst[i]
  • bagd_prv * AGD[i] * PRV[i] + bagd_pox * AGD[i] * pox[i] + bagd_ten * AGD[i] * Ten[i]
  • a_sample[sample[i]] + a_cage[cage[i]];
    }

for ( i in 1:N )
log_lik[i] = ordered_logistic_lpmf( histo[i] | phi[i] , cutpoints );
}
"

Thanks,

Tada

2

Not our design—just took it from Discourse. But I keep losing track of posts like this one, so for me, it’s been a big step backwards. I do like the ability to drop in markdown—it’d be even better if it gave us MathJax.

The main thing you have to worry about is identifiability. If you have too many intercepts floating around, you get additive non-identifiabilities where adding a number to all the cut points and to the intercept leaves the likelihood unchanged.

You can vectorize that calculation of phi now and it should be more efficient. It just requires all the continuous containers to be vectors.

Is this right?

for ( k in 1:N_sample )
SIGMA_Dmat[k,k] = etasq + sigma;
a_sample ~ multi_normal( rep_vector(0,N_sample) , SIGMA_Dmat );

It looks like you set a different element each run through. If that’s the case, there’s not much you can do efficiency-wise (unless you could figure out how this would factor down to Cholesky factors).

1 Like
3

Thanks, Bob. I’m glad to know it’s a valid parameterization for the varying intercepts part.
I think I had too many them in the model though, which gave poor model convergence. After taking the highest level out, the model ran well. Even with the Gaussian process not being so efficient.

Thanks for getting back to me on this.

Cheers,

Tada