Bug in model.matrix.default for higher-order interaction encoding when specific model terms are missing
I think it is not a bug. It is a general property of interactions.
This property is best observed if all variables are factors
(qualitative).
For example, you have three variables (factors). You ask for as many
interactions as possible, except an interaction term between two
particular variables. When this interaction is not a constant, it is
different for different values of the remaining variable. More
precisely: for all values of that variable. In other words: you have a
three-way interaction, with all values of that variable.
An even smaller example is the following script with only two
variables, each being a factor:
df <- expand.grid(X1=c("p","q"), X2=c("A","B","C"))
print(model.matrix(~(X1+X2)^2 ,data=df))
print(model.matrix(~(X1+X2)^2 -X1,data=df))
print(model.matrix(~(X1+X2)^2 -X2,data=df))
The result is:
(Intercept) X1q X2B X2C X1q:X2B X1q:X2C
1 1 0 0 0 0 0
2 1 1 0 0 0 0
3 1 0 1 0 0 0
4 1 1 1 0 1 0
5 1 0 0 1 0 0
6 1 1 0 1 0 1
(Intercept) X2B X2C X1q:X2A X1q:X2B X1q:X2C
1 1 0 0 0 0 0
2 1 0 0 1 0 0
3 1 1 0 0 0 0
4 1 1 0 0 1 0
5 1 0 1 0 0 0
6 1 0 1 0 0 1
(Intercept) X1q X1p:X2B X1q:X2B X1p:X2C X1q:X2C
1 1 0 0 0 0 0
2 1 1 0 0 0 0
3 1 0 1 0 0 0
4 1 1 0 1 0 0
5 1 0 0 0 1 0
6 1 1 0 0 0 1
Thus, in the second result, we have no main effect of X1. Instead, the
effect of X1 depends on the value of X2; either A or B or C. In fact,
this is a two-way interaction, including all three values of X2. In
the third result, we have no main effect of X2, The effect of X2
depends on the value of X1; either p or q.
A complicating element with your example seems to be that your X1 and
X2 are not factors.
Arie
On Thu, Oct 12, 2017 at 7:12 PM, Tyler <tylermw at gmail.com> wrote:
Hi, I recently ran into an inconsistency in the way model.matrix.default handles factor encoding for higher level interactions with categorical variables when the full hierarchy of effects is not present. Depending on which lower level interactions are specified, the factor encoding changes for a higher level interaction. Consider the following minimal reproducible example: --------------
runmatrix = expand.grid(X1=c(1,-1),X2=c(1,-1),X3=c("A","B","C"))> model.matrix(~(X1+X2+X3)^3,data=runmatrix) (Intercept) X1 X2 X3B X3C X1:X2 X1:X3B X1:X3C X2:X3B X2:X3C X1:X2:X3B X1:X2:X3C
1 1 1 1 0 0 1 0 0 0 0 0 0 2 1 -1 1 0 0 -1 0 0 0 0 0 0 3 1 1 -1 0 0 -1 0 0 0 0 0 0 4 1 -1 -1 0 0 1 0 0 0 0 0 0 5 1 1 1 1 0 1 1 0 1 0 1 0 6 1 -1 1 1 0 -1 -1 0 1 0 -1 0 7 1 1 -1 1 0 -1 1 0 -1 0 -1 0 8 1 -1 -1 1 0 1 -1 0 -1 0 1 0 9 1 1 1 0 1 1 0 1 0 1 0 1 10 1 -1 1 0 1 -1 0 -1 0 1 0 -1 11 1 1 -1 0 1 -1 0 1 0 -1 0 -1 12 1 -1 -1 0 1 1 0 -1 0 -1 0 1 attr(,"assign") [1] 0 1 2 3 3 4 5 5 6 6 7 7 attr(,"contrasts") attr(,"contrasts")$X3 [1] "contr.treatment" -------------- Specifying the full hierarchy gives us what we expect: the interaction columns are simply calculated from the product of the main effect columns. The interaction term X1:X2:X3 gives us two columns in the model matrix, X1:X2:X3B and X1:X2:X3C, matching the products of the main effects. If we remove either the X2:X3 interaction or the X1:X3 interaction, we get what we would expect for the X1:X2:X3 interaction, but when we remove the X1:X2 interaction the encoding for X1:X2:X3 changes completely: --------------
model.matrix(~(X1+X2+X3)^3-X1:X3,data=runmatrix) (Intercept) X1 X2 X3B X3C X1:X2 X2:X3B X2:X3C X1:X2:X3B X1:X2:X3C
1 1 1 1 0 0 1 0 0 0 0 2 1 -1 1 0 0 -1 0 0 0 0 3 1 1 -1 0 0 -1 0 0 0 0 4 1 -1 -1 0 0 1 0 0 0 0 5 1 1 1 1 0 1 1 0 1 0 6 1 -1 1 1 0 -1 1 0 -1 0 7 1 1 -1 1 0 -1 -1 0 -1 0 8 1 -1 -1 1 0 1 -1 0 1 0 9 1 1 1 0 1 1 0 1 0 1 10 1 -1 1 0 1 -1 0 1 0 -1 11 1 1 -1 0 1 -1 0 -1 0 -1 12 1 -1 -1 0 1 1 0 -1 0 1 attr(,"assign") [1] 0 1 2 3 3 4 5 5 6 6 attr(,"contrasts") attr(,"contrasts")$X3 [1] "contr.treatment"
model.matrix(~(X1+X2+X3)^3-X2:X3,data=runmatrix) (Intercept) X1 X2 X3B X3C X1:X2 X1:X3B X1:X3C X1:X2:X3B X1:X2:X3C
1 1 1 1 0 0 1 0 0 0 0 2 1 -1 1 0 0 -1 0 0 0 0 3 1 1 -1 0 0 -1 0 0 0 0 4 1 -1 -1 0 0 1 0 0 0 0 5 1 1 1 1 0 1 1 0 1 0 6 1 -1 1 1 0 -1 -1 0 -1 0 7 1 1 -1 1 0 -1 1 0 -1 0 8 1 -1 -1 1 0 1 -1 0 1 0 9 1 1 1 0 1 1 0 1 0 1 10 1 -1 1 0 1 -1 0 -1 0 -1 11 1 1 -1 0 1 -1 0 1 0 -1 12 1 -1 -1 0 1 1 0 -1 0 1 attr(,"assign") [1] 0 1 2 3 3 4 5 5 6 6 attr(,"contrasts") attr(,"contrasts")$X3 [1] "contr.treatment"
model.matrix(~(X1+X2+X3)^3-X1:X2,data=runmatrix) (Intercept) X1 X2 X3B X3C X1:X3B X1:X3C X2:X3B X2:X3C X1:X2:X3A X1:X2:X3B X1:X2:X3C
1 1 1 1 0 0 0 0 0 0 1
0 0
2 1 -1 1 0 0 0 0 0 0 -1
0 0
3 1 1 -1 0 0 0 0 0 0 -1
0 0
4 1 -1 -1 0 0 0 0 0 0 1
0 0
5 1 1 1 1 0 1 0 1 0 0
1 0
6 1 -1 1 1 0 -1 0 1 0 0
-1 0
7 1 1 -1 1 0 1 0 -1 0 0
-1 0
8 1 -1 -1 1 0 -1 0 -1 0 0
1 0
9 1 1 1 0 1 0 1 0 1 0
0 1
10 1 -1 1 0 1 0 -1 0 1 0
0 -1
11 1 1 -1 0 1 0 1 0 -1 0
0 -1
12 1 -1 -1 0 1 0 -1 0 -1 0
0 1
attr(,"assign")
[1] 0 1 2 3 3 4 4 5 5 6 6 6
attr(,"contrasts")
attr(,"contrasts")$X3
[1] "contr.treatment"
--------------
Here, we now see the encoding for the interaction X1:X2:X3 is now the
interaction of X1 and X2 with a new encoding for X3 where each factor level
is represented by its own column. I would expect, given the two column
dummy variable encoding for X3, that the X1:X2:X3 column would also be two
columns regardless of what two-factor interactions we also specified, but
in this case it switches to three. If other two factor interactions are
missing in addition to X1:X2, this issue still occurs. This also happens
regardless of the contrast specified in contrasts.arg for X3. I don't see
any reasoning for this behavior given in the documentation, so I suspect it
is a bug.
Best regards,
Tyler Morgan-Wall
[[alternative HTML version deleted]]
______________________________________________ R-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel