Hierarchical probit model for longitudinal data

I have an issue to code a hierarchical probit model for longitudinal data in stan.

To define the model, let the index i denote the individual, and the index j denote time.
Let Y_{i,j} be the binary reponse of interest of patient i at time j.
Let R_{i,j} be the continuous variable (data) on which depend Y_{i,j}, the higher R_{i,j} is, the higher chance of having Y_{i,j}=1.

The model to code is
Y_{i,j}=\left\{ \begin{array}{l} 0 \quad \text{if } Z_i>R_{i,j}\\ 1 \quad \text{if } Z_i \leq R_{i,j}\\ \end{array} \right.

with Z_i a latent variable, constant for a same individual, with Z_i \sim \mathcal{N}(\mu,\tau^2)

We can for example assume as priors \mu \sim \mathcal{N}(0,5²) and \tau \sim \text{half-Cauchy}(0,5)

I don’t know how to add the constraint that Z_i is constant for a same individual.

For example, I can see how to code the model
Y_{i}=\left\{ \begin{array}{l} 0 \quad \text{if } Z_{i}>R_{i}\\ 1 \quad \text{if } Z_{i} \leq R_{i}\\ \end{array} \right.

with Z_{i} a latent variable Z_{i}\sim \mathcal{N}(\mu,\tau^2) as

data {
    int<lower=0> n;         // number of individuals
    int <lower=0,upper=1> y[n];// binary response
    vector[n] r; ;// continuous variable
  }
  transformed data {
    int<lower=0> N_pos;
    int<lower=1,upper=n> n_pos[sum(y)];
    
    int<lower=0> N_neg;
    int<lower=1,upper=n> n_neg[n - size(n_pos)];
    
    N_pos = size(n_pos);
    N_neg = size(n_neg);
    {
      int i;
      int j;
      i = 1;
      j = 1;
      for (k in 1:n) {
        if (y[k] == 1) {
          n_pos[i] = k;
          i += 1;
        } else {
          n_neg[j] = k;
          j += 1;
        }
      }
    }
  }
  
  parameters {
    real mu;
    real<lower=0> tau;
    vector<upper=0>[N_pos] z_pos;
    vector<lower=0>[N_neg] z_neg;
  }
  
  transformed parameters {
    vector[n] z;
    for (k in 1:N_pos)
      z[n_pos[k]] = z_pos[k] +r[n_pos[k]];
    for (k in 1:N_neg)
      z[n_neg[k]] = z_neg[k] +r[n_neg[k]];
  }
  
  model {
    z ~ normal(mu,tau);
    mu ~ normal(0,5);
    tau ~ cauchy(0, 5) T[0,];

The idea to extend the model for longitudinal data with Z_i common for the same individual

data {
    int<lower=0> n;         // number of individuals
    int<lower=0> N;         // number of individuals x time
    int<lower=0,upper=n> id;// id
    int <lower=0,upper=1> y[N];// binary response
    real r[N]; // continuous response
  }
  transformed data {
    int<lower=0> N_pos;
    int<lower=1,upper=N> n_pos[sum(y)];
    
    int<lower=0> N_neg;
    int<lower=1,upper=N> n_neg[N - size(n_pos)];
    
    N_pos = size(n_pos);
    N_neg = size(n_neg);
    {
      int i;
      int j;
      i = 1;
      j = 1;
      for (k in 1:N) {
        if (y[k] == 1) {
          n_pos[i] = k;
          i += 1;
      } else {
          n_neg[j] = k;
          j += 1;
        }
      }
    }
  }
  
  parameters {
    real mu;
    real<lower=0> tau;
    vector<lower=0>[N_pos] z_pos;
    vector<upper=0>[N_neg] z_neg;
  }
  
  transformed parameters {
    vector[n] z; // n and not N
    for (k in 1:N_pos)
      z[id[n_pos[k]]] = z_pos[k] +r[n_pos[k]];  // how to define z[id[n_pos[k]]]?
    for (k in 1:N_neg)
      z[id[n_neg[k]]] = z_neg[k] +r[n_neg[k]];
  }
  
  model {
    z ~ normal(mu,tau);
    mu ~ normal(0,5);
    tau ~ cauchy(0, 5) T[0,];
  }

The line

z[id[n_pos[k]]] = z_pos[k] +r[n_pos[k]];

was intented to add the constraint Z_i=Z_{i,j} \forall j but it is not working, and it would only redefine Z_i with the last j.

Is there a way to add the constraint Z_i=Z_{i,j} \forall j

Thank you

I think part of the problem is that the model as you have stated it seems weird to me - maybe I am just missing something? If I understand the model correctly, it really just says that Z_i lies in the interval \left(\max\{R_{i,j} | Y_{i,j} = 0\}, \min\{R_{i,j} | Y_{i,j} = 1\} \right) and where exactly in this interval Z_i lies is determined solely by the prior Z_i. This would also have the interesting consequence that not all values of R_{i,j} are allowed, i.e. when Y_{a,b} = 1 and Y_{a,c} = 0 then it follows that R_{a,c} < Z_a < R_{a,b}. Is this correct? If yes, I think the model can be greatly simplified by taking this into account.

Minor notes:
It is generally considered bad practice to have both n and N in the same model as it may cause confusion, consider renaming to e.g. n_individuals. (an exception that seems to be well readable is in for loops of the form for(n in 1:N) { ... } where the relationship between n and N is relatively clear)

The explicit truncation in tau ~ cauchy(0, 5) T[0,]; is not necessary, because when you have bounds on (non-transfromed) parameters, Stan handles the truncation automagically.

Thank you for your answer and your advices!

Yes exactly we have Z_i in the interval (max\{R_{i,j}|Y_{i,j}=0\},min\{R_{i,j}|Y_{i,j}=1\}), and we assume that Z_i follows a normal distribution (priors on the parameters of the normal distribution).

In fact, the R_{i,j} are data, so not all values of Z_i are allowed depending on the data. I don’t know how do add the constraints Z_i \in [max\{R_{i,j}|Y_{i,j}=0\},min\{R_{i,j}|Y_{i,j}=1\}], this is why I added the parameters z_{pos}=Z_i-R_{i,j}\leq 0 for Y_{i,j}=1 and z_{neg}=Z_i-R_{i,j}>0 for Y_{i,j}=1. I corrected some problems in the previous code:

data {
    int<lower=0> n_ind;         // number of individuals
    int<lower=0> N;         // number of individuals x time
    int<lower=0,upper=n_ind> id[N];//  id
    int <lower=0,upper=1> y[N];// binary response
    real r[N];
  }
  transformed data {
    int<lower=0> N_pos;
    int<lower=1,upper=N> n_pos[sum(y)];
    
    int<lower=0> N_neg;
    int<lower=1,upper=N> n_neg[N - size(n_pos)];
    
    N_pos = size(n_pos);
    N_neg = size(n_neg);
    {
      int i;
      int j;
      i = 1;
      j = 1;
      for (k in 1:N) {
        if (y[k] == 1) {
          n_pos[i] = k;
          i += 1;
      } else {
          n_neg[j] = k;
          j += 1;
        }
      }
    }
  }
  
  parameters {
    real mu;
    real<lower=0> tau;
    vector<upper=0>[N_pos] z_pos;
    vector<lower=0>[N_neg] z_neg;
  }
  
  transformed parameters {
    vector[n_ind] z;
    for (k in 1:N_pos)
      z[id[n_pos[k]]] = z_pos[k] +r[n_pos[k]];
    for (k in 1:N_neg)
      z[id[n_neg[k]]] = z_neg[k] +r[n_neg[k]];
  }
  
  model {
    z ~ normal(mu,tau);
    mu ~ normal(0,5);
    tau ~ cauchy(0, 1);
  }

The part

    for (k in 1:N_pos)
      z[id[n_pos[k]]] = z_pos[k] +r[n_pos[k]];
    for (k in 1:N_neg)
      z[id[n_neg[k]]] = z_neg[k] +r[n_neg[k]];
  }

was a try to add the constraints, but it only keeps one constraint, do you have an idea to define the constraint Z_i \in [max\{R_{i,j}|Y_{i,j}=0\},min\{R_{i,j}|Y_{i,j}=1\}]?

I’m short on time, so will just sketch my ideas despite there being some non-trivial stuff to wrap your head around, hope that helps anyway. First, I would use transformed data to compute the bounds, the idea would be to produce something like:

transformed data {
  vector[n_ind] lower_bound; //should be equal to negative_infinity() when no lower bound was found
  vector[n_ind] upper_bound; //should be equal to positive_infinity() when no lower bound was found
  
  //Code to compute this not included :-)
  //Also possibly throw an error via `reject` if some lower bounds are higher/equal to upper bounds
}

Then I would have non-constrained parameters and transform them by hand for the appropriate bounds. You can use the same techniques described in the Stan manual at https://mc-stan.org/docs/2_19/reference-manual/lower-bound-transform-section.html and the following sections. Note that we will need to compute the Jacobian for the variable transform by hand as we will use the transformed variables on the left hand side of a sampling statement (see " User-Transformed Variables" under https://mc-stan.org/docs/2_19/reference-manual/sampling-statements-section.html). Roughly (code not tested, Jacobian computation taken from the manual on bounds, but not double checked with the manual):

parameters {
  real mu;
  real<lower=0> tau;
  vector[n_ind] z_raw;
}

transformed parameters {
  vector[n_ind] z;
  vector[n_ind] log_jacobian;
  for(i in 1:n_ind) {
     if(is_inf(lower_bound[i])) {
       if(is_inf(upper_bound[i])) {
         //No bounds, direct copy
         z[i] = z_raw[i];
         log_jacobian[i] = 0;
       } else {
         //upper bound only
         z[i] = upper_bound[i] - exp(z_raw[i]);
         log_jacobian[i] = z_raw[i];
       }
     } else {
       if(is_inf(upper_bound[i])) {
         //Lower bound only
         z[i] = lower_bound[i] + exp(z_raw[i]);
         log_jacobian[i] = z_raw[i];
       } else {
         //Both bounds
         real inv_logit_z = inv_logit(z_raw[i]);
         z[i] = lower_bound[i] + (upper_bound[i] - lower_bound[i]) * inv_logit_z;
         log_jacobian[i] = 
           log(upper_bound[i] - lower_bound[i])  + log(inv_logit_z)) +
           log(1 - inv_logit_z);
     }
  }
}

model {
    z ~ normal(mu,tau);
    taget += sum(log_jacobian); //Add the log Jacobian to target density
    mu ~ normal(0,5);
    tau ~ cauchy(0, 1);
}