scores[N]; // list of scores
}
parameters {
ordered[3] cutpoints; // 3 cutpoints for 4 scores, keeps in order
}
model {
cutpoints ~ normal(0,10); // our prior for cutpoints (actually like -3-\>p=.05,-2,-1,0-\>p=.50 etc)
for (i in 1:N) {
scores[i] ~ ordered\_logistic(0, cutpoints);
}
}
generated quantities {
vector[3] cumprob; // cumulative probabilities
vector[4] prob; // probabilities
cumprob = inv\_logit(cutpoints);
prob[1] = cumprob[1];
prob[2] = cumprob[2] - cumprob[1];
prob[3] = cumprob[3] - cumprob[2];
prob[4] = 1 - cumprob[3];
}
"
fit \<- stan(model\_code = stanstratrubric\_ord, data = rs, iter = 4000, warmup = 2000,chains = 4,control = list(adapt\_delta = 0.99),refresh = 0)
get\_posterior\_mean(fit,par=c(“prob[1]”,“prob[2]”,“prob[3]”,“prob[4]”))
---
### Author:  [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)
#### Post date: [August 28, 2019, 3:29pm UTC](https://discourse.mc-stan.org/t/ordered-logistic-small-sample-bias/10578/2 "2019-08-28T15:29:28Z")
This is to be expected. If you output the full summary via `print(fit)` the output for the small N case would be something like:
```
Inference for Stan model: 2c592c48b51f3916e2aa109e2a468f2c.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
cutpoints[1] -3.05 0.03 1.31 -6.23 -3.81 -2.86 -2.14 -0.99 2353 1
cutpoints[2] -1.01 0.01 0.72 -2.52 -1.47 -0.98 -0.52 0.32 5584 1
cutpoints[3] 0.53 0.01 0.68 -0.79 0.08 0.51 0.96 1.92 6939 1
cumprob[1] 0.08 0.00 0.07 0.00 0.02 0.05 0.11 0.27 3294 1
cumprob[2] 0.29 0.00 0.13 0.07 0.19 0.27 0.37 0.58 6033 1
cumprob[3] 0.62 0.00 0.14 0.31 0.52 0.63 0.72 0.87 6892 1
prob[1] 0.08 0.00 0.07 0.00 0.02 0.05 0.11 0.27 3294 1
prob[2] 0.21 0.00 0.12 0.04 0.13 0.19 0.28 0.49 5542 1
prob[3] 0.33 0.00 0.13 0.10 0.23 0.32 0.42 0.62 3665 1
prob[4] 0.38 0.00 0.14 0.13 0.28 0.37 0.48 0.69 6892 1
lp__ -13.56 0.03 1.25 -16.73 -14.14 -13.22 -12.66 -12.13 2398 1
Samples were drawn using NUTS(diag_e) at Wed Aug 28 17:23:00 2019.
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).
```
While the means for `prob[i]` jump around the actual values, you may notice the large uncertainty due to the small sample size. I.e. the results tell you that you shouldn’t trust the mean estimates very much, as e.g. for `prob[1]` the full range from `0` to `0.27` lies withing the 95% credible interval. With larger N the uncertainty decreases and means become more informative. I think the take home message is that you cannot only report the mean estimate and reporting the associated uncertainty (95% CI or at least the sd) is very important.
Finally, two minor notes:
- Unless you are aiming for something very special, you might want to use the `brms` package which will give you a well tested implementation of ordinal regression almost for free.
- You can use triple backticks (```) to make code blocks in the post neatly formatted
Hope that helps.
---
### Author:  [@Jim\_Hatton](https://discourse.mc-stan.org/u/Jim_Hatton)
#### Post date: [August 28, 2019, 4:12pm UTC](https://discourse.mc-stan.org/t/ordered-logistic-small-sample-bias/10578/3 "2019-08-28T16:12:48Z")
Thanks Martin for your reply.
As I understand it, the “biased” means are caused by the fact that the underlying searcher is bumping against the end bounds (0 and 1), though it still seems that means should center about the sample values.
I will look into using the backticks.
As far as using brms, I prefer a process of knowledge accretion which in this case, as you can see, I try to build increasingly more complex models.
I really appreciate you taking the time to reply.
---
### Author:  [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)
#### Post date: [August 29, 2019, 10:43am UTC](https://discourse.mc-stan.org/t/ordered-logistic-small-sample-bias/10578/4 "2019-08-29T10:43:33Z")
> [@Jim\_Hatton](#):
>
> As I understand it, the “biased” means are caused by the fact that the underlying searcher is bumping against the end bounds (0 and 1), though it still seems that means should center about the sample values.
Oh, I see, I misunderstood you earlier. So the bias you care about is that when you fit the model repeatedly, the `prob[i]` are always skewed the same? This is IMHO still unsurprising, but for a different reason: the sampler is expected to give very similar results when run on the same data. It should jump around the true values when the data would be simulated from exactly the model you’ve encoded.
But you have one fixed dataset and this dataset is not very typical of your prior. If I simulate the cutpoints according to your prior, I get (please double check the code):
```no-highlight
N_sims <- 1000
cutpoints <- array(NA_real_, c(N_sims, 3))
probs <- array(NA_real_, c(N_sims, 4))
for(i in 1:N_sims) {
cutpoints[i,] <- sort(rnorm(3, 0, 10)) #Enforce the ordering
cumprob <- 1/(1 + exp(-cutpoints[i,]))
probs[i,] <- diff(c(0, cumprob, 1))
}
# Plot implied univariete prior distributions
hist(probs[,1])
hist(probs[,2])
hist(probs[,3])
hist(probs[,4])
# Prior probability probs[1] < 0.08
mean(probs[,1] < 0.08)
```
The implied prior for `probs[1]` is:
And the prior probability that `probs[1] < 0.08` is around `0.8`. The prior actually strongly favors small probabilities for the extreme categories. Since you have a small dataset, the prior has noticeable influence on your inferences (but the wide uncertainty signals that you cannot have a lot of trust in the precise estimates anyway).
Indeed, when I change the prior, the results change. E.g. when I put `cutpoints ~ normal(0, 2);` the results change to something like:
```no-highlight
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
cutpoints[1] -2.15 0.02 0.86 -3.99 -2.72 -2.09 -1.55 -0.63 3033 1
cutpoints[2] -0.71 0.01 0.61 -1.94 -1.11 -0.70 -0.29 0.45 6094 1
cutpoints[3] 0.63 0.01 0.61 -0.53 0.21 0.61 1.03 1.87 7196 1
cumprob[1] 0.13 0.00 0.09 0.02 0.06 0.11 0.17 0.35 3164 1
cumprob[2] 0.34 0.00 0.13 0.13 0.25 0.33 0.43 0.61 6333 1
cumprob[3] 0.64 0.00 0.13 0.37 0.55 0.65 0.74 0.87 7142 1
prob[1] 0.13 0.00 0.09 0.02 0.06 0.11 0.17 0.35 3164 1
prob[2] 0.21 0.00 0.11 0.05 0.13 0.20 0.28 0.46 5054 1
prob[3] 0.30 0.00 0.12 0.10 0.21 0.29 0.38 0.56 4043 1
prob[4] 0.36 0.00 0.13 0.13 0.26 0.35 0.45 0.63 7142 1
lp__ -14.46 0.03 1.27 -17.55 -15.06 -14.14 -13.53 -13.01 2484 1
```
Showing the influence of prior choice.
It is actually quite hard to figure out a prior for the cutpoints that is not very tight AND gives large prior probability to a wide spread of the ordinal responses, Mike Betancourt has an interesting attempt in his case study on the topic ([Ordinal Regression](https://betanalpha.github.io/assets/case_studies/ordinal_regression.html)), but I have no direct experience with it.
Hope I’ve made things clearer for you.
---
### Author:  [@Jim\_Hatton](https://discourse.mc-stan.org/u/Jim_Hatton)
#### Post date: [August 29, 2019, 4:04pm UTC](https://discourse.mc-stan.org/t/ordered-logistic-small-sample-bias/10578/5 "2019-08-29T16:04:06Z")
Hi Martin, Thanks for the explanation. My waking thought this morning was that if a category has a non-zero count, then I know the probability for that category is non-zero and if I have a good sense of the population size, I have a strict lower bound on the probability. I was wondering if there was a way to incorporate that into the model. As you have shown, small sample results are sensitive to the prior. thx, jim
---
### Author:  [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)
#### Post date: [August 30, 2019, 7:08am UTC](https://discourse.mc-stan.org/t/ordered-logistic-small-sample-bias/10578/6 "2019-08-30T07:08:43Z")
> [@Jim\_Hatton](#):
>
> My waking thought this morning was that if a category has a non-zero count, then I know the probability for that category is non-zero and if I have a good sense of the population size, I have a strict lower bound on the probability.
While I am by no means an expert on Bayesian computation, I am not so sure about this: the ordered logistic model already implies that all categories have non-zero (although possibly negligible) probability, regardless of the data you collect. Seeing at least one response in the category implies the probability is not-totally-negligible, but that’s IMHO about it. Same with “strict lower bound” - yes, the posterior probability that the category has very low true probability would tend to be vanishingly small, but there is nothing “strict” about it. E.g. after collecting 100 responses, 10 of which are in category 1, it is still possible that the true probability of category 1 is \<0.001, although it is highly unlikely.
All in all, I would guess that the reasoning you care about is mostly (if not completely) handled by the computed posterior probability - everything that is implied by the model, will be reflected there. Unless I am missing something, your comments do not seem to introduce any new knowledge beyond some informal reasoning about the model.
Hope that helps!