Я работал над трехуровневой моделью RStan, в которой повторяющиеся широкополосные измерения (идентификатор года = yrid) вложены в местные органы власти (идентификатор LA = lay), которые, наконец, вложены в регионы (идентификатор региона = 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)
Вытягивание кода стана для удобочитаемости:
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);
}
}
