15.2 Modelo programando en STAN

El código se divide en varios bloques:

1. El primer bloque especifica las funciones que se utilizarán en el modelo. En este caso, solo hay una función llamada pred_theta, que toma una matriz de covariables, un número de categorías y una matriz de parámetros beta y devuelve una matriz de parámetros theta.

2. El segundo bloque define los datos de entrada que se utilizarán en el modelo. Se especifica el número de dominios, categorías y regresores, y se proporcionan las matrices de datos observados y covariables observadas.

3. El tercer bloque especifica los parámetros del modelo. En este caso, se incluyen la matriz de parámetros beta, las desviaciones estándar de los efectos aleatorios, la matriz de correlación de los efectos aleatorios y la matriz de efectos aleatorios.

4. El cuarto bloque transforma los parámetros para asegurar que los efectos aleatorios tengan la matriz de correlación deseada y calcula los parámetros theta.

5. El quinto bloque define el modelo. En este caso, se incluye una distribución previa para los parámetros beta, las desviaciones estándar y la matriz de correlación, así como la verosimilitud de la distribución multinomial para los datos observados.

6. El sexto bloque genera las cantidades posteriores a partir de los parámetros estimados. En este caso, se incluye la matriz de parámetros theta para los datos de predicción y la matriz de correlación de los efectos aleatorios.

functions {
  matrix pred_theta(matrix Xp, int p, matrix beta){
  int D1 = rows(Xp);
  real num1[D1, p];
  real den1[D1];
  matrix[D1,p] theta_p;
  
  for(d in 1:D1){
    num1[d, 1] = 1;
    num1[d, 2] = exp(Xp[d, ] * beta[1, ]' ) ;
    num1[d, 3] = exp(Xp[d, ] * beta[2, ]' ) ;
    
    den1[d] = sum(num1[d, ]);
  }
  
  for(d in 1:D1){
    for(i in 2:p){
    theta_p[d, i] = num1[d, i]/den1[d];
    }
    theta_p[d, 1] = 1/den1[d];
   }

  return theta_p  ;
  }
  
}

data {
  int<lower=1> D; // número de dominios 
  int<lower=1> P; // categorías
  int<lower=1> K; // cantidad de regresores
  int hat_y[D, P]; // matriz de datos
  matrix[D, K] X_obs; // matriz de covariables
  int<lower=1> D1; // número de dominios 
  matrix[D1, K] X_pred; // matriz de covariables
}
  

parameters {
  matrix[P-1, K] beta;// matriz de parámetros 
  vector<lower=0>[P-1] sigma_u;       // random effects standard deviations
  // declare L_u to be the Choleski factor of a 2x2 correlation matrix
  cholesky_factor_corr[P-1] L_u;
  matrix[P-1, D] z_u;                  
}

transformed parameters {
  simplex[P] theta[D];// vector de parámetros;
  real num[D, P];
  real den[D];
  // this transform random effects so that they have the correlation
  // matrix specified by the correlation matrix above
  matrix[P-1, D] u; // random effect matrix
  u = diag_pre_multiply(sigma_u, L_u) * z_u;
  
  for(d in 1:D){
    num[d, 1] = 1;
    num[d, 2] = exp(X_obs[d, ] * beta[1, ]' + u[1, d]) ;
    num[d, 3] = exp(X_obs[d, ] * beta[2, ]' + u[2, d]) ;
    
    den[d] = sum(num[d, ]);
  }
  
  for(d in 1:D){
    for(p in 2:P){
    theta[d, p] = num[d, p]/den[d];
    }
    theta[d, 1] = 1/den[d];
  }
}

model {
  L_u ~ lkj_corr_cholesky(1); // LKJ prior for the correlation matrix
  to_vector(z_u) ~ normal(0, 10000);
  sigma_u ~ inv_gamma(0.0001, 0.0001);
  
  for(p in 2:P){
    for(k in 1:K){
      beta[p-1, k] ~ normal(0, 10000);
    }
    }
  
  for(d in 1:D){
    target += multinomial_lpmf(hat_y[d, ] | theta[d, ]); 
  }
}

  
generated quantities {
  matrix[D1,P] theta_pred;
  matrix[2, 2] Omega;
  Omega = L_u * L_u'; // so that it return the correlation matrix
  
 theta_pred = pred_theta(X_pred, P, beta);
}