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)
}
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.
Below are the mocked-up complete pooling and no pooling approaches if these help.
Complete pooling approach example:
The no-pooling approach example: