code_linear_dependency.stan

data {
  int<lower=1> D; // Dimensions of the features and coefficient vectors
  int<lower=0> N_obs;  // Number of observed observations (with T = 1)
  int<lower=0> N_cens; // Number of censored observations (with T = 0)
  int<lower=0> M; // Number of judges
  real<lower=0> sigma_tau; // Prior for the variance parameters.
  int<lower=1, upper=M> jj_obs[N_obs];   // judge_ID array
  int<lower=1, upper=M> jj_cens[N_cens]; // judge_ID array
  int<lower=0, upper=1> dec_obs[N_obs];   // Decisions for the observed observations
  int<lower=0, upper=1> dec_cens[N_cens]; // Decisions for the censored observations
  row_vector[D] X_obs[N_obs];   // Features of the observed observations 
  row_vector[D] X_cens[N_cens]; // Features of the censored observations 
  int<lower=0, upper=1> y_obs[N_obs]; // Outcomes of the observed observations
}

parameters {
  // Latent variable
  vector[N_obs] Z_obs;
  vector[N_cens] Z_cens;
  
  // Intercepts and their variance parameters.
  real<lower=0> sigma_T;
  real<lower=0> sigma_Y;

  vector[M] alpha_T_raw;
  // vector[M] alpha_T;
  real alpha_Y_raw;

  // Temporary variables to compute the coefficients
  vector[D] a_XT;
  vector[D] a_XY;
  real<lower=0> a_ZT; // Presume latent variable has a positive coefficient.
  real<lower=0> a_ZY;

  real<lower=0> tau_X; // RL 08DEC2019: Combine taus for X.
  real<lower=0> tau_Z; // RL 24OCT2019: Combine taus for Z.

}

transformed parameters {
  
  // Coefficients
  vector[D] beta_XT;
  vector[D] beta_XY;
  real<lower=0> beta_ZT; // Presume latent variable has a positive coefficient.
  real<lower=0> beta_ZY;
  
  // Intercepts
  vector[M] alpha_T;
  real alpha_Y;
 
  beta_XT = a_XT / sqrt(tau_X); // RL 08DEC2019: Combine taus for X.
  beta_XY = a_XY / sqrt(tau_X);
  beta_ZT = a_ZT / sqrt(tau_Z); // RL 24OCT2019: Combine taus for Z.
  beta_ZY = a_ZY / sqrt(tau_Z);
  
  // beta_XT = a_XT * tau_X; // RL 08DEC2019: Combine taus for X.
  // beta_XY = a_XY * tau_X;
  // beta_ZT = a_ZT * tau_Z; // RL 24OCT2019: Combine taus for Z.
  // beta_ZY = a_ZY * tau_Z;
  
  alpha_T = sigma_T * alpha_T_raw;
  alpha_Y = sigma_Y * alpha_Y_raw;

}

model {
  // Latent variable
  Z_obs ~ normal(0, 1);
  Z_cens ~ normal(0, 1);
  
  // Variance parameter for the intercepts
  sigma_T ~ normal(0, sigma_tau);
  sigma_Y ~ normal(0, sigma_tau);
  
  // Intercepts
  alpha_T_raw ~ normal(0, 1);
  // alpha_T ~ logistic(0, 1);
  alpha_Y_raw ~ normal(0, 1);
  
  // According to 
  // https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations
  // section "Prior for the regression coefficients in logistic regression
  // (non-sparse case)" a good prior is Student's with parameters 3 < nu < 7, 
  // 0 and 1.
  
  // Also, according to
  // https://mc-stan.org/docs/2_19/stan-users-guide/reparameterization-section.html
  // reparametrizing beta ~ student_t(nu, 0, 1) can be done by sampling
  // a ~ N(0, 1) and tau ~ gamma(nu/2, nu/2) so then beta = a * tau^(-1/2).
  // Reparametrizing is done to achieve better mixing.

  a_XT ~ normal(0, 1);
  a_XY ~ normal(0, 1);
  a_ZT ~ normal(0, 1);
  a_ZY ~ normal(0, 1);

  // nu = 6 -> nu/2 = 3
  tau_X ~ gamma(3, 3);
  tau_Z ~ gamma(3, 3);
  // tau_X ~ normal(0, 1);
  // tau_Z ~ normal(0, 1);

  // Compute the regressions for the observed observations
  for(i in 1:N_obs){
    dec_obs[i] ~ bernoulli_logit(alpha_T[jj_obs[i]] + X_obs[i] * beta_XT + beta_ZT * Z_obs[i]);
    y_obs[i] ~ bernoulli_logit(alpha_Y + X_obs[i] * beta_XY + beta_ZY * Z_obs[i]);
  }

  // Compute the regression for the censored observations
  for(i in 1:N_cens)
    dec_cens[i] ~ bernoulli_logit(alpha_T[jj_cens[i]] + X_cens[i] * beta_XT + beta_ZT * Z_cens[i]);
}

generated quantities {
  // Array for predictions
  int<lower=0, upper=1> y_est[N_cens];

  // Generate a draw from the posterior predictive for the outcome.
  for(i in 1:N_cens){
    y_est[i] = bernoulli_logit_rng(alpha_Y + X_cens[i] * beta_XY + beta_ZY * Z_cens[i]);
  }
}