ANOVA with random effects in model matrix format

Modeling
1

Hi all,

I am trying to work through Kery’s (2010) book “Introduction to WinBUGS for Ecologists” by coding everything in matrix formulation in STAN. I am having trouble with Chapter 9, an ANOVA with random effects. I can code a model without the model matrix formulation using code borrowed from Faraway, but can’t quite get it in model format. Below is my code. I am receiving errors that " Rows of left-hand-side (150) and rows of right-hand-side (10) must match in size" at line 18. I understand that it is a matrix multiplication problem, but I just can’t see a path forward. Any help would be appreciated.

library(rstan)
##options
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores(logical = FALSE))
chains = 4
iter = 1000

## set parameters
## MCMC settings
ni <- 2000
nt <- 1
nb <- 1000
nc <- 4

## data generation
npop = 10
nsample = 15
n = npop * nsample
pop.grand.mean = 50
pop.sd = 5
pop.means = rnorm(n=npop,mean=pop.grand.mean,sd=pop.sd)
sigma = 3
eps = rnorm(n,0,sigma)

x = rep(1:npop,rep(nsample,npop))
X = as.matrix(model.matrix(~as.factor(x) - 1)) ### This is fit without intercept
y = as.numeric(X %*% as.matrix(pop.means) + eps)

boxplot(y~x,col="grey",xlab = "Population",ylab="SVL",main="",las=1)
abline(h=pop.grand.mean)

# first firt without matrix notation
data.list = list(N = length(x),J=npop,predictor=as.numeric(x),response = y)

anova_fixed.stan = '
data {
  int<lower=0> N; // number of samples
  int<lower=0> J; // number of random groups
  int<lower=1,upper=J> predictor[N]; // groups as random
  vector[N] response; // response variable
}
parameters {
  vector[J] eta; // random effect
  real mu; // overall mean
  real<lower=0> sigmaalpha; // sd on alpha
  real<lower=0> sigmaepsilon; // residual sd.
}
transformed parameters {
  vector[J] a;
  vector[N] yhat;
  
  a = mu + sigmaalpha * eta;
  
  for (i in 1:N)
    yhat[i] = a[predictor[i]];
}
model {
  eta ~ normal(0, 1);
  
  response ~ normal(yhat, sigmaepsilon);
}'

fit = stan(file = "ANOVA_random_one.stan", data = data.list, iter = ni,
           warmup = nb, chains = nc, thin = nt, refresh = 0)

print(fit)

### now as model matrix
Xmat = model.matrix(~as.factor(x),constrast.arg=(x = "contr.sum"))
data.list = list(N = length(x),Npreds=ncol(Xmat),y = y, Xmat = Xmat)

anova_random.stan = '
data {
  int<lower=0> N; // number of samples
  int<lower=0> Npreds; // number of columns, or number of populations
  vector[N] y;
  matrix [N,Npreds] Xmat; // model matrix;
}
parameters {
  vector[Npreds] beta; // betas
  vector[Npreds] eta; // sigma
  real<lower=0> sigmaalpha; // sd on overall mean
  real<lower=0> sigmaepsilon; // residual sd
}
transformed parameters {
  vector[N] mu;
  vector[N] yhat;

  mu = Xmat*beta;	
  yhat = beta[1] + sigmaalpha * eta;
}
model {
  eta ~ normal(0, 1);
  
  y ~ normal(yhat, sigmaepsilon);
}
'

fit2 = stan(file = "OneWay_ANOVA_rand.stan", data = data.list, iter = ni,
            warmup = nb, chains = 1, thin = nt, refresh = 0)

print(fit)
2

Which line is the error for? Line 18 or so I only see a vector times a scalar but that’ doesn’t match your error message. Could we get the complete error message?

3

Hi sakrejda,

Thanks for the reply. The code below worked for me. The problem in my initial code was that in the matrix model I was treating population as a fixed effect and not as a random effect.

anova_random.stan = '
  data {
    int<lower=0> N; // number of samples
    int<lower=0> Npreds; // number of predictors
    int<lower=1> Ngroups; // number of random groups
    int<lower=1,upper=Ngroups> group[N]; // group identity
    vector[N] y;
    matrix [N,Npreds] Xmat; // model matrix;
  }
  parameters {
    vector[Npreds] beta; // betas
    vector[Ngroups] eta; // random effect
    real<lower=0> sigmaeps; // residual sd
    real<lower=0> sigmaalpha; // sd on overall mean
    
  }
  transformed parameters {
    vector[N] yhat;
    
    for (i in 1:N)
      yhat[i] = Xmat[i]*beta + eta[group[i]];
    
  }
  model {
    eta ~ normal(0, sigmaalpha);
    sigmaeps ~ cauchy(0, 5);
    y ~ normal(yhat, sigmaeps);
  }
'
Xmat = model.matrix(lmer(y~ 1 + (1|x)))
data.list = list(N = length(x),Npreds=ncol(Xmat), Ngroups = npop,y = y, Xmat = Xmat,group = as.numeric(x))

fit2 = stan(file = "anova_random.stan", data = data.list, iter = ni,
            warmup = nb, chains = nc, thin = nt, refresh = 0)