Taking the mean of a matrix in transformed parameters block

Modeling
,
1

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?

2

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.

3

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?

4

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.

5

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

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