In this example, I want to estimate the uncertainty around the beta value for each Cut & Clarity pair. Iโm not sure how to interpret the allpars data. Is it producing a single sigma shared across all values of Cut & Clarity? Do I need to simulate in order to get what I want? I am not trying to get prediction error, I just want the uncertainty around the expected value at each point.
That produces an interval for a single Cut & Clarity parameter, but I was looking for an interval for every value of Cut & Clarity. I want to make a coefficient plot, to show the credible ranges for each value of Cut and Clarity .
Without knowing your inferential and presentation goals and potentially more about your study design, I donโt know whether it makes sense for your purposes. I donโt think a mixed model is the way to analyze the diamonds data, but the presentation you give is a way to get regularized โestimatesโ (technically predictions) of all those things.
Have you looked at MixedModelsMakie? caterpillar and RanefInfo will display similar information using the conditional variances instead of the bootstrap replicates.
I think this is the core thing Iโm getting stuck on โ what do you mean with โbetaโ? Usually, โbetaโ refers to the fixed-effect regression coefficients. The conditional modes / best linear unbiased predictions (BLUPs) / random effects are something else.
Ah yes, beta was not what I meant. I want to show my results in terms of the uncertainty around the expected outcome associated with each value of the categorical parameters. I find this outcome-scale presentation easier to interpret than raw parameter estimates, especially for generalized linear mixed models. It looks like MixedModelsMakie.caterpillar produces estimates in terms of the parameter values, rather than best predictions.
Based on
d[!,:sim] = simulate(modelfit)
let c = combine(groupby(d, [:Cut,:Clarity]), :sim .=> [mean, sem])
layers = (data(c)
* visual(Errorbars)
* mapping(:Cut,:sim_mean,:sim_sem=>(x->2x);layout=:Clarity)
)
draw(layers; axis=(;xticklabelrotation=ฯ/8))
end
Running this code a few times, the predictions jump around quite a bit. E.g.
How should I get reliable estimates/predictions for each level of the categories?
simulate is only a single replication, so things will tend to change a fair amount for any parameter value. For the BLUPs though, they are completely new draws from the random effects distribution. This is one of the critical ways that random effects differ from fixed effects โ youโre estimating the distribution of the random effects, and predicting the individual values. (This is also why the BLUPs arenโt part of the bootstrap summary โ they arenโt parameters and the bootstrap summary only stores parameters.)
Those predictions are relative to the population-level values, i.e. centered at the corresponding fixed effects, so if you want to plot those for a simple intercepts-only model, you can do something like:
A similar result could also be achieved with MixedModelsMakie:
using CairoMakie, MixedModelsMakie
re = ranefinfo(modelfit, grp)
re.ranef[:, 1] .+= fixef(modelfit)[1] # 'absolute' prediction for each group instead of relative to the population effect
caterpillar!(Figure(), re)
Iโm a little bit hesitant to provide any more concrete advice without knowing more about your inferential problem (and I probably donโt have time at the moment to help you with the details of your inferential problem), because this feels like an XY Problem.
