R is creating dummy variables for the color variable (color is either white or red) in my linear regression model. For example, color:pH returns the interaction term colorred:pH. Some questions:

Why is R only returning colorwhite for some interaction terms and colorred with other terms? Why doesn't it do colorwhite and colorred for each term interacting with color?

Next, if I want to rebuild a model using only some of the dummy variables R generated, how can I reference them in a linear model? If I retype the model manually and include colorwhite:pH as a term, R returns an error that the object colorwhite isn't found.

enter image description here

  • 1
    $\begingroup$ Can you supply the formula object you wrote to fit this model? I think you should be using the * notation to generate interaction effects instead of :. $\endgroup$ – AdamO Feb 24 '14 at 20:54
  • $\begingroup$ I did a stepwise regression and passed in the model (quality~.,data=foo) which stepped through all of my variable/model combinations and generated these interaction terms on its own. $\endgroup$ – Info5ek Feb 24 '14 at 22:15
  • $\begingroup$ You still had to supply a formula to the scope or upper/lower commands to stepwise. We need to see that to understand exactly which terms are being added/dropped to this model. $\endgroup$ – AdamO Feb 24 '14 at 22:36
  • $\begingroup$ This model is strange in that it appears selectively to drop parts of categorical variables from the interactions. That must have happened in its specification; the software would issue a warning message if it dropped levels due to collinearity. Normally, such specifications would be considered to be erroneous: when you include a categorical variable in an interaction than (a) you should be including all its levels in that interaction and (b) you should (usually) be including that variable--consisting of all its levels--by itself, too. $\endgroup$ – whuber Feb 24 '14 at 23:04
  • 1
    $\begingroup$ I think there was some sort of error in my input, as after retrying this stepwise test from scratch, these colorred/colorwhite terms are gone and only color remains. $\endgroup$ – Info5ek Feb 25 '14 at 0:37

If you use the proper multiplicative notation, the model coefficients need to account for the interpretation of the intercept term. Assume WLOG that only color and ph are in the model (this is a superfluous example you've provided). In this case, if "color==red" is the default group. There's technically only 1 dummy in the model, 1 if color is white, 0 otherwise.

Then, fitting the pH interaction, the "colorwhite" parameter is interpreted as the expected difference in the outcome comparing white to red having a pH of exactly 0. There is also a pH parameter interpreted as an expected difference in the outcome comparing groups differing by 1 unit in pH having color red. Lastly, the "colorwhite:pH" parameter is interpreted as a difference in differences for those groups, i.e. the incremental change in the pH slope comparing whites to reds.

I think you should rewrite your color formula to remove ":" and replace them with "*"

> set.seed(1)
> a <- sample(letters[1:3], 100, replace=TRUE)
> b <- sample(LETTERS[1:3], 100, replace=TRUE)
> y <- rnorm(100)
> lm(y ~ a * b)

lm(formula = y ~ a * b)

(Intercept)           ab           ac           bB           bC        ab:bB  
     0.1684      -0.3894      -0.2614      -0.2807      -0.3981       0.8720  
      ac:bB        ab:bC        ac:bC  
     0.2099       0.6215       0.4547  

> lm(y ~ a : b) ## wrong

lm(formula = y ~ a:b)

(Intercept)        aa:bA        ab:bA        ac:bA        aa:bB        ab:bB  
   -0.03642      0.20484     -0.18456     -0.05654     -0.07587      0.40675  
      ac:bB        aa:bC        ab:bC        ac:bC  
   -0.12738     -0.19331      0.03883           NA  

> tapply(y, interaction(a, b), mean)
         a.A          b.A          c.A          a.B          b.B          c.B 
 0.168415779 -0.220978904 -0.092958625 -0.112286696  0.370329364 -0.163796519 
         a.C          b.C          c.C 
-0.229732632  0.002406992 -0.036419652 
| cite | improve this answer | |
  • $\begingroup$ Thanks, but I'm not understanding where this shows how to specify the interaction terms listed in my screenshot in a subsequent linear model. $\endgroup$ – Info5ek Feb 24 '14 at 22:30
  • $\begingroup$ You will never have interaction terms with the referent group, parsing factors out into a sequence of dummies will never have an indicator for a referent group, and consequently none of the product terms either. The : notation messes everything up, whereas * is much better. $\endgroup$ – AdamO Feb 24 '14 at 22:35
  • $\begingroup$ I haven't been taught about referent groups, nor indicators for referent groups so I don't understand your comment, but it seems like what you're saying is that there is no way for me to reference these interaction terms in a separate model. $\endgroup$ – Info5ek Feb 25 '14 at 0:08
  • $\begingroup$ Yes, inspect what happens under-the-hood when you parse a factor variable into multiple indicators. y <- factor(sample(letters[1:3], 10, replace=TRUE)); model.matrix(~y) $\endgroup$ – AdamO Feb 25 '14 at 1:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.