9
$\begingroup$

I have a regression model that looks like this: $$Y = \beta_0+\beta_1X_1 + \beta_2X_2 + \beta_3X_3 +\beta_{12}X_1X_2+\beta_{13}X_1X_3+\beta_{123}X_1X_2X_3$$

...or in R notation: y ~ x1 + x2 + x3 + x1:x2 + x1:x3 + x1:x2:x3

Let's say $X_1$ and $X_2$ are categorical variables and $X_3$ is numeric. The complication is that $X_1$ has three levels $X_{1a}, X_{1b}, X_{1c}$ and instead of standard contrasts, I need to test:

  • Whether the intercept for level $X_{1a}$ significantly differs from the average intercept for levels $X_{1b}$ and $X_{1c}$.
  • Whether the response of $X_2$ is significantly different between level $X_{1a}$ and the average of levels $X_{1b}$ and $X_{1c}$.
  • Whether the slope of $X_3$ is significantly different between level $X_{1a}$ and the average of levels $X_{1b}$ and $X_{1c}$.

Based on this post it seems like the matrix I want is...

 2
-1
-1

So I do contrasts(mydata$x1)<-t(ginv(cbind(2,-1,-1))). The estimate of $\beta_1$ changes, but so do the others. I can reproduce the new estimate of $beta_1$ by subtracting the predicted values of the $X_1b$ and $X_1c$ group means (when $X_3=0$ and $X_2$ is at its reference level) from twice the value of $X_1a$ at those levels. But I can't trust that I specified my contrast matrix correctly unless I can also similarly derive the other coefficients.

Does anybody have any advice for how to wrap my head around the relationship between cell means and contrasts? Thanks. Is there a standard name for this type of contrast?


Aha! According to the link posted in Glen_b's answer, the bottom line is, you can convert ANY comparison of group means you want into an R-style contrast attribute as follows:

  1. Make a square matrix. The rows represent the levels of your factor and the columns represent contrasts. Except the first one, which tells the model what the intercept should represent.
  2. If you want your intercept to be the grand mean, fill the first column with all of the same non-zero value, doesn't matter what. If you want the intercept to be one of the level means, put a number in that row and fill the rest with zeros. If you want the intercept to be a mean of several levels, put numbers in those rows and zeros in the rest. If you want it to be a weighted mean, use different numbers, otherwise use the same number. You can even put in negative values in the intercept column and that probably means something too, but it completely changes the other contrasts, so I have no idea what that's for
  3. Fill in the rest of the columns with positive and negative values indicating what levels you want compared to what others. I forget why summing to zero is important, but adjust the values so that the columns do sum to zero.
  4. Transpose the matrix using the t() function.
  5. Use ginv() from the MASS package or solve() to get the inverse of the transposed matrix.
  6. Drop the first column, e.g. mycontrast<-mycontrast[,-1]. You now have a p x p-1 matrix, but the information you put in for your intercept was encoded in the matrix as a whole during step 5.
  7. If you want labels in the summary output more pleasant to read than lm() et al.'s default output, name your matrix's columns accordingly. The intercept will always automatically get named (Intercept) however.
  8. Make your matrix the new contrast for the factor in question, e.g. contrasts(mydata$myfactor)<-mymatrix
  9. Run lm() (and probably many other functions that use formulas) as normal in standard R without having to load glht, doBy, or contrasts.

Glen_b, thank you, and thank you UCLA Statistical Consulting Group. My applied stats prof spent several days handwaving on this topic, and I was still left with no idea how to actually write my own contrast matrix. And now, an hour of reading and playing with R, and I finally think I get it. Guess I should have applied to UCLA instead. Or University of StackExchange.

$\endgroup$
5
$\begingroup$

That comparison of one with the mean of all later variables is (aside from scale), called Helmert coding or Helmert contrasts. The one you give is the first contrast, the other would be a scaled version of $(0, 1, -1)^\top$.

What R calls helmert coding, this calls 'reverse Helmert'. They're equivalent up to a change of variable order.

$\endgroup$
  • $\begingroup$ Should equivalent "up to a change of variable order" should be "given a reversal of the order of levels"? So to get the SAS/SPSS type of helmert (compare to average of remaining levels), one would have to either reverse the order of the levels of the factor concerned, or rev() each column then rev() each row of the matrix returned by contra.helmert ? $\endgroup$ – tim Sep 1 '15 at 20:55
  • $\begingroup$ @tim "a change of variable order" isn't intended to mean "any change of variable order", but quite literally $\:$ a change of variable order (that is, one in particular). $\endgroup$ – Glen_b Sep 1 '15 at 23:41

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.