Значение параметров регрессии BRM

Question

Я использую пакет brms для построения многоуровневой модели с гауссовским процессом на предикторе x. Модель выглядит следующим образом: make_stancode(y ~ gp(x, cov = "exp_quad", by= groups) + (1| groups), data = dat), поэтому gp предиктора x и многоуровневая групповая переменная. В моем случае у меня 5 групп. Я искал код для этого (ниже), и я пытаюсь выяснить значения и размеры некоторых параметров.

Я вижу, что M_1 это количество групп

Мои вопросы:

1. Что означает значение N_1, равно ли оно количеству наблюдений, N? Здесь используется: vector[N_1] z_1[M_1]; // немасштабированные эффекты на уровне группы

2. Для Kgp_1 и Mgp_1 ( int Kgp_1; и int Mgp_1;), если у меня есть 5 групп, оба Kgp_1 и Mgp_1 равны 5? Если так, почему используются две переменные?

   // генерируется с функциями brms 1.10.0 {

     /* compute a latent Gaussian process
      * Args:
      *   x: array of continuous predictor values
      *   sdgp: marginal SD parameter
      *   lscale: length-scale parameter
      *   zgp: vector of independent standard normal variables 
      * Returns:  
      *   a vector to be added to the linear predictor
      */ 
     vector gp(vector[] x, real sdgp, real lscale, vector zgp) { 
       matrix[size(x), size(x)] cov;
       cov = cov_exp_quad(x, sdgp, lscale);
       for (n in 1:size(x)) {
         // deal with numerical non-positive-definiteness
         cov[n, n] = cov[n, n] + 1e-12;
       }
       return cholesky_decompose(cov) * zgp;
     }
   } 
   data { 
     int<lower=1> N;  // total number of observations 
     vector[N] Y;  // response variable 
     int<lower=1> Kgp_1; 
     int<lower=1> Mgp_1; 
     vector[Mgp_1] Xgp_1[N]; 
     int<lower=1> Igp_1[Kgp_1]; 
     int<lower=1> Jgp_1_1[Igp_1[1]]; 
     int<lower=1> Jgp_1_2[Igp_1[2]]; 
     int<lower=1> Jgp_1_3[Igp_1[3]]; 
     int<lower=1> Jgp_1_4[Igp_1[4]]; 
     int<lower=1> Jgp_1_5[Igp_1[5]]; 
     // data for group-level effects of ID 1 
     int<lower=1> J_1[N]; 
     int<lower=1> N_1; 
     int<lower=1> M_1; 
     vector[N] Z_1_1; 
     int prior_only;  // should the likelihood be ignored? 
   } 
   transformed data { 
   } 
   parameters { 
     real temp_Intercept;  // temporary intercept 
     // GP hyperparameters 
     vector<lower=0>[Kgp_1] sdgp_1; 
     vector<lower=0>[Kgp_1] lscale_1; 
     vector[N] zgp_1; 
     real<lower=0> sigma;  // residual SD 
     vector<lower=0>[M_1] sd_1;  // group-level standard deviations 
     vector[N_1] z_1[M_1];  // unscaled group-level effects 
   } 
   transformed parameters { 
     // group-level effects 
     vector[N_1] r_1_1 = sd_1[1] * (z_1[1]); 
   } 
   model { 
     vector[N] mu = rep_vector(0, N) + temp_Intercept; 
     mu[Jgp_1_1] = mu[Jgp_1_1] + gp(Xgp_1[Jgp_1_1], sdgp_1[1], lscale_1[1], zgp_1[Jgp_1_1]); 
     mu[Jgp_1_2] = mu[Jgp_1_2] + gp(Xgp_1[Jgp_1_2], sdgp_1[2], lscale_1[2], zgp_1[Jgp_1_2]); 
     mu[Jgp_1_3] = mu[Jgp_1_3] + gp(Xgp_1[Jgp_1_3], sdgp_1[3], lscale_1[3], zgp_1[Jgp_1_3]); 
     mu[Jgp_1_4] = mu[Jgp_1_4] + gp(Xgp_1[Jgp_1_4], sdgp_1[4], lscale_1[4], zgp_1[Jgp_1_4]); 
     mu[Jgp_1_5] = mu[Jgp_1_5] + gp(Xgp_1[Jgp_1_5], sdgp_1[5], lscale_1[5], zgp_1[Jgp_1_5]); 
     for (n in 1:N) { 
       mu[n] = mu[n] + (r_1_1[J_1[n]]) * Z_1_1[n]; 
     } 
     // priors including all constants 
     target += student_t_lpdf(sdgp_1 | 3, 0, 10)
       - 1 * student_t_lccdf(0 | 3, 0, 10); 
     target += normal_lpdf(lscale_1 | 0, 0.5)
       - 1 * normal_lccdf(0 | 0, 0.5); 
     target += normal_lpdf(zgp_1 | 0, 1); 
     target += student_t_lpdf(sigma | 3, 0, 10)
       - 1 * student_t_lccdf(0 | 3, 0, 10); 
     target += student_t_lpdf(sd_1 | 3, 0, 10)
       - 1 * student_t_lccdf(0 | 3, 0, 10); 
     target += normal_lpdf(z_1[1] | 0, 1); 
     // likelihood including all constants 
     if (!prior_only) { 
       target += normal_lpdf(Y | mu, sigma); 
     } 
   } 
   generated quantities { 
     // actual population-level intercept 
     real b_Intercept = temp_Intercept; 
   } 
   

Answer 1

Если вы используете make_standata(...) по той же формуле вы можете увидеть данные, которые будут переданы в Stan. Отсюда вы можете собрать воедино некоторые переменные. Если я использую lme4::sleepstudy Набор данных в качестве прокси для ваших данных, я получаю:

library(brms)
dat <- lme4::sleepstudy
dat$groups <- dat$Subject
dat$y <- dat$Reaction
dat$x <- dat$Days

s_data <- make_standata(
  y ~ gp(x, cov = "exp_quad", by= groups) + (1| groups), data = dat)
s_data$N_1
#> 18 

За N_1Я получаю 18, что является количеством уровней в groups в этом наборе данных.

Для Kgp_1 и Mgp_1 ( int Kgp_1; и int Mgp_1;), если у меня есть 5 групп, оба Kgp_1 и Mgp_1 равны 5? Если так, почему используются две переменные?

s_data$Mgp_1
#> 1
s_data$Kgp_1
#> 18

Это выглядит как Kgp_1 опять количество групп. Я не уверен, что делает Mgp_1, кроме того, чтобы установить длину вектора vector[Mgp_1] Xgp_1[N];