Regression with random intercept for groups

Interfaces RStan
1

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

2

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;
3

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);
}