Menu
I've been working on a three-level RStan model where repeated broadband measurements (year ID = yrid) are nested within local authorities (LA ID = laid), which are finally nested within regions (region ID = rnid). The (logged)dependent variable is speed and the (logged) predictors are population density (pd) and superfast broadband penetration (sfbb). Currently there are random intercepts at the local authority and regional level (levels 2 & 3).
How do I extend the model to have a random slope at either level 1 or level 2?
Here's a subset of data, the RStan model and the overall R code. Any help would be much appreciated.
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)
Pulling out the stan code for readability:
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);
}
}
409
Assuming you want random effects at both levels, you just need to do the same thing again. So let's say you have an item y[n], who belongs to first level group gg[n] in 1:G, and that group g belongs to second level group hh[g] in 1:H. Then you just have the random intercept parameters
vector[G] c; vector[H] d;
and then the regression just does the index fiddling
for (n in 1:N) mu[n] <- ... + c[gg[n]] + d[hh[gg[n]];
I'd also strongly recommend not including the Stan program as a string --- it messes up the ability to use print statements (or transpositions) and you lose the line numbers. And indenting to make the blocks more easily scannable.
90
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
