Using a group mean of a transformed parameter in a model


#1

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!


#2

Hi, the question is a bit hard to parse as it involves a lot of code that is not really relevant to the question you are asking. You would be more likely to get better answeres if you reduced your code to just a few lines demonstrating the issue you have.

Nevertheless, your proposed solution seems OK.