3
$\begingroup$

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?

$\endgroup$
3
$\begingroup$

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:

enter image description here

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

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ 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?! $\endgroup$ – Stefan Dec 17 '15 at 4:00
  • $\begingroup$ @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! $\endgroup$ – Antoni Parellada Dec 17 '15 at 4:50
  • $\begingroup$ No problem Antoni :) If you found out more about which confidence limits to use, please let me know. Thanks! $\endgroup$ – Stefan Dec 17 '15 at 5:56

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.