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