# Dummy variables in multiple regression

Let's suppose I measure weights of eight male and eight female mice, from four litters, that were subjected to two different treatments: a and b

Here are my data:

set.seed(1)
df <- data.frame(sex = c(rep("m",8),rep("f",8)), treatment = rep(c(rep("a",4),rep("b",4)),2))
df$weight <- c(rnorm(4,1),rnorm(4,2),rnorm(4,3),rnorm(4,4)) df$litter <- rep(c(rep("l1",2),rep("l2",2),rep("l3",2),rep("l4",2)),2)


Now I fit a mixed effects model to the interaction between sex and treatment, which are categorical fixed effects, and litter is defined as a random effect:

df$sex = as.factor(df$sex)
df$treatment = as.factor(df$treatment)
df$litter = as.factor(df$litter)

fit <- lmer(weight ~ sex*treatment + (1|litter), data = df)


And here's the output:

summary(fit)
Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ sex * treatment + (1 | litter)
Data: df

REML criterion at convergence: 38.6

Scaled residuals:
Min       1Q   Median       3Q      Max
-1.47793 -0.54377 -0.01737  0.48467  1.45905

Random effects:
Groups   Name        Variance Std.Dev.
litter   (Intercept) 0.4318   0.6571
Residual             0.7574   0.8703
Number of obs: 16, groups:  litter, 4

Fixed effects:
Estimate Std. Error t value
(Intercept)      3.54300    0.63662   5.565
sexm            -2.46379    0.61541  -4.004
treatmentb       0.01801    0.90031   0.020
sexm:treatmentb  1.08648    0.87032   1.248

Correlation of Fixed Effects:
(Intr) sexm   trtmnt
sexm        -0.483
treatmentb  -0.707  0.342
sxm:trtmntb  0.342 -0.707 -0.483


Obviously the "f" sex and the "a" treatment are set as dummy variables to a zero baseline. My question relates to the interpretation: is there a way to extract the effect size of the dummy variables? Alternatively, if I add an additional sex and treatment variable, for which I set weight to 0, will they therefore serve as dummy variables that really have a zero effect size?

## 1 Answer

You can use the R package {lmerTest} and call on the following function:

summary(fit)
Linear mixed model fit by REML t-tests use Satterthwaite approximations to
degrees of freedom [lmerMod]
Formula: weight ~ sex * treatment + (1 | litter)
Data: df

REML criterion at convergence: 38.6

Scaled residuals:
Min       1Q   Median       3Q      Max
-1.47793 -0.54377 -0.01737  0.48467  1.45905

Random effects:
Groups   Name        Variance Std.Dev.
litter   (Intercept) 0.4318   0.6571
Residual             0.7574   0.8703
Number of obs: 16, groups:  litter, 4

Fixed effects:
Estimate Std. Error       df t value Pr(>|t|)
(Intercept)      3.54300    0.63662  3.34300   5.565  0.00853 **
sexm            -2.46379    0.61541 10.00000  -4.004  0.00250 **
treatmentb       0.01801    0.90031  3.34300   0.020  0.98518
sexm:treatmentb  1.08648    0.87032 10.00000   1.248  0.24033
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) sexm   trtmnt
sexm        -0.483
treatmentb  -0.707  0.342
sxm:trtmntb  0.342 -0.707 -0.483


Under random effects, the litter Std. Dev. of 0.6571 measures the variability between litters, as opposed to the variability that is still unaccounted for (0.8703). This can be visually inspected with boxplots:

Among the fixed effects the differences will simply be in intercepts, also illustrated with boxplots:

Notice on the output that the actual intercepts should be calculated as:

fixef(fit)
(Intercept)            sexm      treatmentb sexm:treatmentb
3.54300434     -2.46379392      0.01800987      1.08647799


so that the (Intercept) in the output (i.e. $\small3.543$) is the average of females treated with a:

> a_fem <- df[df$sex=="f" & df$treatment=="a",]; mean(a_fem$weight) [1] 3.543004  The second coefficient is the average for males treated with a ($\small 3.543 - 2.464 = 1.079$): > a_male <- df[df$sex=="m" & df$treatment=="a",]; mean(a_male$weight) [1] 1.07921


The third coefficient is for females treated with b, amounting to ($\small 3.543 + 0.018 = 3.561$):

> b_fem <- df[df$sex=="f" & df$treatment=="b",]; mean(b_fem$weight) [1] 3.561014  Finally the last coefficient stands for males treated with b or ($\small 3.543−2.464+0.018+1.083=2.187$): > b_male <- df[df$sex=="m" & df$treatment=="b",]; mean(b_male$weight) [1] 2.183698


We see in the output that there is a significant difference in intercepts in females treated with a and males treated with a.

The fixed effects coefficients come with the standard errors calculated. Otherwise, they can be extracted from the variance-covariance matrix:

vcov(fit)
4 x 4 Matrix of class "dpoMatrix"
(Intercept)       sexm treatmentb sexm:treatmentb
(Intercept)       0.4052795 -0.1893622 -0.4052795       0.1893622
sexm             -0.1893622  0.3787245  0.1893622      -0.3787245
treatmentb       -0.4052795  0.1893622  0.8105590      -0.3787245
sexm:treatmentb   0.1893622 -0.3787245 -0.3787245       0.7574489


For example, for females treated with b:

(fem_b_SE <- sqrt(diag(vcov(fit)))[3]) [1] 0.9003105

The confidence intervals can be calculated as follows for example in the case of female treated with b:

(fem_b_SE <- sqrt(diag(vcov(fit)))[3]) [1] 0.9003105
(fem_b_coef <- fixef(fit)[1] + fixef(fit)[3])
3.561014
(fem_b_CI <- fem_b_coef + c(-1, 1) * qnorm(0.975) * fem_b_SE)
[1] 1.796438   5.325590


or in the case of males treated with b:

(male_b_SE <- sqrt(diag(vcov(fit)))[4]) [1] 0.8703154
(male_b_coef <- sum(fixef(fit)[1:4]))
2.183698
(male_b_CI <- male_b_coef + c(-1, 1) * qnorm(0.975) * male_b_SE)
[1] 0.4779114   3.8894852


Notice that the function built into lme4 gives different values:

confint(fit)
Computing profile confidence intervals ...
2.5 %    97.5 %
.sig01           0.0000000  1.308364
.sigma           0.5581087  1.251771
(Intercept)      2.4416191  4.644389
sexm            -3.6591640 -1.268425
treatmentb      -1.5518958  1.589257
sexm:treatmentb -0.6040305  2.776985


are significantly different, and give a bunch of error messages.

• Antoni,in your last step where you derive the mean value for sexm:treatmentb, you are missing the estimate for treatmentb,so instead of $\small 3.543−2.46+1.086=2.165$, it should be $\small 3.543−2.464+0.018+1.083=2.187$, which also matches your calculated mean value. If that is correct than the example confidence interval should also reflect this.Also if you ran lsmeans::lsmeans(fit, pairwise~sex:treatment), you will get the same mean values you calculated and lower and upper CLs.However they are different from yours and also different from confint(fit).Not sure which way to go here?! – Stefan Dec 17 '15 at 4:00
• @Stefan Thank you very much. I didn't realize the calculated values were different, which had tipped me off to the coefficient I had forgotten to include. It's now fixed. I'm looking into your other comments. It's great you paid attention to potentially misleading info! – Antoni Parellada Dec 17 '15 at 4:50
• No problem Antoni :) If you found out more about which confidence limits to use, please let me know. Thanks! – Stefan Dec 17 '15 at 5:56