Regression with random intercept for groups

Hi all,

I am new to using Stan but have some experience with WinBUGS from a few years ago. I am trying to work though the examples in the Marc Kery (2010) book “Introduction to Winbugs for Ecologists”, but am having trouble recoding his example into Stan and would appreciate some help.

The example looks at estimating snake mass as a function of snake length, with population as a random effect. The data is simulated and I am trying to recover the parameters, but am having trouble with sigma for intercept. In addition, I am getting the error message “There were X divergent transition after warm-up”.

Here is the R code

``````n.groups = 56
n.sample = 10
n = n.groups * n.sample
pop = gl(n = n.groups,k = n.sample)

original.length = runif(n,45,70) # body length
mn = mean(original.length)
sd = sd(original.length)
length = (original.length - mn) / sd
hist(length,col="grey")

Xmat = model.matrix(~ pop * length-1-length) ### design matrix without intercept

intercept.mean = 230 # mu_alpha
intercept.sd = 20 # sigma alpha
slope.mean = 60 # mu_beta
slope.sd = 30 # sigma_beta

intercept.effects = rnorm(n = n.groups, mean = intercept.mean, sd = intercept.sd)
slope.effects = rnorm(n = n.groups, mean = slope.mean, sd = slope.sd)
all.effects = c(intercept.effects, slope.effects)

lin.pred = Xmat[,] %*% all.effects
eps = rnorm(n = n,mean = 0,sd=30)
mass = as.vector(lin.pred + eps)
hist(mass,col="grey")

### put data together
data = data.frame(pop=pop,length = length,mass = mass)

##data.list <- list(y = as.vector(mass), X = Xmat, nX = ncol(Xmat), N = nrow(Xmat))
data.list = list(y = as.vector(mass),x = Xmat,pop = as.numeric(pop),length = length,Nobs = length(as.vector(mass)),
Ngroups = n.groups)

fit = stan(file = "longitudinal.stan",data = data.list,chains = 3,
iter = 5000, warmup = 500, thin = 3)

print(fit,pars=c("alpha","beta","sigmaint","sigmaeps"))
``````

Here is the Stan file

``````data {
int<lower=0> Nobs; // Number of observed data
int<lower=0> Ngroups;
int pop[Nobs]; // populations
vector[Nobs] y;
vector[Nobs] length; // new code
}

parameters {
real beta; // common slope
vector[Ngroups] alpha; // random intercept
real<lower=0> sigmaint; // corresponds to mu.int, random intercept sd
real<lower=0> sigmaeps; // estimate error

}

transformed parameters {
vector[Nobs] yhat;

yhat = alpha[pop] + beta * length;

}
model {
//alpha ~ normal(0,sigmaint);
y ~ normal(yhat, sigmaeps);
}
``````

I am sure that I am doing it wrong, but just can’t figure out what I am doing wrong. Any help would be appreciated.

To start with, you are going to need proper priors on the parameters, but it looks as if `sigmaint` and `etaint` are not even being used in your model block. Also, you can write

``````vector[Nobs] yhat;
for (i in 1:Nobs)
yhat[i] = alpha[pop[i]] + beta * length[i];
``````

as

``````vector[Nobs] yhat = alpha[pop] + beta * length;
``````

Thanks,

I rewrote the code based on an example located here.

The model is doing a better job recover the parameters, but it is overestimating sigmaeps, which should be 30.

Here is the code.

``````data {
int<lower=0> Nobs; // Number of observed data
int<lower=0> Ngroups;
int pop[Nobs]; // populations
vector[Nobs] y;
vector[Nobs] length;
}

parameters {
real beta; // common slope
real alpha; // common intercept

// random effect
vector[Ngroups] alpha2; // random intercept
real<lower=0> sigmaint; // corresponds to mu.int, random intercept sd

real<lower=0> sigmaeps; // estimate error
}

transformed parameters {
vector[Nobs] yhat;

// varying intercepts
real alpha3[Ngroups];

// varying intercepts definition
for (k in 1:Ngroups)
alpha3[k] = alpha + alpha2[k];

for (i in 1:Nobs)
yhat[i] = alpha3[pop[i]]+ beta*length[i];

}
model {
// random effects distribution
alpha2 ~ normal(0,sigmaint);
y ~ normal(yhat, sigmaeps);
}``````