When running a mixed-effects model with categorical predictors, you may wish to test the fixed effects of the model. When your model includes categorical variables with three or more levels or interactions, this requires a multiple degrees of freedom test. In other software packages like SAS, Type III tests of fixed effects are presented with the regression output. In R, this is not the case. However, we can use contrast and ANOVA-type commands to extract these effects.
We will use the dataset hsbdemo and the R packages foreign (to read in the data) and nlme (to run a mixed-effects model). We will run a random-intercept model where read is predicted by female, prog, and math.
library(foreign) library(nlme)d1 <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta") attach(d1)
In order to use contrast commands, we will need to convert our string factor variables to integer factor variables starting at 1.
female12 <- as.factor(1*(female=="male") + 2*(female=="female"))prog123 <- as.factor(1*(prog=="general") + 2*(prog=="academic") + 3*(prog=="vocational"))
Next, before running the model, we define the contrasts for our factor variables. For a simple model without interactions, we can use the contr.treatment command which also allows us to specify a baseline (which is helpful if you are attempting to match model results seen in a publication or another software). We provide the number of levels and the level we wish to use as the baseline. This command generates a corresponding matrix that we then indicate as the contrast matrix for a given variable using the contrasts command. These matrices correspond to “dummy coding” the categorical variables. For further information on dummy coding, see FAQ: What is dummy coding?.
a <- contr.treatment(2, base = 2, contrasts = TRUE) class(a) [1] "matrix" a 1 1 1 2 0 b <- contr.treatment(3, base = 3, contrasts = TRUE) b 1 2 1 1 0 2 0 1 3 0 0 contrasts(female12) <- a contrasts(prog123) <- b
We can now run our model and see that the coefficients of the model are reported according to the baselines indicated in these contrasts.
m1.lme <- lme(read ~ female12 + prog123 + math, random =~ 1|cid, method = "ML")
summary(m1.lme)
Linear mixed-effects model fit by maximum likelihood
Data: NULL
AIC BIC logLik
1289.254 1312.342 -637.6271
Random effects:
Formula: ~1 | cid
(Intercept) Residual
StdDev: 11.00646 4.809769
Fixed effects: read ~ female12 + prog123 + math
Value Std.Error DF t-value p-value
(Intercept) 70.08683 4.776200 176 14.674184 0.0000
female121 1.35292 0.731138 176 1.850427 0.0659
prog1231 -2.10673 1.006942 176 -2.092208 0.0379
prog1232 -2.04756 0.968723 176 -2.113675 0.0360
math -0.33022 0.074898 176 -4.408968 0.0000
Correlation:
(Intr) fml121 pr1231 pr1232
female121 -0.035
prog1231 -0.231 -0.022
prog1232 -0.223 -0.001 0.415
math -0.843 -0.039 0.195 0.188
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.5019691 -0.6206737 0.0354358 0.7148373 2.3331462
Number of Observations: 200
Number of Groups: 20
To see the Type III tests of fixed effects, we use the anova.lme command. To indicate Type III tests, we provide type = “marginal”.
anova.lme(m1.lme, type = "marginal", adjustSigma = F)numDF denDF F-value p-value (Intercept) 1 176 220.85299 <.0001 female12 1 176 3.51188 0.0626 prog123 2 176 3.20534 0.0429 math 1 176 19.93743 <.0001
If we have a more complex model that includes interactions of categorical variables, dummy coding no longer captures the main effects of these categorical variables but rather the simple effects. To test for main effects in our Type III tests, we will need to use a different contrast command that defines our categorical variables using “effect coding”. The contrast matrices under effect coding are no longer composed of only zeroes and ones; each column in the contrast matrix now sums to zero. We can create these effect coding matrices in R with the contr.sum function.
contrasts(female12) <- contr.sum contrasts(female12) [,1] 1 1 2 -1 contrasts(prog123) <- contr.sum contrasts(prog123) [,1] [,2] 0 1 0 1 0 1 2 -1 -1
For more information on effect coding, see FAQ: What is effect coding?. We can run our model to estimate coefficients based on this coding, then generate the desired tests after running the model.
m1.lme2 <- lme(read ~ female12*prog123 + math, random =~ 1|cid, method = "ML")
summary(m1.lme2)
Linear mixed-effects model fit by maximum likelihood
Data: NULL
AIC BIC logLik
1288.429 1318.114 -635.2145
Random effects:
Formula: ~1 | cid
(Intercept) Residual
StdDev: 10.88791 4.751459
Fixed effects: read ~ female12 * prog123 + math
Value Std.Error DF t-value p-value
(Intercept) 68.77214 4.674240 174 14.713011 0.0000
female121 0.94010 0.386053 174 2.435143 0.0159
prog1231 -0.69864 0.613890 174 -1.138046 0.2567
prog1232 -0.62045 0.590308 174 -1.051065 0.2947
math -0.31881 0.075398 174 -4.228375 0.0000
female121:prog1231 -0.02636 0.566363 174 -0.046538 0.9629
female121:prog1232 0.98816 0.584777 174 1.689812 0.0929
Correlation:
(Intr) fml121 pr1231 pr1232 math f121:1231
female121 0.035
prog1231 -0.064 0.001
prog1232 -0.076 0.009 -0.581
math -0.844 -0.033 0.093 0.110
female121:prog1231 0.129 0.123 0.137 -0.097 -0.149
female121:prog1232 -0.127 0.163 -0.072 0.086 0.148 -0.644
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.40201967 -0.66416223 0.05083885 0.71444086 2.32978676
anova.lme(m1.lme2, type = "marginal", adjustSigma = F)
numDF denDF F-value p-value
(Intercept) 1 174 224.32402 <.0001
female12 1 174 6.14499 0.0141
prog123 2 174 2.96372 0.0542
math 1 174 18.52762 <.0001
female12:prog123 2 174 2.44197 0.0900
The table output from the anova.lme command contains the Type III tests of fixed effects for this model.