11
$\begingroup$

If a factor variable (e.g. gender with levels M and F) is used in the glm formula, dummy variable(s) are created, and can be found in the glm model summary along with their associated coefficients (e.g. genderM)

If, instead of relying on R to split up the factor in this way, the factor is encoded in a series of numeric 0/1 variables (e.g. genderM (1 for M, 0 for F), genderF (1 for F, 0 for M) and these variables are then used as numeric variables in the glm formula, would the coefficient result be any different?

Basically the question is: does R use a different coefficient calculation when working with factor variables versus numeric variables?

Follow-up question (possibly answered by the above): besides just the efficiency of letting R create dummy variables, is there any problem with re-coding factors as a series of numeric 0,1 variables and using those in the model instead?

$\endgroup$
  • 2
    $\begingroup$ The coefficients will be the same, R's default coding for factors is exactly as you described (this is called "dummy" coding). If doing this manually, you will need to aware of the the "dummy variable trap" - you cannot include all $N$ created variables, but only can include $N-1$ (alternately, exclude the intercept); otherwise your model is overidentified. There are also other methods of encoding factor variables as well. $\endgroup$ – Affine Apr 16 '14 at 15:34
  • $\begingroup$ @Affine my guess is that if I fed in both genderM and genderF, one of them would return NA for coefficients (with the message that a variable was excluded because of singularities). This makes sense because they perfectly linearly related in this case. But you say that I cannot include all N; does that mean, even though genderF is set to NA, it would result in differences/problems for the genderM coefficient? Or, more simply, if GLM/LM excludes variables because of singularities, is using an overidentified model a problem? (I agree with your point - just querying the practical ramifications) $\endgroup$ – Bryan Apr 16 '14 at 15:45
20
$\begingroup$

Categorical variables (called "factors" in R) need to be represented by numerical codes in multiple regression models. There are very many possible ways to construct numerical codes appropriately (see this great list at UCLA's stats help site). By default, R uses reference level coding (which R calls "contr.treatment"), and which is pretty much the default statistics-wide. This can be changed for all contrasts for your entire R session using ?options, or for specific analyses / variables using ?contrasts or ?C (note the capital). If you need more information about reference level coding, I explain it here: Regression based for example on days of the week.

Some people find reference level coding confusing, and you don't have to use it. If you want, you can have two variables for male and female; this is called level means coding. However, if you do that, you will need to suppress the intercept or the model matrix will be singular and the regression cannot be fit as @Affine notes above and as I explain here: Qualitative variable coding leads to singularities. To suppress the intercept, you modify your formula by adding -1 or +0 like so: y~... -1 or y~... +0.

Using level means coding instead of reference level coding will change the coefficients estimated and the meaning of the hypothesis tests that are printed with your output. When you have a two level factor (e.g., male vs. female) and you use reference level coding, you will see the intercept called (constant) and only one variable listed in the output (perhaps sexM). The intercept is the mean of the reference group (perhaps females) and sexM is the difference between the mean of males and the mean of females. The p-value associated with the intercept is a one-sample $t$-test of whether the reference level has a mean of $0$ and the p-value associated with sexM tells you if the sexes differ on your response. But if you use level means coding instead, you will have two variables listed and each p-value will correspond to a one-sample $t$-test of whether the mean of that level is $0$. That is, none of the p-values will be a test of whether the sexes differ.

set.seed(1)
y    = c(    rnorm(30), rnorm(30, mean=1)         )
sex  = rep(c("Female",  "Male"          ), each=30)
fem  = ifelse(sex=="Female", 1, 0)
male = ifelse(sex=="Male", 1, 0)

ref.level.coding.model   = lm(y~sex)
level.means.coding.model = lm(y~fem+male+0)

summary(ref.level.coding.model)
# ...
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.08246    0.15740   0.524    0.602    
# sexMale      1.05032    0.22260   4.718 1.54e-05 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ...
summary(level.means.coding.model)
# ...
# Coefficients:
#      Estimate Std. Error t value Pr(>|t|)    
# fem   0.08246    0.15740   0.524    0.602    
# male  1.13277    0.15740   7.197 1.37e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ...
$\endgroup$
  • 1
    $\begingroup$ Thanks for the addition of the code: this clearly demonstrates that the intercept in reference cell coding = sexFemale in cell means coding (I was additionally confused by "what happened" to the N-1 variable...) Maybe needs another question, but this is easy to understand with one variable, what about 2 or more? e.g. age: "old" "young". In ref cell coding, would coefficients show for sexMale, ageYoung (for example) and the intercept account for BOTH sexFemale and ageOld? $\endgroup$ – Bryan Apr 16 '14 at 17:31
  • 1
    $\begingroup$ If you have a factor w/ 3 levels, the intercept is the mean of the ref level, & the other 2 will be represented in the output. Their coefficients will both be the difference b/t them & the ref level & their ps will be the significance of those difs. If you have 2 factors, both will have ref levels, & the intercept will be the mean of those people who are in both ref groups (eg, young F) & the other levels will be dif b/t the given level of factor 1 w/ the ref level of the other factor & the both ref levels group. Eg old is old F - ` young F, & M` is young M - young F. $\endgroup$ – gung Apr 16 '14 at 18:19
  • 1
    $\begingroup$ I played a bit with this and experienced substantial R^2 difference between both approaches. I know it's only an R^2 , but is there an explanation for that? $\endgroup$ – hans0l0 Jul 17 '15 at 18:46
  • $\begingroup$ @hans0l0, I have no idea. There should be no difference. $\endgroup$ – gung Jul 17 '15 at 18:56
  • $\begingroup$ @gung but there is a difference. no doubt about that. Have a look at your own reproducible example. Somehow you just did not print it. I mean it's also a different F statistic and degrees of freedom do also vary between both setups. So it's only mildly surprising but I struggle to come up with a good explanation against those who would simply prefer the setup with the higher R^2 $\endgroup$ – hans0l0 Jul 18 '15 at 1:19
2
$\begingroup$

The estimated coefficients would be the same subject to the condition that you create your dummy variables (i.e. the numerical ones) consistent to R. For example: lets' create a fake data and fit a Poisson glm using factor. Note that gl function creates a factor variable.

> counts <- c(18,17,15,20,10,20,25,13,12)
> outcome <- gl(3,1,9)
> outcome
[1] 1 2 3 1 2 3 1 2 3
Levels: 1 2 3
> class(outcome)
[1] "factor"
> glm.1<- glm(counts ~ outcome, family = poisson())
> summary(glm.1)

Call:
glm(formula = counts ~ outcome, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   3.0445     0.1260  24.165   <2e-16 ***
outcome2     -0.4543     0.2022  -2.247   0.0246 *  
outcome3     -0.2930     0.1927  -1.520   0.1285    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

Since outcome has three levels, I create two dummy variables (dummy.1=0 if outcome=2 and dummy.2=1 if outcome=3) and refit using these numerical values:

> dummy.1=rep(0,9)
> dummy.2=rep(0,9)
> dummy.1[outcome==2]=1
> dummy.2[outcome==3]=1
> glm.2<- glm(counts ~ dummy.1+dummy.2, family = poisson())
> summary(glm.2)

Call:
glm(formula = counts ~ dummy.1 + dummy.2, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   3.0445     0.1260  24.165   <2e-16 ***
dummy.1      -0.4543     0.2022  -2.247   0.0246 *  
dummy.2      -0.2930     0.1927  -1.520   0.1285    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

As you can see the estimated coefficients are the same. But you need to be careful when creating your dummy variables if you want to get the same result. For example if I create two dummy variables as (dummy.1=0 if outcome=1 and dummy.2=1 if outcome=2) then the estimated results are different as follow:

> dummy.1=rep(0,9)
> dummy.2=rep(0,9)
> dummy.1[outcome==1]=1
> dummy.2[outcome==2]=1
> glm.3<- glm(counts ~ dummy.1+dummy.2, family = poisson())
> summary(glm.3)

Call:
glm(formula = counts ~ dummy.1 + dummy.2, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.7515     0.1459   18.86   <2e-16 ***
dummy.1       0.2930     0.1927    1.52    0.128    
dummy.2      -0.1613     0.2151   -0.75    0.453    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

This is because when you add outcome variable in glm.1, R by default creates two dummy variables namely outcome2 and outcome3 and defines them similarly to dummy.1 and dummy.2 in glm.2 i.e. the first level of outcome is when all other dummy variables (outcome2 and outcome3) are set to be zero.

$\endgroup$
  • $\begingroup$ Thanks for the code proof of the estimated coefficients being the same. Also the warning on creating my own is useful: I wanted to create my own because then a model variable would tie back directly to a database column by name (which could be useful downstream), but looks like I need to understand the issues with how I go about doing this. $\endgroup$ – Bryan Apr 16 '14 at 17:36

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.