Question about array multiplication in JAGS

Question

I am working with race-stratified population estimates and I want to integrate race-stratified populations from three different data sources (census, PEP, and ACS). I developed a model to use information from all these three sources and estimate the true population which is defined as gamma.ctr for county c time t and race r (1=white and 2 for non-white).

The problem is that PEP data is not race-stratified and I need to find a way to estimate race-stratified pep data.

Before, I used one of the other two sources (census or ACS) to estimate ethnicity proportions and multiply PEP data by these proportions to obtain race-stratified PEP population as input data to the model.

Now I decided to do this multiplication within the model based on ethnicity proportions that are defined by gamma.ctr (true pop in county c, year t, and race r) which is updated by all data sources not one of them.

So I considered the input PEP data as peppop.ct (the population for county c and time t, not race-stratified). Then I defined ethnicity proportion as:

prob[c,t]=gamma.ctr[c,t,1]/(gamma.ctr[c,t,1]+gamma.ctr[c,t,2])

I want to multiply PEP data by these proportions to find race-stratified estimates within the JAGS model:

for (c in 1:Narea){
      for (t in 1:nyears){
      prob.ct[c,t]<-gamma.ctr[c,t,1]/(gamma.ctr[c,t,1]+gamma.ctr[c,1,2])
      peppop.ctr[c,t,1]<-peppop.ct[c,t] * prob.ct[c,t]
      peppop.ctr[c,t,2]<-peppop.ct[c,t] * (1-prob.ct[c,t])
      }
}

I want to use this peppop.ctr as a response varaible later like this:

for (t in 1:nyears){
      peppop.ctr[c,t,r] ~ dnorm(gamma.ctr[c,t,r], taupep.ctr[c,t,r])
      }

But I receive this error: Attempt to redefine node peppop.cpr[1,1,1]

It think the reason for this error is the fact that peppop.ctr are defined twice in left hand side of the equation and the error is related to redefining peppop.ctr in line:

peppop.ctr[c,t,1]<-peppop.ct[c,t] * prob.ct[c,t]

Is it possible to help me to solve this error. I need to estimate peppop.ctr first and then use these estimates to update gamma.ctr parameters. Any help is really appreciated.

Answer 1

You can use the zeros trick to both define a variable (e.g., y below) and then also use that variable as the dependent variable in some subsequent analysis. Here's an example:

library(runjags)


x <- rnorm(1000)
y <- 2 + 3 * x + rnorm(1000)
p <- runif(1000, .1, .9)
w <- y*p
z <- y-w

datl <- list(
  x=x, 
  w=w, 
  z=z, 
  zeros = rep(0, length(x)), 
  N = length(x)
)

  
mod <- "model{
  y <- w + z 
  C <- 10000 # this just has to be large enough to ensure all phi[i]'s > 0
  for (i in 1:N) {
    L[i] <- dnorm(y[i], mu[i], tau)
    mu[i] <- b[1] + b[2]*x[i]
    phi[i] <- -log(L[i]) + C
    zeros[i] ~ dpois(phi[i])
    
  }
  #sig ~ dunif(0, sd(y))
  #tau <- pow(sig, -2)
  tau ~ dgamma(1,.1)
  b[1] ~ dnorm(0, .0001)
  b[2] ~ dnorm(0, .0001)
}
"

out <- run.jags(model = mod, data=datl, monitor = c("b", "tau"), n.chains = 2)

summary(out)
#>        Lower95   Median  Upper95     Mean         SD Mode        MCerr MC%ofSD
#> b[1] 1.9334538 1.991586 2.051245 1.991802 0.03026722   NA 0.0002722566     0.9
#> b[2] 2.9019547 2.963257 3.023057 2.963190 0.03057184   NA 0.0002744883     0.9
#> tau  0.9939587 1.087178 1.183521 1.087744 0.04845667   NA 0.0004280217     0.9
#>      SSeff        AC.10      psrf
#> b[1] 12359 -0.010240572 1.0000684
#> b[2] 12405 -0.006480322 0.9999677
#> tau  12817  0.010135609 1.0000195
summary(lm(y ~ x))
#> 
#> Call:
#> lm(formula = y ~ x)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.2650 -0.6213 -0.0032  0.6528  3.3956 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  1.99165    0.03034   65.65   <2e-16 ***
#> x            2.96340    0.03013   98.34   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.9593 on 998 degrees of freedom
#> Multiple R-squared:  0.9065, Adjusted R-squared:  0.9064 
#> F-statistic:  9671 on 1 and 998 DF,  p-value: < 2.2e-16

^{Created on 2022-05-14 by the reprex package (v2.0.1)}