Oxide example model fit may be an inappropriate model

[dmbates]

One of the tests in test/pls.jl uses the :oxide models defined in test/modelcache.jl

    :oxide => [@formula(Thickness ~ 1 + (1|Lot/Wafer)),
               @formula(Thickness ~ 1 + Source + (1+Source|Lot) + (1+Source|Lot&Wafer))],

The second model is notoriously hard to fit. Different optimizers give very different parameter estimates but with similar values of the objective. I think this is because the model is ill-defined as Source is constant within each Lot.

julia> groupby(DataFrame(MixedModels.dataset(:oxide)), [:Source, :Lot])
GroupedDataFrame with 8 groups based on keys: Source, Lot
First Group (9 rows): Source = "1", Lot = "1"
 Row │ Source  Lot     Wafer   Site    Thickness 
     │ String  String  String  String  Float64   
─────┼───────────────────────────────────────────
   1 │ 1       1       1       1          2006.0
   2 │ 1       1       1       2          1999.0
   3 │ 1       1       1       3          2007.0
   4 │ 1       1       2       1          1980.0
   5 │ 1       1       2       2          1988.0
   6 │ 1       1       2       3          1982.0
   7 │ 1       1       3       1          2000.0
   8 │ 1       1       3       2          1998.0
   9 │ 1       1       3       3          2007.0
⋮
Last Group (9 rows): Source = "2", Lot = "8"
 Row │ Source  Lot     Wafer   Site    Thickness 
     │ String  String  String  String  Float64   
─────┼───────────────────────────────────────────
   1 │ 2       8       1       1          1996.0
   2 │ 2       8       1       2          1989.0
   3 │ 2       8       1       3          1996.0
   4 │ 2       8       2       1          1997.0
   5 │ 2       8       2       2          1993.0
   6 │ 2       8       2       3          1996.0
   7 │ 2       8       3       1          1990.0
   8 │ 2       8       3       2          1989.0
   9 │ 2       8       3       3          1992.0

To me this means that you can't expect to fit a random-effects term like (1 + Source|Lot).

Am I confusing myself?

[dmbates]

If you use a contrast like EffectsCoding() for Source then the conditional means of the random effects for the term (1 + Source | Lot & Wafer) end up being close to multiples of [1, -1] for the first level of Source and [1, 1] for the second level.

julia> first(m2.b)
2×24 Matrix{Float64}:
  3.70686  -5.46891   2.67088  -0.979013  -1.71899  -2.75497  …  -0.997343  -0.39882   -0.0995575  -0.178385  0.56977   -1.67469
 -3.72693   5.49852  -2.68534   0.984313   1.7283    2.76988     -1.00222   -0.400768  -0.100044   -0.179257  0.572553  -1.68288

I'm not sure if that is a consequence of an unstable model or of only having 3 observations for each Lot & Wafer combination.