My understanding is that I can get the mean of f1High f2B by summing up the first three coefficients

> 2.467899
 b<-mean(dat[dat$f1=="High" & dat$f2=="B","y"])
> 2.992012

It's a perfectly balanced design, can someone please give me a clear explanation why a and b are different?

  • $\begingroup$ Perhaps if you created a minimal reproducible example, that would permit you to inspect modmat, the coefficients, and any other details, thereby revealing what's going on. You should have no trouble reducing modmat to a $4\times 4$ matrix, which will be accessible. In particular, take a look at model.matrix(mod) and compare it to modmat. $\endgroup$ – whuber Feb 7 at 22:10
  • $\begingroup$ @whuber, unfortunately my R skills are not as good as I wish they were (the example I provide was taken from the Internet). Anyway, I can tell that the difference between modmat and model.matrix(mod) is the inclusion of interactions in the former. As suggested in Greg Snow's answer, the addition of the interactions saturates the model, giving the same means. But I'm still confused why. The data are the same regardless of whether or not I include interactions in my model. This must relate to a statistical property that I'm failing to understand, but I suppose that is what I'm asking. $\endgroup$ – locus Feb 8 at 1:12

Getting the exact mean will only happen if the model is saturated, or the relationship is purely additive. Since you don't include an interaction in the fitted model, it is not saturated. Since you add random error that is not orthogonal to your model terms and your generating coefficients are not purely additive, it is not purely additive.

To see what you expect, include an interaction in the model fit, then add the appropriate coefficients (or create a perfectly additive model with orthogonal errors).

  • $\begingroup$ thank you that was helpful. I was aware that including interactions would get me the same mean, but wasn't sure why. So if I understood correctly the reason is that the model is not saturated. I'm relatively new to all this, could you please clarify what it means for a model to be saturated? The data are the same whether or not interactions are included, why would the means be different in a non-saturated model? Is the random error orthogonal to the model terms in the saturated model? $\endgroup$ – locus Feb 8 at 1:01
  • 1
    $\begingroup$ @locus, the simple definition of saturated is that there are as many coefficients being fit as there are means being predicted. If you fit the model including the interactions then you have 6 means to predict and 6 coefficients to fit to predict those means. But fitting without the interaction essentially forces 2 parameters to be 0, meaning that the method has to adjust the other 4 parameters to get your predictions as close as possible, but this moves them away from the exact fit. $\endgroup$ – Greg Snow Feb 8 at 16:21
  • $\begingroup$ thanks @Greg Snow. I read here in CV that a saturated model is in most cases a bad outcome. Does that mean that I shouldn't model every single interaction like in the example above in order to avoid a saturated model? $\endgroup$ – locus Feb 17 at 2:58
  • 1
    $\begingroup$ @locus, for regression models a saturated model will pretty much always be overfitting and should be avoided. A model predicting just 6 means is a bit different and saturation may be justified (but not automatically, you should look to see if the simpler model is "good enough"). $\endgroup$ – Greg Snow Feb 19 at 21:44

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.