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.

Wade