**SAS and R**, and kindly contributed to R-bloggers)

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.

**leave a comment**for the author, please follow the link and comment on his blog:

**SAS and R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...