The lme4 package in R has the function
lmList, which I use heavily to examine shrinkage in mixed-models.
MixedModels.jl does not seem to have the function.
MixedModelsMakie.jl is great to visualize shrinkage, but it does not give the me the subgroup data analyses.
Is there a equivalent method in Julia to fit
glm to different subgroups of the data?
There’s no comparable function in Julia. Using the same technique that MixedModelsMakie uses, you could do:
unshrunken_raneftables(m::LinearMixedModel, θref=10^4 .* m.θ; uscale::Bool=false)
Compute the conditional modes at θref.
This function is **not** thread safe because it temporarily mutates
the passed model before restoring its original form.
function unshrunken_raneftables(m::LinearMixedModel, θref=10^4 .* m.θ; uscale::Bool=false)
vv = try
raneftables(updateL!(setθ!(m, θref)); uscale)
@error "Failed to compute unshrunken values with the following exception:"
updateL!(setθ!(m, m.optsum.final)) # restore parameter estimates and update m
model = fit(MixedModel, @formula(reaction ~ 1 + days + (1+days|subj)), MixedModels.dataset(:sleepstudy))
# returns a named tuple of tables, where each entry in the named tuple is a grouping factor
unshrunken = unshrunken_raneftables(model)
# look at the entries for grouping factor
df = DataFrame(unshrunken[:subj])
# add in the fixed effects to get approximate OLS fits --
# not quite OLS because there is still some minimal shrinkage going on,
# but that's necessary to make this work for more complicated models with
# multiple grouping factors or fixed effects without associated random effects
pseudo_ols = transform(df, "(Intercept)" => ByRow(x -> x + fixef(model)),
"days" => ByRow(x -> x + fixef(model)); renamecols=false)
We can also use DataFrame’s convenient
groupby and do real OLS for this simple example:
# now we do the same thing using DataFrames and GLM
sleepstudy = DataFrame(MixedModels.dataset(:sleepstudy))
ols = DataFrame()
for gdf in groupby(sleepstudy, :subj)
group = only(unique(gdf.subj))
m = fit(LinearModel, @formula(reaction ~ 1 + days), gdf)
row = NamedTuple(Symbol.(coefnames(m)) .=> coef(m))
row = merge(row, (; subj=group))
To double check you get similar answers:
all(isapprox.(ols[!, "(Intercept)"], pseudo_ols[!, "(Intercept)"]; atol=0.001))
all(isapprox.(ols[!, "days"], pseudo_ols[!, "days"]; atol=0.001))
Thanks a lot,
unshrunken_raneftables is what I was look for. That’s much more efficient than
ols on grouped data.
for shake of completeness: