I’m modelling data from a (simulated) experiment that has two fully crossed random effects factors: subject and item. In other words, I have N subjects, each of which was tested on M items. Items belong to one of two conditions (although that is not relevant here, I think).
If I model the data as follows (pseudocode):
response ~ grand_mean + grand_slope * condition + subject_means[s] + subject_slopes[s] * condition
and then sample using NUTS (default parameters and initialization), then I get very highly negatively correlated estimates for grand_mean and subject_means, and grand_slope and subject_slopes. This makes sense to me, since of course this model is overspecified. The Bayesian solution to this is partial pooling, where I have
subject_means ~ Normal(mu_a, sd_a)
subject_slopes ~ Normal(mu_b, sd_b)
response ~ subject_means[s] + subject_slopes[s] * condition
This works perfectly fine, removes the very high correlation (except of course between subject_means and mu_a, but those are expected and healthy), and gives a much cleaner posterior (and is also much faster to sample). So far so good.
However, remember that I have two random factors here, the other being item. If I add the item effect (slope effect only) back in:
subject_means ~ Normal(mu_a, sd_a)
subject_slopes ~ Normal(mu_b, sd_b)
response ~ subject_means[s] + subject_slopes[s] * condition + item_slopes[i] * condition
then I end up with the same problem: poor and slow sampling due to extremely high (negative) correlations between item_slopes and subject_slopes. Now I’m at a loss as to how to incorporate the same partial pooling trick for the item factor. I tried the following:
subject_means ~ Normal(mu_a, sd_a)
subject_slopes ~ Normal(mu_b, sd_b)
item_slopes ~ Normal(mu_c, sd_c)
response ~ subject_means[s] + subject_slopes[s] * condition + item_slopes[i] * condition
but that still leaves the model overspecified. Note that the design here is crossed, and not nested: the same items are presented to all the subjects. Therefore, it doesn’t (I think) make sense to have something like item_slopes ~ Normal(mu=subject_slopes)
or so. Especially since each item is only presented exactly once to each subject.
Does anyone have any ideas on how to proceed here? I don’t necessarily need/want to use partial pooling across the items, the only reason I put it in here is to try to avoid the sampling issues (which worked very well when adding it in across the subjects). I guess there might be a way to keep the overspecified model, yet get good samples?