Example 9.24: Changing the parameterization for categorical predictors

In our book, we discuss the important question of how to assign different parameterizations to categorical variables when fitting models (section 3.1.3). We show code in R for use in the lm() function, as follows:

lm(y ~ x, contrasts=list(x,"contr.treatment")

This works great in lm() and some other functions, notably glm(). But for functions from contributed packages, the contrasts option may not work.

Here we show a more generic approach to setting contrasts in R, using Firth logistic regression, which is discussed in Example 8.15, to demonstrate. This approach is also shown in passing in section 3.7.5.

R
We'll simulate a simple data set for logistic regression, then examine the results of the default parameterization.

n = 100
j = rep(c(0,1,2), each = n)
linpred = -2.5 + j
y = (runif(n*3) < exp(linpred)/(1 + exp(linpred)) )
library(logistf)
flrfactor = logistf(y ~ as.factor(j))
summary(flrfactor)
                    coef  se(coef)      Chisq            p
(Intercept)   -2.1539746 0.3276441        Inf 0.000000e+00
as.factor(j)1  0.3679788 0.4343756  0.7331622 3.918601e-01
as.factor(j)2  1.7936917 0.3855682 26.2224650 3.042623e-07

To see what R is doing, use the contrasts() function:

> contrasts(as.factor(j))
  1 2
0 0 0
1 1 0
2 0 1

R made indicator ("dummy") variables for two of the three levels, so that the estimated coefficients are the log relative odds for these levels vs. the omitted level. This is the "contr.treatment" structure (default for unordered factors). The defaults can be changed with options("contrasts"), but this is a sensible one.

But what if we wanted to assess whether a linear effect was plausible, independent of any quadratic effect? For glm() objects we could examine the anova() between the model with the linear term and the model with the linear and quadratic terms. Or, we could use the syntax shown in the introduction, but with "contr.poly" in place of "contr.treatment". The latter approach may be preferable, and for the logistf() function (and likely many other contributed functions) the contrasts = option does not work. In those cases, use the contrasts function:

jfactor = as.factor(j)
contrasts(jfactor) = contr.poly(3)
flrfc = logistf(y ~ jfactor)
summary(flrfc)
                  coef  se(coef)     Chisq            p
(Intercept) -1.4334177 0.1598591       Inf 0.000000e+00
jfactor.L    1.2683316 0.2726379 26.222465 3.042623e-07
jfactor.Q    0.4318181 0.2810660  2.472087 1.158840e-01

Not surprisingly, there is no need for a quadratic term, after the linear trend is accounted for. The canned contrasts available in R are somewhat limited--effect cell coding is not included, for example. You can assign contrasts(x) a matrix you write manually in such cases.

SAS
In SAS, the class statement for the logistic procedure allows many parametrizations, including "orthpoly", which matches the "contr.poly" contrast from R. However, most modeling procedures do not have this flexibility, and you would have to generate your contrasts manually in those cases, typically by creating new variables with the appropriate contrast values. Here we show the reference cell coding that is the default in R. Perversely, it is not the the default in proc logistic despite it being the only option in most procedures. On the other hand, it does allow the user to specify the reference category.

data test;
do i = 1 to 300;
  j = (i gt 100) + (i gt 200);
  linpred = -2.5 + j;
  y = (uniform(0) lt exp(linpred)/(1 + exp(linpred)) );
  output;
  end;
run;

title "Reference cell";
proc logistic data = test;
class j (param=ref ref='0');
model y(event='1') = j / firth clparm = pl;
run;

title "Polynomials";
proc logistic data = test;
class j (param=orthpoly);
model y(event='1') = j;
run;

With the results:

Reference cell
                                Standard         Wald
  Parameter     DF   Estimate      Error   Chi-Square   Pr > ChiSq
  Intercept      1    -2.6110     0.1252     434.6071       <.0001
  j         1    1     1.2078     0.1483      66.3056       <.0001
  j         2    1     2.2060     0.1409     245.1215       <.0001


Polynomials
                                  Standard        Wald
  Parameter         DF  Estimate     Error  Chi-Square  Pr > ChiSq
  Intercept          1   -1.4761    0.0540    746.6063      <.0001
  j         OPOLY1   1    0.9032    0.0577    245.3952      <.0001
  j         OPOLY2   1   -0.0502    0.0501      1.0029      0.3166

An unrelated note about aggregators
We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers and PROC-X with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.