Shrinkage Plots

[dmbates]

• For this and other Makie-based plots I will use the MixedModels v4.0.0-DEV formulation. We can back-port to a v3.6.0 or later, if necessary.
• For shrinkage plots I plan to start with a single random-effects term so that the grouping for each point in the plot is clearly defined.
  • in that case the single term(i.e. only(m.reterms)) contains:
    • the refs vector to determine the groups
    • a levels field to create the labels
    • a cnames field to determine which coefficients are to be shown
• the Xymat field will be split according to refs for fitting the within-group models.
  • Probably drop any columns whose names are not in cnames (names come fromm.feterm)
  • May need to skip levels for which there are too few observations to fit the within-group model
• ranef() returns the conditional modes as a NamedTuple of matrices (levels are columns, for easier display)
  • retrieve fixed-effects estimates and names and add them to the ranef() matrix

[dmbates]

I have committed a shrinkage function to a newly created db/shrinkage branch. This method returns the levels of the grouping factor (currently restricted to a single random-effects term), the mixed-effects coefficients (conditional modes plus fixed effects) for the groups and the within-group estimates.

julia> using DataFrames, MixedModels, MixedModelsMakie
[ Info: Precompiling MixedModelsMakie [b12ae82c-6730-437f-aff9-d2c38332a376]

julia> m = fit(MixedModel, @formula(reaction ~ 1+days+(1+days|subj)), MixedModels.dataset(:sleepstudy))
Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 + days | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -875.9697  1751.9393  1763.9393  1764.4249  1783.0971

Variance components:
            Column    Variance Std.Dev.   Corr.
subj     (Intercept)  565.51068 23.78047
         days          32.68212  5.71683 +0.08
Residual              654.94145 25.59182
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405      6.63226  37.91    <1e-99
days          10.4673     1.50224   6.97    <1e-11
──────────────────────────────────────────────────

julia> DataFrame(shrinkage(m))
18×3 DataFrame
 Row │ levels  within                             mixed                             
     │ String  NamedTup…                          NamedTup…                         
─────┼──────────────────────────────────────────────────────────────────────────────
   1 │ S308    ((Intercept) = 244.193, days = 2…  ((Intercept) = 254.221, days = 1…
   2 │ S309    ((Intercept) = 205.055, days = 2…  ((Intercept) = 211.357, days = 1…
   3 │ S310    ((Intercept) = 203.484, days = 6…  ((Intercept) = 212.972, days = 4…
   4 │ S330    ((Intercept) = 289.685, days = 3…  ((Intercept) = 274.237, days = 5…
   5 │ S331    ((Intercept) = 285.739, days = 5…  ((Intercept) = 272.955, days = 7…
   6 │ S332    ((Intercept) = 264.252, days = 9…  ((Intercept) = 260.221, days = 1…
   7 │ S333    ((Intercept) = 275.019, days = 9…  ((Intercept) = 267.847, days = 1…
   8 │ S334    ((Intercept) = 240.163, days = 1…  ((Intercept) = 244.408, days = 1…
   9 │ S335    ((Intercept) = 263.035, days = -…  ((Intercept) = 250.368, days = -…
  10 │ S337    ((Intercept) = 290.104, days = 1…  ((Intercept) = 286.071, days = 1…
  11 │ S349    ((Intercept) = 215.112, days = 1…  ((Intercept) = 226.847, days = 1…
  12 │ S350    ((Intercept) = 225.835, days = 1…  ((Intercept) = 239.071, days = 1…
  13 │ S351    ((Intercept) = 261.147, days = 6…  ((Intercept) = 255.679, days = 7…
  14 │ S352    ((Intercept) = 276.372, days = 1…  ((Intercept) = 272.027, days = 1…
  15 │ S369    ((Intercept) = 254.968, days = 1…  ((Intercept) = 254.664, days = 1…
  16 │ S370    ((Intercept) = 210.449, days = 1…  ((Intercept) = 226.695, days = 1…
  17 │ S371    ((Intercept) = 253.636, days = 9…  ((Intercept) = 252.128, days = 9…
  18 │ S372    ((Intercept) = 267.045, days = 1…  ((Intercept) = 263.524, days = 1…

I would appreciate if @jkrumbiegel and others could advise me on code style with regard to the plotting.

Instead of producing a column table (in the sense of the Tables package) I could define a struct, say MixedShrinkage, and an AbstractPlotting recipe for that type. For this example the desired plot is a scatter plot on one Axis of the "within" points, the "mixed" points and arrows from the "within" to the "mixed" values for each level. See the examples directory in the db/shrinkage branch.

If there were more than two coefficients for each level then I think the plot should have multiple Axis regions in a scatterplot matrix arrangement (possibly just one triangle of the matrix). I am trying to decide if I should write a recipe for MixedShrinkage that just calls scatter! for each set of coordinates ("within" and "mixed") then calls arrows! and perhaps adds a single point at the global fixed-effects estimate. If the coefficients are other than two-dimensional then a pair of names would be passed as a named argument.

Is the "create a Type and a recipe for said Type" route the general approach to extending CairoMakie?

[jkrumbiegel]

Makie currently only has plot recipes that create nested plot objects. So you build a more complex plot type out of primitives. In your case, that could be scatter + arrow, to form a new plot recipe type called Shrinkage maybe. You would then get plotting functions shrinkage and shrinkage! that work exactly the same as, say, lines or scatter!.

But if you want to have a plot where an axis is created with data-dependent labels, and you maybe want a facet layout with multiple axes, then that's not something currently possible with Makie recipes. We just haven't decided yet, what such an API should look like. So I would simply write functions for that kind of thing now.

You could also split your approach into a plot recipe for the single-axis shrinkage plot (just the scatters connected by arrows, for example), and a function that handles axis / layout creation and uses the shrinkage plot recipe internally.

If I write modular plotting functions these days that create more complex "infrastructure" like Axes, Colorbars, etc., I usually write them so they take a Figure or Figure(Sub)Position as an argument. And because of the nested layouts it will usually "just work", no matter if your object is the top figure, or a nested position inside it. That would enable you to place your whole facet plot in a figure with other content.

Here's a small example, where I write a function that creates a mock facet plot from a matrix.

my_facet_plot(matrix) = my_facet_plot(Figure(), matrix)

function my_facet_plot(f, matrix)
    s = size(matrix)

    axs = map(CartesianIndices(matrix)) do idx
        ax = Axis(f[Tuple(idx)...])
        barplot!([matrix[idx]])
        tightlimits!(ax, Bottom())
        ax
    end

    for col in 1:s[2]
        Box(f[1, col, Top()], color = :gray90)
        Label(f[1, col, Top()], "$col", padding = (5, 5, 5, 5))
    end
        
    for row in 1:s[1]
        Box(f[row, s[2], Right()], color = :gray90)
        Label(f[row, s[2], Right()], "$row", padding = (5, 5, 5, 5), rotation = pi/2)
    end

    linkaxes!(axs...)

    f
end

So now, if I don't pass a figure as an argument, I get a standalone plot (it would be nice to be able to pass figure keywords though, I've left that out here):

my_facet_plot(rand(3, 3))

But I can also pass a Figure, or a FigurePosition, as the first argument, and the code will still work to create a nested version. This way I could easily use this function to create a larger figure:

f = Figure(resolution = (1200, 1200), fontsize = 14)
my_facet_plot(f[1, 1], rand(2, 2))
my_facet_plot(f[1, 2], rand(3, 3))
my_facet_plot(f[2, 1:2], rand(1, 6))
f

[dmbates]

Thanks, @jkrumbiegel I think I have the idea now.

[dmbates]

Just to close this out, #21 provides a shrinkage plot that avoids going through fitting fixed-effects models to each subgroup of observations. When you get into complicated designs with crossed grouping factors it is not easy to decide what should be compared. Instead we use the existing LinearMixedModel struct and insert a very large value of the parameter vector, corresponding to a covariance factor, Λ, that is a large multiple of the identity. This has the effect of removing the shrinkage on the values of b while still allowing them to be evaluated. (If you remove the shrinkage entirely then [Z X] will be rank deficient and the conditional means are undefined.)