data{ int N; vector[N] x; int y[N]; int is_missing[N]; int J; int K; vector[K] thresh; } parameters{ real alpha; real alpha_miss; real beta; real beta_miss; real sigma; } model{ vector[J] theta; // vector to hold probabilities of ordinal categories vector[N] eta; vector[N] pi; vector[3] logP_y_eta; real logPy; // priors alpha ~ normal(0, 5); alpha_miss ~ normal(0, 1); beta ~ normal(0, 1); beta_miss ~ normal(0, 1); sigma ~ student_t(3, 0, 1); // the likelihood statement for(n in 1:N){ // first, model the ordinal category probabilities from the observed data eta[n] = alpha + beta * x[n]; theta[1] = normal_cdf(thresh[1], eta[n], sigma); for( l in 2:K ){ theta[l] = fmax(0, normal_cdf(thresh[l], eta[n], sigma) - normal_cdf(thresh[l-1], eta[n], sigma) ); } theta[J] = 1 - normal_cdf(thresh[K], eta[n], sigma); // then, compute the likelihood if y is observed if(y[n] > 0){ y[n] ~ categorical(theta); } else{ // we want the log probabilities of each ordinal category given eta // i.e. P(y | eta) = P(y, eta) / P(eta) // here is log(P(y,eta)) for each category logP_y_eta[1] = log(theta[1]) + categorical_lpmf(1 | theta); // P(y=1)P(y=1|eta) logP_y_eta[2] = log(theta[2]) + categorical_lpmf(2 | theta); // P(y=2)P(y=2|eta) logP_y_eta[3] = log(theta[3]) + categorical_lpmf(3 | theta); // P(y=3)P(y=3|eta) // now we marginalise over the log probabilities to get P(y) logPy = log_sum_exp(logP_y_eta); // add it to target target += logPy; } // finally, we want to add a missing data model to explore systematic missingness pi[n] = alpha_miss + beta_miss * x[n]; is_missing[n] ~ bernoulli_logit( pi[n] ); } } generated quantities{ matrix[N,J] y_ord_impute_probs; vector[N] y_ord_impute; vector[N] is_missing_hat; vector[J] theta; vector[N] eta; vector[N] pi; vector[3] logP_y_eta; real logPy; for(n in 1:N){ // first, model the ordinal category probabilities eta[n] = alpha + beta * x[n]; theta[1] = normal_cdf(thresh[1], eta[n], sigma); for( l in 2:K ){ theta[l] = fmax(0, normal_cdf(thresh[l], eta[n], sigma) - normal_cdf(thresh[l-1], eta[n], sigma) ); } theta[J] = 1 - normal_cdf(thresh[K], eta[n], sigma); // then, compute the likelihood // if y is observed if(y[n] > 0){ y_ord_impute_probs[n, 1:3] = [0,0,0]; // set probs to 0 y_ord_impute_probs[n, y[n]] = 1; // set probs to 1 y_ord_impute[n] = y[n]; } else{ // we want the log probabilities of each ordinal category given eta // i.e. P(y | eta) = P(y, eta) / P(eta) logP_y_eta[1] = log(theta[1]) + categorical_lpmf(1 | theta); logP_y_eta[2] = log(theta[2]) + categorical_lpmf(2 | theta); logP_y_eta[3] = log(theta[3]) + categorical_lpmf(3 | theta); // now we marginalise over the log probabilities to get P(y) logPy = log_sum_exp(logP_y_eta); // now we compute P(y,eta) / P(eta) y_ord_impute_probs[n, 1] = exp( logP_y_eta[1] - logPy ); y_ord_impute_probs[n, 2] = exp( logP_y_eta[2] - logPy ); y_ord_impute_probs[n, 3] = exp( logP_y_eta[3] - logPy ); y_ord_impute[n] = categorical_rng(y_ord_impute_probs[n,1:3]'); } // finally, we want to add a missing data model pi[n] = alpha_miss + beta_miss * x[n]; is_missing_hat[n] = bernoulli_logit_lpmf( is_missing[n] | pi[n] ); } }