2

Rstan で 3 レベルのネストされた線形モデルを実行しようとしていますが、エラー メッセージが表示され続けます。

この 3 レベルのネストされたモデルからインスピレーションを得ました。

http://rstudio-pubs-static.s3.amazonaws.com/64315_bc3a395edd104095a8384db8d9952f43.html

本質的に、私は 174 の地方自治体内にネストされ、最終的に 11 の地域内にネストされた、4 つの別々の年 (合計 696 の観察) からの平均ブロードバンド速度の測定データを繰り返しました。

データは次のとおりです。

yrid laid rnid speed
1     1    1     2.3
2     1    1     8.4
3     1    1     9.2
4     1    1     5.8
5     2    1     1.5
6     2    1     9.8
7     2    1     4.7

コード：

#load and prepare data
total<-read.csv("total.csv")

names(total)[names(total)=="X"] <- "yrid"
names(total)[names(total)=="lad.num"] <- "laid"
names(total)[names(total)=="region"] <- "rnid"

## Sort by laid for easier handling in Stan
total <- arrange(total, 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"]

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

## Combine as a stan dataset
dat <- (list(Ni           = length(unique(Ni)),
             Nj           = length(unique(Nj)),
             Nk           = length(unique(Nk)),
             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;
   ## 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
   ## Population intercept (a real number)
   real beta_0;
   ## Population slope
   ## real beta_1;
   ## 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_0 + 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]];
   }
   }
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);
   }
   }"
#run model
   resStan <-stan(model_code = stanmodelcode, data=dat, 
   iter=3000,    chains=3,   warmup=500, thin=10)}

エラー：

COMPILING THE C++ CODE FOR MODEL 'stanmodelcode' NOW.
Error : mismatch in dimension declared and found in context; 
processing stage=data initialization; variable name=laid; position=0; 
dims declared=(1); dims found=(696) failed to create the sampler; 
sampling not done
4

1 に答える 1

1

エラー メッセージには、変数の次元が 1 であると想定されていることlaidが示されていますが、実際に stan に渡されるデータの次元は 696 です。

あなたのコードでは

dat <- (list(Ni           = length(unique(Ni)),

つまり、Ni の一意の値の長さの値を Ni に割り当てます。ニとは何？あなたが持っている

Ni <- length(unique(total$yrid))

私が信じているのは 696 です。したがって、エラーの原因である がNi = length(unique(696))得られます。Ni = 1

dat変数を次のように変更してみてください。

dat <- (list(Ni           = Ni,
             Nj           = Nj,
             Nk           = Nk,
             laid         = laid,
             rnid         = rnid,
             regionLookup = regionLookupVec,
             speed        = speed))
于 2015-09-03T18:16:00.223 に答える