My data are nested as individuals inside families inside of communities. I use three measures of wealth at the family level that can replace each other to create a portfolio measure within my transformed parameters block. I would like to calculate the mean at the community level of the portfolios and use it in the model.

Here is my code in the transformed parameters block without the community level mean.

```
transformed parameters {
vector[f] port; //portfolio values
vector[n] theta; //means in the normal distribution
vector<lower=0>[3] gw; //weightings for wealth measures in portfolio
for (i in 1:3)
gw[i] = exp(phi[i])/(exp(phi[1])+exp(phi[2])+exp(phi[3]));
for (i in 1:f) //for each family
port[i] = gw[1]*inc[i] + gw[2]*acre[i] + gw[3]*tlu[i];
for (i in 1:n)//for each individual
theta[i] = bp*log(port[fam[i]]) + xv[vil[i]] + ba*log(age[i]) + bs*sex[i]; //+ xf[fam[i]]
}
model {
//Priors
bp ~ normal(0,5);
bs ~ normal(0,5);
ba ~ normal(0,5);
gs ~ cauchy(0,2.5);
vs ~ cauchy(0,2.5);
phi[1] ~ normal(0,1);
phi[2] ~ normal(0,1);
phi[3] ~ normal(0,1);
//Likelihood
for (V in 1:v)
xv[V] ~ normal(0,vs);
y ~ normal(theta, gs);
}
```

This code performs exactly as expected. The code below clearly does not work, but it should indicate what I hope to do.

```
transformed parameters {
vector[f] port; //portfolio values
vector[n] theta;
vector<lower=0>[3] gw; //weightings for wealth measures in portfolio
vector[v] vMean; //Hold the village means for v villages
real count;
for (i in 1:3)
gw[i] = exp(phi[i])/(exp(phi[1])+exp(phi[2])+exp(phi[3]));
for (i in 1:f)
port[i] = gw[1]*inc[i] + gw[2]*acre[i] + gw[3]*tlu[i];
for (i in 1:v)//calculate the village means
{
vMean[i] <- mean(port[which(vil==i)] //R code that would do what I want
}
}
for (i in 1:n)
theta[i] = bp*log(port[fam[i]]-vMean[vil[i]]) + xv[vil[i]] + ba*log(age[i]) + bs*sex[i]; //+ xf[fam[i]]
}
model {
//Priors
bp ~ normal(0,5);
bs ~ normal(0,5);
ba ~ normal(0,5);
gs ~ cauchy(0,2.5);
vs ~ cauchy(0,2.5);
phi[1] ~ normal(0,1);
phi[2] ~ normal(0,1);
phi[3] ~ normal(0,1);
//Likelihood
for (V in 1:v)
xv[V] ~ normal(0,vs);
y ~ normal(theta, gs);
}
```

My attempt for the “calculate the village means” section (given my rudimentary knowledge of c++) was

```
for (i in 1:v) //calculate the village means
{
for (j in 1:n)
{
if(vil[j]==i)
{
vMean[i] = vMean[i] + port[fam[j]];
count = count+1;
}
}
vMean[i] = vMean[i]/count;
}
```

I know this is more of a programming question than it is a statistical question, but any assistance you might give would be awesome!