Writing a hierarchical model - unsure where to assign population-level parameters

[josephlewis]

Hi @florianhartig,

I'm trying to write a hierarchical model where I can estimate both population-level and individual-level parameters but I'm not sure how to do this. Namely, where do I assign the population-level prior? Before the likelihood function and then call the global variable in the likelihood function or within the function? Since this population-level prior will influence the individual-level prior I'm having difficulty with how to have the individual-level prior also be estimated and returned. Apologies for the lack of a reproducible example. Hopefully the provided information is sufficient. Thank you.

For example, this is my hierarchical model attempt. Here I can estimate the b_rate but I'm not sure how to also estimate and return the individual b_exp.

# this is the population-level prior for sd
b_rate <- param[1]
b_exp <- rexp(n = n, rate = b_rate)

# I also want to return the value for each individual, i.e. b1 and b2. so I know how each individual varies whilst also knowing the population-level parameter values so I can predict new groups.
b1 <- truncnorm::rtruncnorm(n = n, a = 0, mean = 3.5, sd = b_exp)
b2 <- truncnorm::rtruncnorm(n = n, a = 0, mean = 3.5, sd = b_exp)

Below are the mocked-up complete pooling and no pooling approaches if these help.

Complete pooling approach example:

density_1 = function(par){
  d1 = dlnorm(par[1], meanlog = log(3.5), sdlog = log(2), log = TRUE)
  d2 = dbeta(x = par[2], shape1 = 5, shape2 = 50, log = TRUE)
  return(d1+d2)
}

sampler_1 = function(n = 1){
  d1 = rlnorm(n, meanlog = log(3.5), sdlog = log(2))
  d2 = rbeta(n, shape1 = 5, shape2 = 50)
  return(cbind(d1, d2))
}

likelihood_1 <- function(param) { 
  
  b = param[1]
  c = param[2]
  
  ## process model and calculate log-likelihood
  # The same b and c parameters are used for modelling all the dataset. The process_model() function loops over individual subsets an sums the log-likelihoods.
  ll <- process_model(b,c)
  
  return(ll)
}

The no-pooling approach example:

density_2 = function(par){
  d1 = dlnorm(par[1], meanlog = log(3.5), sdlog = log(2), log = TRUE)
  d2 = dbeta(x = par[2], shape1 = 5, shape2 = 50, log = TRUE)
  d3 = dlnorm(par[3], meanlog = log(3.5), sdlog = log(2), log = TRUE)
  d4 = dbeta(x = par[4], shape1 = 5, shape2 = 50, log = TRUE)
  return(d1+d2+d3+d4)
}

sampler_2 = function(n = 1){
  d1 = rlnorm(n, meanlog = log(3.5), sdlog = log(2))
  d2 = rbeta(n, shape1 = 5, shape2 = 50)
  d3 = rlnorm(n, meanlog = log(3.5), sdlog = log(2))
  d4 = rbeta(n, shape1 = 5, shape2 = 50)
  return(cbind(d1, d2, d3, d4))
}

likelihood_2 <- function(param) { 
  
  b1 = param[1]
  c1 = param[2]
  
  b2 = param[1]
  c2 = param[2]

  ## process model and calculate log-likelihood
  # Different b and c parameters are used for modelling each 'group' of the dataset. The log-likelihoods are summed at the end. The process_model here differs from the complete pooling approach by only looping over one group of the dataset rather than all the groups.
  ll_1 <- process_model(b1,c1)
  ll_2 <- process_model(b2,c2)
  
  return(ll+ll2)
}

[FelixNoessler]

Hey,

if you are still interest how to write the hierarchical model (parameter with partial pooling)... I will show it for the parameter b. In the same way you can do it for other parameters (even variance parameters). You have to specify a prior for the poulation $b_{population}$ and the adjustments from the population to the individuals $b_{ind,adj}$. This adjustment prior is usually a normal distribution with zero mean, for the population prior you can choose any prior. For example:

$$ \begin{align} b_{population} &\sim N(50, 5) \\ b_{ind,adj} &\sim N(0, 2) \\ b_{ind,i} &= b_{population} + b_{ind,adj,i} \end{align} $$

Then we can define the prior and assume five individuals:

log_prior <- function(par){
    b_population = par[1]
    b_ind_adj = par[2:6]

    # if your parameter for the individual has to be positive
    if (any((b_population + b_ind_adj) <= 0)) {return(-Inf)}

    logprior_population <- dnorm(b_population, mean = 50, sd = 5, log = T)
    logprior_individual_adjustment <- sum(dnorm(b_ind_adj, mean = 0, sd = 2, log = T))
   
    return(logprior_population + logprior_individual_adjustment)
}

and define the likelihood:

log_likelihood <- function(par){
    b_population = par[1]
    b_ind_adj = par[2:6]
    b_ind = b_population + b_ind_adj
    
    ll = 0
    for (i in 1:5) {
        ll = ll + process_model(b = b_ind[i])
    }

   return(ll)
}

If you are interested to get the parameter b for all individuals from the final MCMC chains, you have to add to the population parameter $b_{population}$ the adjustment for each individual $b_{ind,adj}$.

You can also check the example random effect model here in the Bayesian Tools package for a statistical hierarchical model.