RStan:指定三级随机斜率模型？

Question

我一直在研究一个三级 RStan 模型，其中重复的宽带测量(年份 ID = yrid)嵌套在地方当局(LA ID = 铺设)内，最终嵌套在区域内(区域 ID = rnid)。 (记录的)因变量是速度，(记录的)预测变量是人口密度 (pd) 和超高速宽带渗透率 (sfbb)。目前在地方当局和区域级别(2 级和 3 级)有随机拦截。

如何扩展模型以在 1 级或 2 级具有随机斜率？

这是数据的子集、RStan 模型和整个 R 代码。任何帮助将非常感激。

 library(rstan)
 ###Data
  yrid = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
     21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,
     41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,
     61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,
     81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,
     101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,
     121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,
     141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,
     161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,
     181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199)
  laid <- c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,
     6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10,
     11,11,11,11,12,12,12,12,13,13,13,13,14,14,14,14,15,15,15,15,
     16,16,16,16,17,17,17,17,18,18,18,18,19,19,19,19,20,20,20,20,
     21,21,21,21,22,22,22,22,23,23,23,23,24,24,24,24,25,25,25,25,
     26,26,26,26,27,27,27,27,28,28,28,28,29,29,29,29,30,30,30,30,
     31,31,31,31,32,32,32,32,33,33,33,33,34,34,34,34,35,35,35,35,
     36,36,36,36,37,37,37,37,38,38,38,38,39,39,39,39,40,40,40,40,
     41,41,41,41,42,42,42,42,43,43,43,43,44,44,44,44,45,45,45,45,
     46,46,46,46,47,47,47,47,48,48,48,48,49,49,49,49,50,50,50)
 rnid <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,
     3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,
     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
     4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
     5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
     6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,
     7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8)
pd <- c(7.59262,7.59875,7.6027,7.60375,7.5301,7.53444,7.53604,7.54136,8.378,8.3936,
   8.40061,8.41183,7.36682,7.36992,7.37607,7.38268,7.20065,7.2011,7.20162,7.20578,
   7.78846,7.79947,7.80743,7.81992,7.71797,7.72011,7.72396,7.73026,7.66336,7.66561,
   7.66744,7.66833,7.66973,7.67587,7.68327,7.69321,7.4326,7.43449,7.43762,7.44167,
   7.43053,7.43053,7.43189,7.43396,8.33459,8.34315,8.34548,8.35036,7.15921,7.16325,
   7.16379,7.16943,7.4898,7.48869,7.48689,7.48796,7.61918,7.62046,7.62075,7.62261,
   6.55763,6.56541,6.57438,6.58286,6.27777,6.27833,6.28133,6.28339,6.80184,6.8045,
   6.80572,6.81113,7.31315,7.32324,7.32804,7.33446,7.24893,7.24843,7.24744,7.24993,
   7.80751,7.81927,7.83475,7.84514,7.80045,7.80147,7.80543,7.80792,7.74119,7.74253,
   7.74323,7.74457,7.6027,7.6042,7.60564,7.60852,8.29695,8.30721,8.31356,8.32186,
   8.07527,8.09465,8.11516,8.13795,8.06994,8.07091,8.07347,8.07788,8.19141,8.19883,
   8.20841,8.21603,7.05652,7.05893,7.06613,7.07089,7.85991,7.86511,7.8699,7.87721,
   8.18894,8.19332,8.19572,8.20125,7.26382,7.26669,7.2701,7.27351,6.32972,6.33505,
   6.34036,6.34529,6.94235,6.94832,6.95483,6.96111,7.21575,7.22504,7.23006,7.23648,
   6.87109,6.87472,6.8811,6.88623,6.89163,6.89264,6.89811,6.897,7.85077,7.85294,
   7.85438,7.85582,6.31409,6.31264,6.31192,6.31264,6.84119,6.84428,6.84843,6.85309,
   6.28171,6.27796,6.27983,6.27983,5.43764,5.44025,5.44328,5.44674,4.14155,4.14155,
   4.13996,4.14155,5.76142,5.76519,5.76676,5.77082,7.37092,7.37092,7.37331,7.37651,
   7.02322,7.02811,7.035,7.04132,5.88444,5.88666,5.88915,5.89275,6.98296,6.98296,
   6.98091,6.97616,8.31179,8.3111,8.30687,8.30048,8.18502,8.1893,8.19085)
 sfbb <- c(4.41884,4.47506,4.53903,4.52179,4.49981,4.51196,4.55071,4.48864,4.35671,4.39938,
     4.46245,4.46591,4.47734,4.54966,4.57883,4.54329,4.06044,4.40305,4.56643,4.5326,
     4.31749,4.44265,4.51196,4.48864,4.54329,4.56643,4.59107,4.55388,4.45435,4.52287,
     4.5401,4.49981,4.46591,4.51525,4.53903,4.49981,4.46591,4.4613,4.54116,4.51086,
     4.33073,4.33598,4.45435,4.45435,4.34381,4.41159,4.45435,4.45435,4.41884,4.42365,
     4.46476,4.40672,4.21951,4.32281,4.50756,4.51086,4.09434,4.12552,4.26127,4.34381,
     4.27667,4.18662,4.25277,4.26268,4.00733,4.10264,4.33205,4.34381,4.00733,4.01638,
     4.28497,4.33073,4.2485,4.27528,4.39815,4.33073,4.21951,4.23989,4.2683,4.31749,
     4.14313,4.25277,4.31615,4.41884,4.07754,4.21951,4.37701,4.40672,4.34381,4.34381,
     4.50976,4.48864,4.00733,4.27528,4.51305,4.51086,4.45435,4.4762,4.51961,4.5326,
     4.18965,4.21951,4.46476,4.49981,4.52179,4.51525,4.56017,4.54329,4.35671,4.34899,
     4.38701,4.44265,4.48864,4.48639,4.51086,4.51086,4.46591,4.49088,4.51852,4.49981,
     4.5326,4.53475,4.57471,4.56435,4.39445,4.45783,4.52721,4.46591,4.02535,4.0656,
     4.20469,4.11087,4.04305,4.27805,4.37952,4.38203,4.40672,4.38203,4.4613,4.44265,
     3.98898,4.29046,4.43912,4.40672,4.47734,4.52829,4.53582,4.51086,4.51086,4.5326,
     4.55808,4.52179,4.46591,4.47734,4.50866,4.46591,4.43082,4.48639,4.51196,4.47734,
     4.41884,4.41643,4.46014,4.47734,3.09104,4.0448,4.13517,4.17439,2.19722,3.59731,
     3.90399,4.15888,4.20469,4.38701,4.4128,4.34381,4.18965,4.19117,4.42843,4.44265,
     4.33073,4.34251,4.37827,4.44265,3.7612,4.08092,4.1987,4.2485,4.12713,4.12228,
     4.22098,4.30407,4.14313,4.13677,4.29456,4.45435,-2.30259,1.52606,2.49321)
   speed <- c(2.10413,2.76632,2.95491,3.29953,1.96009,2.57261,2.81541,3.11795,2.10413,2.56495,
      2.8792,3.17805,1.94591,2.7213,2.96011,3.30689,1.93152,2.38876,2.82731,3.19867,
      2.20827,2.82731,3.03495,3.35341,2.17475,2.82138,3.10459,3.44362,2.04122,2.5416,
      2.83321,3.15274,2.14007,2.68102,3.06805,3.36384,2.15176,2.63189,3.03495,3.36384,
      2.05412,2.49321,2.91235,3.28091,2.25129,2.67415,2.96011,3.32504,2.07944,2.55723,
      2.98568,3.31054,2.17475,2.61007,2.83321,3.2884,2.11626,2.5416,2.79728,3.15274,
      1.93152,2.5337,2.70805,3.02529,1.97408,2.52573,2.7213,3.06805,1.97408,2.4248,
      2.63189,2.97553,1.97408,2.52573,2.73437,3.03975,2.17475,2.57261,2.96527,3.22684,
      2.17475,2.68785,2.93386,3.22684,2.09186,2.71469,2.91777,3.24649,2.10413,2.58022,
      3.03495,3.39786,1.93152,2.38876,2.89037,3.21888,2.18605,2.70805,3.03013,3.36038,
      2.11626,2.50144,2.89591,3.2308,2.14007,2.59525,3.03013,3.36384,2.16332,2.61007,
      2.9069,3.20275,2.11626,2.83321,3.09104,3.40784,2.10413,2.59525,2.96011,3.26957,
      2.20827,2.95491,3.10009,3.46574,2.11626,2.58776,2.94969,3.24649,1.91692,2.50144,
      2.63906,2.97041,1.94591,2.451,2.78501,3.11795,2.06686,2.5416,2.9069,3.2308,
      1.91692,2.41591,2.66026,2.98568,2.06686,2.93386,3.14415,3.47816,2.19722,2.91777,
      3.15274,3.44362,2.15176,2.8679,3.10459,3.421,2.11626,2.74084,3.13983,3.45947,
      2.05412,2.57261,3.03013,3.38439,1.84055,2.24071,2.4681,2.74727,1.84055,2.11626,
      2.35138,2.66723,1.90211,2.50144,2.76001,3.054,2.00148,2.41591,2.8679,3.24259,
      2.0149,2.5096,2.91235,3.26194,1.90211,2.37024,2.6174,2.92316,1.98787,2.52573,
      2.69463,3.05871,2.25129,2.80336,2.90142,3.2581,1.98787,2.2192,2.34181)


 total <- data.frame(speed, pd, sfbb, yrid, laid, rnid)

 ## Create a vector of school IDs where j-th element gives school ID for class ID j
 regionLookupVec <- unique(total[c("laid","rnid")])[,"rnid"]

 ## Design matrix for model 
 desMat <- model.matrix(object = ~ 1 + pd + sfbb , data = total)

 ## Combine as a stan dataset
 Ni <- length(unique(total$yrid))
 Nj <- length(unique(total$laid))
 Nk <- length(unique(total$rnid))              
 p  <- ncol(desMat)
 desMat <-(desMat)
 laid <- (total$laid)
 rnid <- (total$rnid)
 regionLookup <- (regionLookupVec)
 speed <- (total$speed)

 ## Combine as a stan dataset
 dat <- (list(    Ni           = Ni,
                  Nj           = Nj,
                  Nk           = Nk,
                  p            = p,
                  desMat       = desMat,
                  laid         = laid,
                  rnid         = rnid,
                  regionLookup = regionLookupVec,
                  speed        = speed))
 -------------------------------------------------------------------------------------
   #load model
   stanmodelcode <- "data {
 ##Define variables in data
 ##Number of level-1 observations (an integer)
 int<lower=0> Ni;
 ## Number of level-2 clusters
 int<lower=0> Nj;
 ## Number of level-3 clusters
 int<lower=0> Nk;
 ##Number of fixed effect parameters
 int<lower=0> p;
 // Design matrix
 real desMat[Ni,p];
 ## Cluster IDs
 int<lower=1> laid[Ni];
 int<lower=1> rnid[Ni];
 ## Level 3 look up vector for level 2
 int<lower=1> regionLookup[Nj];
 ## Continuous outcome
 real speed[Ni];
 ## Continuous predictor
 ## real X_1ijk[Ni];
 }
 parameters {
 ## Define parameters to estimate
 ## Fixed effects (a real number)
 real beta[p];
 ## Level-1 errors
 real<lower=0> sigma_e0;
 ## Level-2 random effect
 real u_0jk[Nj];
 real<lower=0> sigma_u0jk;
 ## Level-3 random effect
 real u_0k[Nk];
 real<lower=0> sigma_u0k;
 }

 transformed parameters  {
 ## Varying intercepts
 real beta_0jk[Nj];
 real beta_0k[Nk];
 ## Individual mean
 real mu[Ni];
 ## Varying intercepts definition
 ## Level-3 (level-3 random intercepts)
 for (k in 1:Nk) {
 beta_0k[k] <- beta[1] + u_0k[k];
 }
 ## Level-2 (level-2 random intercepts)
 for (j in 1:Nj) {
 beta_0jk[j] <- beta_0k[regionLookup[j]] + u_0jk[j];
 }
 ## Individual mean
 for (i in 1:Ni) {
 mu[i] <- beta_0jk[laid[i]] +
 desMat[i,2]*beta[2] + 
 desMat[i,3]*beta[3];
 }
 }
 model {
 ## Prior part of Bayesian inference
 ## Flat prior for mu (no need to specify if non-informative)
 ## Random effects distribution
 u_0k  ~ normal(0, sigma_u0k);
 u_0jk ~ normal(0, sigma_u0jk);
 ## Likelihood part of Bayesian inference
 ## Outcome model N(mu, sigma^2) (use SD rather than Var)
 for (i in 1:Ni) {
 speed[i] ~ normal(mu[i], sigma_e0);
 }
 }"
-------------------------------------------------------------------------------------

 resStan <-stan(model_code = stanmodelcode, data=dat, iter=100, chains=2, warmup=10, 
                  thin=1, control=list(adapt_delta = 0.9), verbose = TRUE)  

拉出 stan 代码以提高可读性:

data {
     ##Define variables in data
     ##Number of level-1 observations (an integer)
     int<lower=0> Ni;
     ## Number of level-2 clusters
     int<lower=0> Nj;
     ## Number of level-3 clusters
     int<lower=0> Nk;
     ##Number of fixed effect parameters
     int<lower=0> p;
     // Design matrix
     real desMat[Ni,p];
     ## Cluster IDs
     int<lower=1> laid[Ni];
     int<lower=1> rnid[Ni];
     ## Level 3 look up vector for level 2
     int<lower=1> regionLookup[Nj];
     ## Continuous outcome
     real speed[Ni];
     ## Continuous predictor
     ## real X_1ijk[Ni];
   }
   parameters {
     ## Define parameters to estimate
     ## Fixed effects (a real number)
     real beta[p];
     ## Level-1 errors
     real<lower=0> sigma_e0;
     ## Level-2 random effect
     real u_0jk[Nj];
     real<lower=0> sigma_u0jk;
     ## Level-3 random effect
     real u_0k[Nk];
     real<lower=0> sigma_u0k;
   }

   transformed parameters  {
     ## Varying intercepts
     real beta_0jk[Nj];
     real beta_0k[Nk];
     ## Individual mean
     real mu[Ni];
     ## Varying intercepts definition
     ## Level-3 (level-3 random intercepts)
     for (k in 1:Nk) {
       beta_0k[k] <- beta[1] + u_0k[k];
     }
     ## Level-2 (level-2 random intercepts)
     for (j in 1:Nj) {
       beta_0jk[j] <- beta_0k[regionLookup[j]] + u_0jk[j];
     }
     ## Individual mean
     for (i in 1:Ni) {
       mu[i] <- beta_0jk[laid[i]] +
       desMat[i,2]*beta[2] + 
       desMat[i,3]*beta[3];
     }
   }
   model {
     ## Prior part of Bayesian inference
     ## Flat prior for mu (no need to specify if non-informative)
     ## Random effects distribution
     u_0k  ~ normal(0, sigma_u0k);
     u_0jk ~ normal(0, sigma_u0jk);
     ## Likelihood part of Bayesian inference
     ## Outcome model N(mu, sigma^2) (use SD rather than Var)
     for (i in 1:Ni) {
       speed[i] ~ normal(mu[i], sigma_e0);
     }
   }

Answer 1

假设你想要两个级别的随机效果，你只需要再次做同样的事情。因此，假设您有一个项目 y[n]，它属于 1:G 中的第一级组 gg[n]，而该组 g 属于 1:H 中的第二级组 hh[g]。然后你只有随机截取参数

向量[G] c;
向量[H] d;

然后回归只是做索引摆弄

for (n in 1:N)
mu[n] <- ... + c[gg[n]] + d[hh[gg[n]];

我还强烈建议不要将 Stan 程序作为字符串包含在内——它会破坏使用打印语句(或换位)的能力，并且会丢失行号。并缩进使块更易于扫描。