Tyler
2017-Oct-12 17:12 UTC
[Rd] Bug in model.matrix.default for higher-order interaction encoding when specific model terms are missing
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:X3C1 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:X3C1 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:X3C1 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:X3C1 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]]
Arie ten Cate
2017-Oct-15 05:49 UTC
[Rd] 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
Tyler
2017-Oct-15 17:05 UTC
[Rd] Bug in model.matrix.default for higher-order interaction encoding when specific model terms are missing
You could possibly try to explain away the behavior for a missing main effects term, since without the main effects term we don't have main effect columns in the model matrix used to compute the interaction columns (At best this is undocumented behavior--I still think it's a bug, as we know how we would encode the categorical factors if they were in fact present. It's either specified in contrasts.arg or using the default set in options). However, when all the main effects are present, why would the three-factor interaction column not simply be the product of the main effect columns? In my example: we know X1, we know X2, and we know X3. Why does the encoding of X1:X2:X3 depend on whether we specified a two-factor interaction, AND only changes for specific missing interactions? In addition, I can use a two-term example similar to yours to show how this behavior results in a singular covariance matrix when, given the desired factor encoding, it should not be singular. We start with a full factorial design for a two-level continuous factor and a three-level categorical factor, and remove a single row. This design matrix does not leave enough degrees of freedom to determine goodness-of-fit, but should allow us to obtain parameter estimates.> design = expand.grid(X1=c(1,-1),X2=c("A","B","C")) > design = design[-1,] > designX1 X2 2 -1 A 3 1 B 4 -1 B 5 1 C 6 -1 C Here, we first calculate the model matrix for the full model, and then manually remove the X1 column from the model matrix. This gives us the model matrix one would expect if X1 were removed from the model. We then successfully calculate the covariance matrix.> mm = model.matrix(~(X1+X2)^2,data=design) > mm(Intercept) X1 X2B X2C X1:X2B X1:X2C 2 1 -1 0 0 0 0 3 1 1 1 0 1 0 4 1 -1 1 0 -1 0 5 1 1 0 1 0 1 6 1 -1 0 1 0 -1> mm = mm[,-2] > solve(t(mm) %*% mm)(Intercept) X2B X2C X1:X2B X1:X2C (Intercept) 1 -1.0 -1.0 0.0 0.0 X2B -1 1.5 1.0 0.0 0.0 X2C -1 1.0 1.5 0.0 0.0 X1:X2B 0 0.0 0.0 0.5 0.0 X1:X2C 0 0.0 0.0 0.0 0.5 Here, we see the actual behavior for model.matrix. The undesired re-coding of the model matrix interaction term makes the information matrix singular.> mm = model.matrix(~(X1+X2)^2-X1,data=design) > mm(Intercept) X2B X2C X1:X2A X1:X2B X1:X2C 2 1 0 0 -1 0 0 3 1 1 0 0 1 0 4 1 1 0 0 -1 0 5 1 0 1 0 0 1 6 1 0 1 0 0 -1> solve(t(mm) %*% mm)Error in solve.default(t(mm) %*% mm) : system is computationally singular: reciprocal condition number = 5.55112e-18 I still believe this is a bug. Best regards, Tyler Morgan-Wall On Sun, Oct 15, 2017 at 1:49 AM, Arie ten Cate <arietencate at gmail.com> wrote:> 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 >[[alternative HTML version deleted]]
Reasonably Related Threads
- Bug in model.matrix.default for higher-order interaction encoding when specific model terms are missing
- Bug in model.matrix.default for higher-order interaction encoding when specific model terms are missing
- Bug in model.matrix.default for higher-order interaction encoding when specific model terms are missing
- Bug in model.matrix.default for higher-order interaction encoding when specific model terms are missing
- Bug in model.matrix.default for higher-order interaction encoding when specific model terms are missing