model{
  ## Model Estimates
  for(m in 1:M_count){ # method loop
    for(p in 1:P_count) { # province loop matched to C
      for (t in 1:n_years) { 
        z[m,p,t] <- alpha_pms[m,p] + inprod(Bik[p,t,],beta.k[m,p,])
      } # end time loop
    } # end M loop 
  } # end P loop


## Parameter Estimates 
  for(p in 1:P_count) { # province loop
    alpha_pms[1:M_count,p] <- alpha_cms[1:M_count,matchcountry[p]] + err_alpha_pms[1:M_count,p] # NCP parametrization
    for(m in 1:M_count){ # method loop
      beta.k[m,p,kstar[p]] <- 0
      for(j in 1:(kstar[p]-1)){ # before kstar[p] 
        beta.k[m,p,(kstar[p] - j)] <- beta.k[m,p,(kstar[p] - j)+1] - delta.k[m,p,(kstar[p] - j)]
      } # end K1 loop
      for(j in (kstar[p]+1):(K-1)){ # after kstar[p] 
        beta.k[m,p,j] <- beta.k[m,p,(j-1)] + delta.k[m,p,(j-1)]
      } # end K2 loop
      beta.k[m,p,K] <- -(alpha_pms[m,p] + sum(beta.k[m,p,1:(K-1)])) # hard sum-to-0 constraint on spline coefficients
      for(j in 1:H){
        delta.k[m,p,j] ~ dnorm(0, tau_delta)
      } # end H loop
    } # end m loop
  } # end P loop

## Estimating all the Categories here (including total private)
  for(p in 1:P_count){ # province loop
    for(m in 1:M_count){ # method loop
      for (t in 1:n_years) { # years loop  
        P[1,m,p,t] <- 1/(1+exp(-(z[m,p,t]))) # modelling this as before assuming that z[m,p,t] is log(pi_public/pi_privat)
        P[2,m,p,t] <- 1-P[1,m,p,t] # this then gives you the total private sector
      }
    }
  }

## Likelihood
  for (k in 1:n_obs) {
    y[k] ~ dnorm(z[matchmethod[k],matchsubnat[k],matchyears[k]],  tau_y[k]) # Logit Normal data model
    y.sim[k] ~ dnorm(z[matchmethod[k],matchsubnat[k],matchyears[k]],  tau_y[k]) 
    tau_y[k] <- 1/(se_prop[k]^2) 
  }
    
## Priors
  # Inverse-Wishart prior for country-level params
  inv.Sigma.alpha_cms ~ dwish(Omega, M_count + 1)

  # Inverse-Wishart prior for province-level params
  inv.Sigma.alpha_pms ~ dwish(Omega, M_count + 2)
  
  # Rates of change
  sigma_delta ~ dnorm(0, 2^-2)T(0,) 
  tau_delta <- sigma_delta^-2
  
  # Intercepts
  for(m in 1:M_count){ # method loop
    mu_alpha_cms[m] ~ dnorm(0, pow(10, -2))
  }
  for(c in 1:C_count) { # country loop
    err_alpha_cms[1:M_count,c] ~ dmnorm(rep(0,M_count), inv.Sigma.alpha_cms[1:M_count,1:M_count])
    alpha_cms[1:M_count,c] <- mu_alpha_cms + err_alpha_cms[1:M_count,c]
  }
  for(p in 1:P_count) { # province loop
    err_alpha_pms[1:M_count,p] ~ dmnorm(rep(0,M_count), inv.Sigma.alpha_pms[1:M_count,1:M_count])
  }
}
