Categorical variable coding to compare all levels to all levels

Question

I am trying to determine the best coding system for my categorical variables to use in a regression with categorical and continuous variables. I have been using this page as a resource but none of the coding methods seem to match what I want to do: http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm

I have reason to expect that there could be differences among all four levels of a categorical variable and therefore would like to code it to perform all pairwise comparisons (1vs2, 1vs3, 1vs4, 2vs3, 2vs4, 3vs4,). I thought that a Sum (deviation) method would accomplish this as it says it compares each level to the grand mean. But in my model I only see terms for levels 1, 2, 3, and do not have a term for level 4 of the variable. Is this level incorporated into the intercept? If so is there another coding system I should use to get at the comparisons I am looking for?

Answer 1

You want all possible pairwise comparisons of levels, but there are much more pairs than there are degrees of freedom in the factor. Say the factor has five levels, then you need 4 parameters to code it, but there are $\binom{5}{2}$ pairs, that is, 10 pairs. So it is impossible to find a coding with one parameter for each comparison.

The solution is to use whatever coding you wants, and then compute the 10 pairwise contrasts afterwards, after estimating the model, from the model output. In R, for instance, this could be done many ways , either "by hand", or with the use of packages like contrast or multcomp.

Below an R example, done "by hand", for confidence intervals of all pairwise comparisons:

xfac  <-  factor(rep(1:5, each=10))
y     <-  rnorm(50, mean=c(rep(0, 20), rep(1, 30)), sd=2)

mod   <-  lm( y ~ 0 + xfac)

# generating a hypothesis contrasts matrix with 10 rows:
# each row is one contrast:      
cmat  <-  matrix(0, 10, 5)
nam   <-  character(length=10)   
row   <-  0
for (i in 1:4) for (j in (i+1):5)   {
                row  <-  row+1
                nam[row]  <-  paste("x[", i, "]-x[", j, "]", 
                                     sep="")
                cmat[row, c(i, j)]  <-  c(1, -1)
               }
rownames(cmat)  <-  nam  

# We write a contrast testing function by hand:
my.contrast  <-  function(mod,  cmat)  {
    co  <-  coef(mod)
    CV  <-  vcov(mod)
    se  <-  sqrt( diag( cmat %*% CV %*% t(cmat) ))
    df  <-  mod$df.residual
    contr  <-  cmat %*% co
    ul  <-  qt(0.975, df=df)
    ci  <-  cbind(contr-ul*se, contr+ul*se)
    ci
}

And then using it gives the result:

> my.contrast(mod, cmat)
               [,1]      [,2]
x[1]-x[2] -1.946376 1.7921298
x[1]-x[3] -3.044916 0.6935897
x[1]-x[4] -2.136283 1.6022227
x[1]-x[5] -2.301393 1.4371135
x[2]-x[3] -2.967793 0.7707130
x[2]-x[4] -2.059160 1.6793460
x[2]-x[5] -2.224269 1.5142368
x[3]-x[4] -0.960620 2.7778861
x[3]-x[5] -1.125729 2.6127769
x[4]-x[5] -2.034362 1.7041439

Answer 2

In the case of ANOVA regression, when you have all the categorical variables, usually one of those will be represented by the intercept (the one you choose). You can verify this: the intercept should be the mean of the variable that was "left out" from your model.

the logic of the variable that is represented in the regression is due to the following logic.

• Dummy(b) = 1 or 0
• Dummy(c) = 1 or 0
• Dummy(d) = 1 or 0

then

• Dummy(a) = b0 since all others are zero.

Therefore y = b0 + b1 * b + b2 * c + b3 * d; if all others are zero the y = b0, where b0 is the intercept and the mean of the first variable.

Hope this helps.

J