Hi There

I am having trouble taking the mean of cells in a matrix within the transformed parameters block. I need to do this to perform a sum-to-zero procedure (see pp.430, 555, 587 of *Doing Bayesian Data Analysis* by John Kruschke).

Here is the toy data, a mixed design, with both within- and between-subject factors. Each individual gets three conditions, but there are three between-subjects groups as well. There should be differences across condition (cond1 = 100, cond2 = 75, cond3 = 125), but not across groups. So no main effect of group, no main effect of condition, and no interaction. The model doesn’t work super well and takes a while to run (sorry!) but it is closest to my real data and illustrates my problem so I am using it as my example.

The model runs without the transformed parameters block ok, but with it throws an error. In that block I am trying to calculate the mean of all the cells in the matrix of interaction terms `aGxC`

```
# create data
df <- data.frame(id = factor(rep(1:60, 3)),
group = factor(rep(c("g1", "g2", "g3"), each = 20, length.out = 180)),
condition = factor(rep(c("cond1", "cond2", "cond3"), each = 60)),
score = c(ceiling(rnorm(60, 100, 15)), ceiling(rnorm(60, 75, 15)), ceiling(rnorm(60, 125, 15))))
##### Step 1: put data into a list
mixList <- list(N = nrow(df),
nSubj = nlevels(df$id),
nGroup = nlevels(df$group),
nCond = nlevels(df$condition),
nGxC = nlevels(df$group)*nlevels(df$condition),
sIndex = as.integer(df$id),
gIndex = as.integer(df$group),
cIndex = as.integer(df$condition),
score = df$score,
grandMean = 102.0056,
stdDev = 25.81834)
###### Step 2: build model
write("
data{
int<lower=1> N;
int<lower=1> nSubj;
int<lower=1> nGroup;
int<lower=1> nCond;
real grandMean;
real stdDev;
int<lower=1,upper=nSubj> sIndex[N];
int<lower=1,upper=nGroup> gIndex[N];
int<lower=1,upper=nCond> cIndex[N];
real score[N];
}
parameters{
real a0;
vector[nSubj] aSubj;
vector[nGroup] aGroup;
vector[nCond] aCond;
real aGxC [nGroup, nCond];
real<lower=0> sigma;
real<lower=0> sigma_a;
real<lower=0> sigma_s;
real<lower=0> sigma_g;
real<lower=0> sigma_c;
real<lower=0> sigma_gc;
}
transformed parameters{
real mG;
mG = mean(aGxC[1:nGroup,1:nCond]);
}
model{
vector[N] mu;
//hyper-priors
sigma_s ~ normal(0,10);
sigma_g ~ normal(0,10);
sigma_c ~ normal(0,10);
sigma_gc ~ normal(0,10);
//priors
sigma ~ cauchy(0,2);
a0 ~ normal(grandMean, stdDev);
aSubj ~ normal(0, sigma_s);
aGroup ~ normal(0,sigma_g);
aCond ~ normal(0,sigma_c);
for (i in 1:nGroup) {
for (j in 1:nCond) {
aGxC[i,j] ~ normal(0, sigma_gc);
}
}
// likelihood
for (i in 1:N) {
mu[i] = a0 + aGroup[gIndex[i]] + aCond[cIndex[i]] + aSubj[sIndex[i]] + aGxC[gIndex[i], cIndex[i]];
}
score ~ normal(mu, sigma);
}
", file = "temp.stan")
##### Step 3: generate the chains
mix <- stan(file = "temp.stan",
data = mixList,
iter = 4e3,
warmup = 1e3,
cores = 1,
chains = 1)
```

This throws the error

```
No matches for:
mean(real[,])
Available argument signatures for mean:
mean(real[])
mean(vector)
mean(row vector)
mean(matrix)
```

Now I know I am supposed to use a vector or vectors rather than a matrix but I am coming from JAGS where you can take the means of matrices, so am finding it a bit tough to get my head around to do this.

**Can anyone show me a technique for taking the mean of a matrix, or rather how to achieve the same result using a vector of vectors?**