In some cases when levels of one factor all “inside” levels of another we have not full-rank model matrix from StatsModels. For example - calssical bioequivlence crossover study (2X2) with factors: Var ~ Formulation + Sequence + Period + Subject - levels of Subject distributed in Sequnce. So we get result bellow.
Should we get operator for @formula to make coding inside each level of top level factor?
For example, we have data and code bellow. Result will be simmilar to posted above. Subject: 1 not included in model matrix (this is a base level), but Subject: 15 included and it make collinearity. So, if make DummyCoding separetly for subjects in TR sequence and in RT sequence there will be no problems.
Subject Sequence Period Formulation Var
3 TR 1 T 225.95
1 RT 1 R 181.09
2 RT 1 R 114.48
4 RT 1 R 176.91
5 TR 1 T 147.01
6 TR 1 T 97.53
7 RT 1 R 146.60
8 TR 1 T 45.58
9 RT 1 R 109.20
10 RT 1 R 125.61
11 TR 1 T 92.26
12 RT 1 R 237.95
13 TR 1 T 145.46
14 TR 1 T 179.96
15 TR 1 T 173.86
16 RT 1 R 144.00
17 RT 1 R 185.10
18 TR 1 T 117.99
1 RT 2 T 210.14
2 RT 2 T 98.72
3 TR 2 R 241.09
4 RT 2 T 186.65
5 TR 2 R 139.56
6 TR 2 R 124.77
7 RT 2 T 137.62
8 TR 2 R 57.71
9 RT 2 T 139.36
10 RT 2 T 120.43
11 TR 2 R 116.10
12 RT 2 T 228.63
13 TR 2 R 165.09
14 TR 2 R 181.09
15 TR 2 R 206.66
16 RT 2 T 143.25
17 RT 2 T 192.22
18 TR 2 R 125.50
and I get ok (as you will see, there are some warnings because some of your categories have no observation, which is probably because you have extracted only a part of the dataset). I also added the VIF to have something about the potential collinearity. I haven’t looked in much detail, but the lack of observations for some categories probably gives weird results in GLM. I would recommend using FreqTables.jl from @nalimilan to get a more precision where the problem lies.
Thought if you want to be specific about which category is the base you can specify it with the contrast (same format as in GLM). you can also find an example here: regression example with contrasts
Hi! Thanks a lot for LinearRegressionKit.jl, it works perfect!
But this problem is not because data is incomplete. It is complete. This happens because we have in same time Sequence factor and Subject, and Subject is nested with Sequence - One level of Subject corresponds to only one level of Sequence. In some software (IBM SPSS for example) for this case used syntax like this: /DESIGN=Sequence Period Formulation Subject(Sequence) - this means that even Subject have same levels they coded separately for each Sequence, and also this means that DummyCoding (for n-1 factor levels) used for each Sequence (so if we have 2 sequences - we have n - 2 subject levels, if 3 sequences we have n - 3 subject levels ets.).
Using DataFrames, CSV, GLM
model = lm(@formula(log(Var) ~ Sequence+Period+Formulation+Subject), dfrds; dropcollinear = true,)
using MixedAnova
glm_init()
MixedAnova.anova(model, type = 3)
ANOVA results:
Analysis of Variance
Type 3 test / F test
:(log(Var)) ~ 1 + Sequence + Period + Formulation + Subject
Table:
───────────────────────────────────────────────────────────
DOF Exp.SS Mean Square F value Pr(>|F|)
───────────────────────────────────────────────────────────
(Intercept) 1 49.86 49.86 7795.7027 <1e-22
Sequence 1 1.7847 1.7847 279.0437 <1e-10
Period 1 0.0453 0.0453 7.0905 0.0170
Formulation 1 0.0228 0.0228 3.5725 0.0770
Subject 17
(Residuals) 16 0.1023 0.0064
───────────────────────────────────────────────────────────
We see: wrong DOF for Subject, absence SS and MS for Subject, wrong F value for Sequence (because Sequence is inter-subject factor and F should be = MSsequence / MSsubject), not residual. It is not problem of coding but it can help packages like MixedAnova.jl correctly handle this data structure.
And the question: do we need in StatsModels some operator like var ~ factor1 + factor2 in factor1 to handle this models correctly?
Thanks for providing more details.
If I understand correctly, you’re looking for nested effects.
I believe they are defined in MixedModels.jl; there is a documentation here I believe they use the / as it builds on the syntax from R.
Yep! I use it when possible MixedModels, Metida, GLM, LsqFit, FixedEffectsModels, ets… But sometime by statistical plan I should use ANOVA. And there is no problem for me - for some task I can use R or SPSS or write self-implemented anova. But even in self implemented tasks I found that I can’t make full-rank model matrix with StatsModels for some common cases. So maybe additional functionality for StatsModels help developers to make appropriate implementation of ANOVA…
I think part of the idea is to automatically drop columns corresponding to unestimable contrasts, e.g. when certain combinations of levels from different categorical predictors never co-occur – missing or empty cells in a quasi-factorial design. Currently, MixedModels.jl and GLM.jl both handle this case as part of more general pivoting / rank deficiency detection, but it would be good to support dropping clearly unestimable columns (e.g., constant and constant multiples of other columns) from the model matrix or even prevent the allocation of such columns. I’ve talked to @dave.f.kleinschmidt a bit about this and part of the problem at the moment IIRC is that this falls in awkwardly between apply_schema (schema-time) and modelcols (data-time) operations.
Yes! Thank you. That’s exactly what i wanted to say. Second part of suggestion - decide - do we need operator in StatModels or in RegressionFormulae that make model constructing more clear for nested factors. I look at RegressionFormulae it gives ^ and / - and it just expand term to others. Idea is make some operator for example ∈ with following mechanism:
To unpack this a bit more: currently, we have three stages between @formula(...) and actual matrices that get passed to GLM.jl etc.:
The formula expr is parsed into Terms (purely syntactic operations, nothing known about types of variables)
A schema is extracted from the data: categorical (i.e. non-numeric) variables are detected adn their unique values are extracted, one column at a time. This schema is then applied to the parsed Terms, replacing Term(:a) with e.g. CategoricalTerm(:a, ...).
These “concrete terms” can now be fed a table of data and will generate the numeric matrices
The issue is that, in order to support this kind of operator (create different contrasts given some kind of special syntax), we’d need to know about that syntax BEFORE the schema is computed; when the schema is applied, the full dataset is not available (by design). But special syntax is (at the moment) supported by operations that take place during schema application.
In the short term, I’m afraid the best solution is to create a new variable where you just string the pairs together; then the contrasts shoudl be computed correctly.
I try to intercept in function term to make specific term with this function for better nested contrast :
function StatsModels.apply_schema(
t::FunctionTerm{typeof(∈)},
sch::Schema,
Mod::Type{<:RegressionModel},
)
length(t.args_parsed) == 2 ||
throw(ArgumentError("malformed term: $t (Exactly two arguments required)"))
#=
How to get data here?
=#
#Make specific term NestedTerm(sym, contrasts, base)
end
But how to get access to data here? There are too many apply_schema in StstsModels…
Or should it be done in schema stage? Or should I reload concrete_term?
StstsModels is a very confusing package - so I can’t to trace the way of FunctionTerm{typeof(∈)}
You’re lookin’ at it (in the issue ). It’s a pretty disruptive change so it needs some serious time to sit down and get it right and there’s a number of higher priority things. Plus I write julia all day long at my day job now so there’s not as much bandwidth to push this stuff forward as much as I’d like!
RE “how to get data here”: you can’t, and that’s by design. apply_schema does just that: applies teh information in the schema to transform terms into different terms. If the information you need isn’t present in the schema, you can’t do what you’re trying to do I’m afraid.
The reason that constraint is important (IMO) is that we want to support streaming/online/out of core stuff too, or at least not prevent support for it in the future. So lots of stuff that seems reasonable (like having access during apply_schema to the full dataset, which are usually pretty small data frames in teh grand scheme of things) are not allowed.
If we had a separate “apply context” stage, you could hook in there to create your NestedTerm, and then we could specialize schema extraction on these custom term types. I do think that’s a better design overall but like I said above, don’t really have the bandwidth to really think through how to make it work without adding TOO much extra complexity (which, as you already pointed out, already makes StatsModels.jl an intimidating package for many developers)
Second part of suggestion - decide - do we need operator in StatModels or in 
). It’s a pretty disruptive change so it needs some serious time to sit down and get it right and there’s a number of higher priority things. Plus I write julia all day long at my day job now so there’s not as much bandwidth to push this stuff forward as much as I’d like!