# Taking the mean of a matrix in transformed parameters block

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?

In Stan is an array. Easiest you could define:

``````matrix [nGroup, nCond] aGxC;
...
real mG = mean(aGxC);
``````

BTW:

If you define `aGxC` as a matrix. One can write `to_vector(aGxC) ~ normal(0, sigma_gc);` and
omit the loops.
The likelihood loop can be vectorized too.

If one wants the keep the array. There are `to_matrix` functions to convert array(s) into a matrix.

Thanks @andre.pfeuffer. How do you vectorize the likelihood?

I used

`score ~ normal(a0 + aGroup[gIndex] + aCond[cIndex] + aSubj[sIndex] + aGxC[gIndex, cIndex], sigma);`

but got

``````SYNTAX ERROR, MESSAGE(S) FROM PARSER:

Unknown variable: to_vector
No matches for:

vector + matrix
``````

and

``````Expression is ill formed
error in 'model287b3c8f7fd3_temp' at line 53, column 96
-------------------------------------------------
51:
52:       // likelihood
53:       score ~ normal(a0 + aGroup[gIndex] + aCond[cIndex] + aSubj[sIndex] + aGxC[gIndex][cIndex], sigma);
^
54:
-------------------------------------------------
``````

So how do I vectorise it?

You could use something like `to_vector(aGxC)[cPrecalculated];`
where `cPrecalculated` is the matrix to vector transform index. Well it doesn’t win a price, but
it should be faster.

Sorry for the follow-up @andre.pfeuffer but what is a “matrix to vector transform index”?

``````> M<-matrix(1:4,2,2)
> as.vector(M)[4]
[1] 4
> M
[,1] [,2]
[1,]    1    3
[2,]    2    4
``````

Matrix position [2,2] transforms to vector position [4].

One can define a calculation put this is `transformed data` section:

int cPrecalculated[N];
for(i in 1:N)
cPrecalculated[i] = gIndex[i] + nCond * (cIndex[i]-1);