# ANOVA with random effects in model matrix format

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 + 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)``````

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?

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)``````