Usually, when I do hierarchical modeling, the problem is relatively simple. For instance, let us say I want to know the average weight of frogs in an area, and I collect data on frogs from different ponds. In order to control for non-independence of observations within a pond, one can use a hierarchical model. The hierarchical model 'shrinks' the weight estimate for each pond toward the grand mean, and when computing the grand mean, it takes into account the fact that some ponds may have been under- or over-sampled. This is how one might do it with R for instance:
pond <- c(1,1,1,1,1,2,2,3,3,3,3,3)
weight <- c(23,25,26,34,24,9,10,20,21,22,20,21)
hierarchical_model <- lmer(weight ~ 1 + (1|pond))
This is easy in the sense that the estimate of the parameter 'weight' for each pond is simply an adjusted average of the weight of each frog in the pond. But what if my data are not like that?
For a current project, I am interested in people's beliefs about the distribution of physical strength in the population. I assume that people have an internal mental representation of strength as being normally distributed, and that a person's belief about the distribution of strength can be summarized by a vector (mu, sigma) that is specific to that person.
I have an experiment in which I ask people "If we randomly took someone from the population, what is the probability that this person would lift a stone of X kg?", for various values of X. Then, for each participant, I can find the values of mu and sigma that best describe the participant's answers, using a simple Maximum-Likelihood algorithm.
To find the average mu and the average sigma in my population of participants, I could just average the parameter estimates for each participant, but this feels wrong: some parameter estimates are more precise than others, and I would like to do something like hierarchical modeling in order to take this into account. However, this is not a context where I can use a simple function like lmer. How should I proceed?