Epred within Stan

Modeling
1

Hi,

I would like to replicate on Stan the brms model example from this excellent post

// Beta regression model 
// Example from: https://www.andrewheiss.com/blog/2021/11/08/beta-regression-guide/

data {
  int<lower=1> N;                      // sample size
  int<lower=1> K;                      // K predictors
  vector<lower=0,upper=1>[N] Y;        // response 
  matrix[N,K] X;                       // predictor matrix
}

parameters {
  vector[K] theta;                     // reg coefficients
  vector<lower=0>[K] phi;                   // dispersion parameter
}

transformed parameters{
  vector[K] beta;
  vector[K] beta_phi;

  beta = theta * 5;                    // same as beta ~ normal(0, 5); fairly diffuse
  beta_phi = phi * 5;
  
}

model {
  // model calculations
  vector[N] LP;                        // linear predictor
  vector[N] LP_phi;
  vector[N] mu;                        // transformed linear predictor
  vector[N] mu_phi; 
  vector[N] A;                         // parameter for beta distn
  vector[N] B;                         // parameter for beta distn

  LP = X * beta;
  LP_phi = X * beta_phi;
  
  for (i in 1:N) { 
    mu[i] = inv_logit(LP[i]);   
  }
  
  for (i in 1:N) { 
    mu_phi[i] = exp(LP_phi[i]);   
  }

  A = mu .* mu_phi;
  B = (1.0 - mu) .* mu_phi;

  // priors
  theta ~ normal(0, 1);   
  phi ~ cauchy(0, 5);                  // different options for phi  
  //phi ~ inv_gamma(.001, .001);
  //phi ~ uniform(0, 500);             // put upper on phi if using this

  // likelihood
  Y ~ beta(A, B);
}

generated quantities{
  vector[N] y_rep;
  
  for (i in 1:N) { 
    real mu;
    real mu_phi;
    real A;
    real B;
    
    mu = inv_logit(X[i] * beta);   
    mu_phi = exp(X[i] * beta_phi);  
    
    A = mu * mu_phi;
    B = (1.0 - mu) * mu_phi;
    
    y_rep[i] = beta_rng(A, B); 

  }
}

The above Stan code (found here) produces similar estimates with the brms model brm( bf(prop_fem ~ quota, phi ~ quota) but the predicted values from the beta_rng are similar to posterior_predict. Instead I would like to get something similar to posterior_epred.

Any ideas of how I can get the posterior_epred within Stan?

Thanks for your time.

2

In this model, the expectation of the posterior predictive distribution is simply mu.

2 Likes