Ordered logistic regression with latent responses

I am trying to run an ordered logistic regression model, using quite special data.

In my data, I am measuring one outcome per respondent, which can be coded as 1, 2 and 3 (i.e., the survey responses “no”, “I am not sure” and “yes”).
Then, I have two different input variables, with which the respondents were treated. However, quite a big subset of my participants did not receive only one treatment with the two input variables but several treatments at the same time (i.e., 974 respondents had one unique treatment, 637 respondents had two treatments, 268 respondents had three treatments,…).

As my aim is, to estimate the effect of my two input variables on the probability of giving a particular response, I thought about the following model specification:

Actually, each of the treatments we observe (i.e., a combination of the two input variables) would lead to a specific response, i.e., either yes, I am not sure or no. Or to put it differently, each of the treatments is associated with a particular probability of being in one of the three categories.
However, we only observe one outcome per panelist and not one outcome per treatment (which we would be actually interested in).
So, if we have a respondent who received two treatments but gave only one response, the response we eventually measure is actually already a combination of the probabilities of being in one of the categories per treatment.
If the respondent received two treatments, I would assume that the probability of being in the yes category is increased, wheras the probability of being in the no category is decreased (this is actually an assumption that can, of course, be discussed - however, I should be able to support this assumption with theory).

So, I am estimating a probability (theta) of being in one of the categories for each treatment. In a next step, I combine those distinct probabilities to an overall probability (THETA) per respondent, which is a combination of the distinct treatment probabilities per respondent.
For example, a respondent has probabilities theta of being in the “no”, “i am not sure” and “yes” category of (0.6, 0.2, 0.2) for treatment 1 and (0.1, 0.5, 0.4) for treatment 2.

The combined probabilities THETA would then be calculated as follows:
no-category: 0.6 * 0.1 = 0.06
yes-category: 1 - ((1 - 0.2) * (1 - 0.4)) = 0.52
not sure-category: 1- 0.06 - 0.52 = 0.42

This is the code I am using at the moment:

data {
  int<lower=2> K;
  int<lower=0> N;
  int<lower=1> L;
  int<lower=1,upper=K> CrR[L];
  vector[N] var1;
  vector[N] var2;
  int<lower=1,upper=L> ll[N];
  int<lower=1,upper=N> indices1[L];
  int<lower=1,upper=N> indices2[L];
parameters {
  real b1;
  real b2;
  ordered[K-1] c;
model {
  vector[N] theta[K];
  vector[K] THETA;
  for (n in 1:N) {
    theta[1,n] = 1 - (inv_logit(var1[n]*b1 +var2[n]*b2 - c[1]));
    for(k in 2:K-1)
      theta[k,n] = (inv_logit(var1[n]*b1 + var2[n]*b2 - c[k-1])) - (inv_logit(var1[n]*b1 + var2[n]*b2 - c[k]));
    theta[K,n] = (inv_logit(var1[n]*b1 + var2[n]*b2 - c[K-1]));
    for (l in 1:L) {
      THETA[1] = theta[1,indices1[l]];
      THETA[3] = theta[3,indices1[l]];
      for (m in 1:(indices2[l]-indices1[l])) {
        if (m > 0) {
        THETA[1] = THETA[1] * theta[1,(indices1[l]+m)];
        THETA[3] = 1 - ((1-THETA[3]) * (1-theta[3,(indices1[l]+m)]));
      THETA[2] = 1 - THETA[1] - THETA[3];
      CrR[l] ~ categorical(THETA); 

Please note that I only implemented the model code for three categories at the moment. I would still need to generalize the code to a flexible number of categories.

K… number of categories
N… overall number of observations (one observation for each treatment)
L… overall number of respondents
CrR… measured responses
var1… input variable 1
var2… input variable 2
ll… respondent ID
indices1… indices variable for each respondent, stating which treatment is the first treatment belonging to respondent ll
indices2… indices variable for each respondent, stating which treatment is the last treatment belonging to respondent ll

So, I first tried running the model with those respondents only, who received only one treatment, such that I can directly compare it to the normal ordered logistic model.
And actually, I get the same estimates as in the normal ordered logistic model.

So, what is your overall oppinion about my modelling idea? Is my way of calculating the overall probabilities per respondent proabably not a valid approach, i.e., influencing final results to much?

I am happy about any comments and feedback you can give me!
Thank you in advance,

maybe I did not understand your model description well, but I think what you describe is a relatively non-standard way to code an ordinal model. You might want to check out Paul Bürkner’s Ordinal Regression Models in Psychology: A Tutorial which describes how these models are handled in the brms package (which uses Stan as backend). The brms package (for R only) could also be sufficient for you to code the model via a generalized R formula syntax which should be easier than debugging your own Stan model. If brms would not be flexible enough for your case, you can also let it to produce Stan code for a similar and then tweak it to suit your needs.

Hope that helps.