Mixed Models: Multiple Random Parameters
This article is part of the Mixed Model series. For a list of topics covered by this series, see the Introduction article. If you're new to mixed models we highly recommend reading the articles in order.
Overview
This article includes an example models which contain both a random intercept and random slope. The article also introduces nested and crossed random effects and provides an example of each.
Random slopes
When a slope is random, the intercept may or may not be random as well. The gmm model, from prior articles, includes a random intercept which we accepted as significant. We will add a random slope for the x2 variable to the gmm model. The following code does a summary of the gmm to remind us of the details of the gmm model.
Enter the following command in your script and run it.
summary(gmm, correlation = FALSE, show.resids = FALSE)The results of the above command are shown below.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod] Family: binomial ( logit ) Formula: bin ~ x1 + x2 + (1 | g1) Data: pbDat AIC BIC logLik deviance df.resid 113.0 123.4 -52.5 105.0 96 Random effects: Groups Name Variance Std.Dev. g1 (Intercept) 4.255 2.063 Number of obs: 100, groups: g1, 12 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.1651 0.6618 0.249 0.803024 x1 0.4996 0.2983 1.675 0.093919 . x2 -1.5155 0.3949 -3.837 0.000124 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The gmm model has 12 groups associated with the g1 variable. The g1 variable is random, which results in a mean intercept and a standard deviation for the intercept. There are also two fixed continuous variables, x1 and x2. This provides a fixed slope for each, although the slope for x1 may be 0.
Adding a random slope for x2 will allow for different x2 slopes for each group in g1. These random slopes may or may not be correlated with the random intercepts already associated with g1. We will start the investigation of the random slope for x2 with the assumption that it is uncorrelated with the random intercept. The following code adds the uncorrelated random slope.
Enter the following commands in your script and run them.
gmmSx2 <- glmer(bin~x1 + x2 + (1|g1) + (0 + x2|g1), family=binomial, data=pbDat) print(summary(gmmSx2, correlation = FALSE), show.resids = FALSE)The results of the above commands are shown below.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod] Family: binomial ( logit ) Formula: bin ~ x1 + x2 + (1 | g1) + (0 + x2 | g1) Data: pbDat AIC BIC logLik deviance df.resid 109.3 122.4 -49.7 99.3 95 Random effects: Groups Name Variance Std.Dev. g1 (Intercept) 3.001 1.732 g1.1 x2 5.040 2.245 Number of obs: 100, groups: g1, 12 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.1713 0.6245 0.274 0.7839 x1 0.6853 0.3410 2.009 0.0445 * x2 -1.7150 0.9071 -1.891 0.0587 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The mean x2 slope for the g1 groups is estimated as -1.715 with a standard deviation of 5.0399. The mean x2 slope is just out side the normal significance level of .05. This loss of significance is due to the estimated standard error for x2 more than doubling in this model. To know if adding the group level slope for x2 improved the model, we need a p-value for the significance of the between group variance of the x2 slope. We will use a parametric bootstrap to estimate the p-value for this variance being greater than zero.
Enter the following commands in your script and run them.
pbnmSx2 <- pbnm(gmmSx2,gmm,nsim=1000,tasks=10,cores=2,seed=7894523) summary(pbnmSx2)The results of the above commands are shown below.
Parametric bootstrap testing: x2 | g1 = 0 from: glmer(formula = bin ~ x1 + x2 + (1 | g1) + (0 + x2 | g1), data = pbDat, family = binomial) 1000 samples were taken Thu Aug 18 14:44:28 2016 67 samples had warnings, 60 in alternate model 20 in null model 67 unused samples. 0.004 <= P(abs(x2 | g1) > |5.039862|) <= 0.071
As was seen in prior articles, there are a number of the simulation runs which fail to converge. This is possibly due to not having enough observations for the model form being fit. The data set contains only 100 observations and there are 12 groups. This results in a small number of observations to use for any group level estimate. We will continue our investigation of the random slope for x2, accepting that there is a concern that there are not enough observations. The p-value for the test of significance is greater than or equal to 0.004 and less than or equal to 0.071. The test is inconclusive.
We will now include the correlation between the random slope of x2 and the random intercept in the model.
Enter the following commands in your script and run them.
gmmCSx2 <- glmer(bin~x1 + x2 + (x2|g1), family=binomial, data=pbDat) print(summary(gmmCSx2, correlation = FALSE), show.resids = FALSE)The results of the above commands are shown below.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod] Family: binomial ( logit ) Formula: bin ~ x1 + x2 + (x2 | g1) Data: pbDat AIC BIC logLik deviance df.resid 111.0 126.7 -49.5 99.0 94 Random effects: Groups Name Variance Std.Dev. Corr g1 (Intercept) 3.162 1.778 x2 5.547 2.355 -0.27 Number of obs: 100, groups: g1, 12 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.2952 0.6773 0.436 0.6630 x1 0.6900 0.3444 2.004 0.0451 * x2 -1.7664 0.9490 -1.861 0.0627 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The estimated correlation between the intercepts and x2 slopes is -0.2657. We will use a parametric bootstrap to estimate the p-value for this correlation being different than zero.
Enter the following commands in your script and run them.
pbnmCSx2 <- pbnm(gmmCSx2,gmmSx2,nsim=1000,tasks=10,cores=2,seed=13573745) summary(pbnmCSx2)The results of the above commands are shown below.
Parametric bootstrap testing: (Intercept)#x2 | g1 = 0 from: glmer(formula = bin ~ x1 + x2 + (x2 | g1), data = pbDat, family = binomial) 1000 samples were taken Thu Aug 18 14:49:33 2016 2 samples had errors, 2 in alternate model 2 in null model 74 samples had warnings, 56 in alternate model 21 in null model 76 unused samples. 0.653 <= P(abs((Intercept)#x2 | g1) > |-0.2656907|) <= 0.729
The p-value is greater than 0.653. We retain the null assumption that the g1 group means and slopes for x2 are not correlated. The models to consider are gmm with no random x2 slope and gmmSx2 with a random slope for each g1 group. The evidence is not clear which is the better representation of the data. One would need to rely on science of the problem, which we do not have here, to make an informed judgment on selecting between these two models. Regardless of which of these two models is selected, comments on the small amount of data available would need to be reported as a strong caution with the model results.
Exercises
Use the models from the sleepstudy data set for the following exercises.
Add a random slope for the Days variable to your linear sleep study model. Is the random slope for Days significant?
Add a correlation parameter for the random Days and intercept. Is the correlation parameter significant?
Nested and crossed effects
A categorical variable, say L2, is said to be nested with another categorical variable, say, L3, if each level of L2 occurs only within a single level of L3. variables are crossed if the levels of of one random variable, say R1, occur within multiple levels of a second random variable, say R2. As an example, consider boxes of products packaged on shipping pallets. If all the boxes on a shipping pallet contain the same product, the shipping pallet identifier variable is nested within product type. If instead the boxes on a shipping pallet contained different products, then the shipping pallet identifier variable is crossed with product type.
The nesting level is the number of levels of nesting associated with a variable. Nesting level 1 is the individual observations. Nesting level 2 is associated with a grouping variable when multiple observations are associated with at least one of its levels. Nesting level 3 is associated with a grouping variable which has a nesting level 2 variable nested within its levels. Thus the variable levels of a nesting level 3 variable has the observations associated with the nesting level 2 grouping levels nested within its levels as well as the nesting level 2 variable's levels. Nesting levels continue to count up for each higher grouping of the data. The variable at each nesting level maybe modeled as fixed or random effects. There is no effect associated with nesting level 1. There are effects associated with higher nesting levels.
The table below provides an example of nested and crossed variables. The Lev2 variable is nested within the Lev3 variable. Observations 1, 2, and 3 are collected in effect level A in the Lev2 variable and observations 4, 5, and 6 in level B. Effect levels A and B are collected in effect level J of the Lev3 variable. Lev3 is at nesting level 3 since it has 2 nesting levels below it, observations and Lev2's levels. The variable Cross is crossed with Lev2 and Lev3. It is crossed since Cross's level P is contained in Lev2's levels A, B, C, and D and in J and K of Lev3. Cross is a nesting level 2 variable since it has observations nested within it's levels.
| obs | Lev2 | Lev3 | Cross |
|---|---|---|---|
| 1 | A | J | P |
| 2 | A | J | Q |
| 3 | A | J | R |
| 4 | B | J | P |
| 5 | B | J | Q |
| 6 | B | J | R |
| 7 | C | K | P |
| 8 | C | K | Q |
| 9 | C | K | R |
| 10 | D | K | P |
| 11 | D | K | Q |
| 12 | D | K | R |
Some caution is needed when determining if variables are crossed or not. A variable r1 may take a fixed set of values, such as A, B, C and D. These values may be repeated within the levels of variable r2. The levels of r1 may represent four levels which are applied consistently across the levels of r2, but the levels of r1 may instead represent four samples for each of the levels of L2 which are labeled the same for convenience. These four samples may not have any connection from one level of r2 to another level. Although r1 and r2 may appear crossed in this case, they are hierarchical. You must understand the data to make an informed choice on nested verses crossed in these situations.
As will be seen later in this article, mis-specifying a variable as crossed or nested can cause estimation errors.
The following is an example of specifying nested random effects.
The example will use the following variables.
- A: factor with 15 levels
- B: factor with 25 levels
- C: numeric
- y: numeric
- y ~ C + (1|A) + (1|A:B) results in the following model parameters
- (intercept) (mean intercept associate with the groups of A and A:B)
- slope effect associated with C
- variance of intercept associated with the groups of A
- variance of intercept associated with the groups of A:B
- variance of residuals
The random effect B is nested in the random effect A. The population is the unique levels of A interacted with B.
Crossed random effect example
The pbDat data set does not contain crossed and nested random effects. We will generate a data set which contains three random variables, r1, r2, and r3. The data set will also contain two response variables, yc (effects of r1 and r2 crossed) and yn (effects of r3 nested within the effects of r1.).
Enter the following commands in your script and run them.
set.seed(9873754) n <- 2000 n1 <- 20 n2 <- 25 x1 <- rnorm(n) r1 <- cut(runif(n), breaks = c(-1,seq(1/n1,(n1-1)/n1,1/n1),2), labels = paste0("a", seq(1,n1) ) ) r2 <- cut(runif(n), breaks = c(-1,seq(1/n2,(n2-1)/n2,1/n2),2), labels = paste0("a", seq((n1+1),(n1+n2)) ) ) r3 <- factor(r1:r2) error <- rnorm(n,0,.1) effectsR1 <- rnorm(length(levels(r1)))[as.numeric(r1)] effectsR2c <- rnorm(length(levels(r2)))[as.numeric(r2)] effectsR2n <- rnorm(length(levels(r3)))[as.numeric(r3)] yc <- 3 - 1*x1 + error + .6*effectsR1 - .3*effectsR2c yn <- 3 - 1*x1 + error + .6*effectsR1 - .3*effectsR2n datUnsorted <- data.frame(yc,yn,x1,r1,r2,r3,error) dat <- datUnsorted[order(r1,r2),]There are no display results for the above commands.
Crossed effects are assumed to be independent of each other and are specified as separate random terms in the model formula. The following code builds the crossed model.
Enter the following commands in your script and run them.
rc <- lmer(yc ~ x1 + (1|r1) + (1|r2), REML=FALSE, data=dat) print(summary(rc, correlation = FALSE), show.resids = FALSE)The results of the above commands are shown below.
Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: yc ~ x1 + (1 | r1) + (1 | r2) Data: dat AIC BIC logLik deviance df.resid -3141.0 -3113.0 1575.5 -3151.0 1995 Random effects: Groups Name Variance Std.Dev. r2 (Intercept) 0.08729 0.2954 r1 (Intercept) 0.25677 0.5067 Residual 0.01036 0.1018 Number of obs: 2000, groups: r2, 25; r1, 20 Fixed effects: Estimate Std. Error t value (Intercept) 2.99384 0.12781 23.4 x1 -0.99821 0.00233 -428.5
The standard deviations used for r1 and r2 when generating yc were .3 and .6 respectively. The rc model calculate these standard deviations as 0.5067 and 0.2954. The rc model does not provide standard errors for either the estimated variance or standard deviation for r1 and r2. So there is no immediate means to determine if the true effect is likely to be greater than zero. A p-value could be calculated by use of a parametric bootstrap.
Nested random effect example
Nested effects, like crossed effects, are specified as separate random terms in the model formula. Nested effects are assumed to have an association between them. Each effect level of the lower nesting level is associated with only one effect level of the next higher nesting level. It is this assignment of the group identifiers that creates the association between the levels, not the formula specification. The following code builds the nested model fo yn.
Enter the following commands in your script and run them.
rn <- lmer(yn ~ x1 + (1|r1) + (1|r3), data=dat) print(summary(rn, correlation = FALSE), show.resids = FALSE)The results of the above commands are shown below.
Linear mixed model fit by REML ['lmerMod'] Formula: yn ~ x1 + (1 | r1) + (1 | r3) Data: dat REML criterion at convergence: -1680.5 Random effects: Groups Name Variance Std.Dev. r3 (Intercept) 0.08818 0.2970 r1 (Intercept) 0.26895 0.5186 Residual 0.01033 0.1016 Number of obs: 2000, groups: r3, 491; r1, 20 Fixed effects: Estimate Std. Error t value (Intercept) 3.003918 0.116765 25.7 x1 -0.997991 0.002621 -380.7
The rn model calculate these standard deviations of r1 and r3 as 0.5186 and 0.297. As with the crossed model, a parametric bootstrap can be used to generate a p-value for these standard deviations.
The group variable r3 could have also been specified in the formula as r1:r3. The following code generates the model for yn with this alternate coding of the nested effect.
Enter the following commands in your script and run them.
rn3 <- lmer(yn ~ x1 + (1|r1) + (1|r1:r3), data=dat) print(summary(rn3, correlation = FALSE), show.resids = FALSE)The results of the above commands are shown below.
Linear mixed model fit by REML ['lmerMod'] Formula: yn ~ x1 + (1 | r1) + (1 | r1:r3) Data: dat REML criterion at convergence: -1680.5 Random effects: Groups Name Variance Std.Dev. r1:r3 (Intercept) 0.08818 0.2970 r1 (Intercept) 0.26895 0.5186 Residual 0.01033 0.1016 Number of obs: 2000, groups: r1:r3, 491; r1, 20 Fixed effects: Estimate Std. Error t value (Intercept) 3.003918 0.116765 25.7 x1 -0.997991 0.002621 -380.7
The models have the same estimates for the effects. The only difference between the rn3 model and the rn model is the name of the grouping variable used for the nested effect.
Specifying the random term as r1:r3 has the advantage of making it clear to the readers of your code that the r3 effects are nested within the r1 effects. The group variable r3 could have also been specified in the formula as r1:r2. The r1:r2 coding is how the r3 variable was generated. This coding of the r3 groups would be used whenever the lower nesting level groups identify a group relative to the higher nesting level effect. For example in the shipping pallet example where each pallet has boxes of the same product, there would be a first box for each shipping pallet. Each of these first boxes are a different box and there is no relationship between the first boxes on any two different shipping pallets. The box grouping variable would need to be coded as interacted with the pallet grouping variable.
Mis-specifying a crossed or nested effect
There are times when it is not clear if a random effect is nested or crossed. If the variables' true relationship is crossed and the relationship is coded as nested, extra groups will be created for the lower nesting level and may result in an over-fit model. If the variables' true relationship is nested and the relationship is coded as crossed, there will be fewer groups than should be created for the lower nesting level. To see what this might do to a model, we will fit the nested model with the random variables crossed. The following code uses the yn response with the random terms 1|r1 and 1|r2 instead of 1|r1 and 1|r1:r2.
Enter the following commands in your script and run them.
rnc <- lmer(yn ~ x1 + (1|r1) + (1|r2), data=dat) print(summary(rnc, correlation = FALSE), show.resids = FALSE)The results of the above commands are shown below.
Linear mixed model fit by REML ['lmerMod'] Formula: yn ~ x1 + (1 | r1) + (1 | r2) Data: dat REML criterion at convergence: 939 Random effects: Groups Name Variance Std.Dev. r2 (Intercept) 0.004309 0.06564 r1 (Intercept) 0.264841 0.51463 Residual 0.086370 0.29389 Number of obs: 2000, groups: r2, 25; r1, 20 Fixed effects: Estimate Std. Error t value (Intercept) 3.012074 0.116010 25.96 x1 -0.991309 0.006718 -147.56
This model's estimate for r1 is almost the same. The variance associated with the other random variable is now much smaller and the residual variance has increase by about 2.5. Since we generated the data we know that the residual variance should be .01 (.1 squared). The nested model (rn) has the correct residual variance. In the model assuming crossed terms (rnc) some of the variance of yn could not be identified as being associated with r2 and instead the variance was included with the residuals.
The following code tests the significance of the r2 random variable in the rnc model.
Testing the significance of the random effects
Enter the following commands in your script and run them.
rnNr2 <- lmer(yn ~ x1 + (1|r1), data=dat) pbnmrnNr2 <- pbnm(rnc,rnNr2,nsim=1000,tasks=10,cores=2,seed=570445) summary(pbnmrnNr2)The results of the above commands are shown below.
Parametric bootstrap testing: (Intercept) | r2 = 0 from: lmer(formula = yn ~ x1 + (1 | r1) + (1 | r2), data = dat) 1000 samples were taken Thu Aug 18 14:50:39 2016 P(abs((Intercept) | r2) > |0.00430897|) = 0
The r2 random variable is significant even though it is much smaller than its true value. There is no indication in the model that there is an error in the specification of the random terms. It is the science of the project which is needed to determine what the correct relationship for the random model terms.
Previous: Diagnostics and Other Inferences
Last Revised: 6/24/2016