Multilevel hierarchical linear regression model with nested predictors

Modeling
1

Hello Everyone,

About my specific model, i’m kind of struggling right now, and i am wondering if maybe you could help me out. Specifically i want to fit a multilevel hierarchical linear regression model that has nested predictors. It’s actually a time series but we do not treat it as so because some of the predictors quantify seasonality effects. There are 3 important levels in our model, and we call them P, C, and R. There are respectively nR, nC and nP (typically ~ 10) possible discrete values at each specific level and in our setting, P is “over” C, and C is “over” R. In the end, I want to get estimates of the regression coefficients at the P level, at the P x C level, and at the P x C x R level. Here’s the simplified hierarchical structure:

image001
image001.png856×276 19.5 KB

I have trouble understanding 1) how we should specify the code for that in Stan and 2) how the input data should look like. Attached is the Stan code we’ve written so far in which you can see the whole structure and the choice of priors. Ifeel like we’re almost there but also that something is missing to make it work.

test_forum.R (1.9 KB)

Any advice would be amazing.

Thanks a lot ,

2

That doesn’t make it not a time series, just not purely a classical AR(k) or whatever.

Things that have to be downloaded are harder to answer. It’s multiple steps for me from the browser to editor and then I have to fish the Stan model out of the R. Here’s the model indented with whitespace removed:


data {
  int<lower=0> nObservation; // The number of observations
  int<lower=1> nRCP  ; //The number of RCP
  int<lower=1> nCP ; // The number of CP
  int<lower=1> nP ; // The number of P
  int<lower=1> RCP_iD[nRCP] ; 
  int<lower=1> CP_iD[nCP] ; 
  matrix [nObservation,nRCP ] X ; // The model matrix
  real Y[nObservation] ; //  variable to be explicit
}
parameters {
  real beta0;
  real<lower=0> sigma ;
  real <lower=0> sigmaM;
  real <lower=0> sigmaP;
  real <lower=0> sigmaC;
  real <lower=0> sigmaR;
}
transformed parameters {
  matrix [nObservation,1] mu;
  real<lower=0> muM;
  vector<lower=0>[nP] muP ;
  vector<lower=0>[nCP] muCP;
  vector<lower=0>[ nRCP] betaRCP ;
  for(k in 1:nP){
    muP[k]= muM + sigmaP;
  }
  for(p in 1:nCP){
    muCP[p] = muP[CP_iD[p]] + sigmaC;
  }
  for(j in 1:nRCP){
    betaRCP[j] = muCP[RCP_iD[j]] + sigmaR;
  }
  for (i in 1:nObservation){
    mu = beta0 + betaRCP[1:nRCP]*X[i,1:nRCP]  ;
  }
}
model {
  muM ~ normal(0,2) ;
  sigmaP ~ cauchy(0,2);
  muP ~ normal(muM,sigmaP) ;
  sigmaC ~cauchy(0,2) ;
  muCP ~ normal(muP,sigmaC) ;
  sigmaR ~ cauchy(0,2)  ;
  betaRCP ~ normal(muCP,sigmaR) ;
  //mu <- beta0 + x*betaRCP ;
  sigma ~ cauchy(0,2) ;
  // likelihood
  for (i in 1:nObservation){
    Y[i]~ normal(mu[i],sigma) ;
  } 
}

Is ther esomething specific you’re asking about here? You seem to have commented out the predictor calculation for the likelihood. You also want to use a non-centered parameterizations for those random effects.