About a decade ago I complained that it was difficult to do contrast analysis following a factorial linear model using Stata. I cannot emphasize more strongly how this has changed with -contrast- and -margins, contrast-. It's now fantastic.
For example, I had a 2x2x4 3-way AxBxC factorial linear model, with some covariates. In Stata I got all 24 needed contrasts in 4 contrast commands using something like:
regress y i.A_B##i.C
contrast (r2.A_B@C) (r2.A_B#ar.C)
contrast (r3.A_B@C) (r3.A_B#ar.C)
contrast (r4.A_B@C) (r4.A_B#ar.C)
regress y i.A##i.B##i.C
contrast A#B@C
Because I wanted the contrast p-values to come from a simultaneous estimation procedure, I needed to use R (package multcomp, function glht()). However, common packages in R for contrast estimation stop at all pairwise or many-to-one contrasts (e.g., package lsmeans). Other contrast packages for interaction analysis don't do simultaneous estimation (package phia). In order to use R glht() for these partial interaction contrasts, I had to express my contrasts in terms of the regression coefficients. If you have ever done this by hand in a regression with16 coefficients, it's a little mind bending and error prone. You can produce exactly the same results for multiplicity unadjusted p-values in R and Stata as I did, but it required 50+ lines of R code to get there.
Kudos to Stata Corp for excellent contrast implementation!
For example, I had a 2x2x4 3-way AxBxC factorial linear model, with some covariates. In Stata I got all 24 needed contrasts in 4 contrast commands using something like:
regress y i.A_B##i.C
contrast (r2.A_B@C) (r2.A_B#ar.C)
contrast (r3.A_B@C) (r3.A_B#ar.C)
contrast (r4.A_B@C) (r4.A_B#ar.C)
regress y i.A##i.B##i.C
contrast A#B@C
Because I wanted the contrast p-values to come from a simultaneous estimation procedure, I needed to use R (package multcomp, function glht()). However, common packages in R for contrast estimation stop at all pairwise or many-to-one contrasts (e.g., package lsmeans). Other contrast packages for interaction analysis don't do simultaneous estimation (package phia). In order to use R glht() for these partial interaction contrasts, I had to express my contrasts in terms of the regression coefficients. If you have ever done this by hand in a regression with16 coefficients, it's a little mind bending and error prone. You can produce exactly the same results for multiplicity unadjusted p-values in R and Stata as I did, but it required 50+ lines of R code to get there.
Kudos to Stata Corp for excellent contrast implementation!
