In the previous two posts, I applied custom hypothesis matrix (and corresponding contrast matrix) in linear models with between-subjects designs (\(2\times2\) design and \(2\times2\times2\) design). In those models, only the fixed effects can be estimated. In this post, I will try to apply priori contrasts in linear mixed models where both fixed and random effects can be estimated.
1 Preparation
library(tidyverse)
library(MASS)
library(hypr)
library(lme4)
library(optimx)
set.seed(20201114)Different from the previous two posts, a \(2\times3\) within-subjects design is used for simulating the dataset here. The full explanation of the simulation codes is available in this post. In short, for this dataset, thirty participants viewed 30 faces with different spatial frequency (low, medium, vs. high) in different orientations (upright vs. inverted), while the event-related potentials (ERP), i.e., the dependent variable, are recorded.
devtools::source_url("https://github.com/HaiyangJin/Website-Shared-Code/blob/master/simulate_data_lmm_v2.R?raw=TRUE")
# randomly set sd of two random effects as .001 (then we should get a singular fit warning for the maximal model)
P <- list(beta_ori_u_sd = .001,
beta_hm_u_sd = .001)
# use the default values for other parameters
simulation <- simudata_lmm2(P)
(df_lmm <- simulation$simudata)## # A tibble: 5,400 x 5
## subject stimulus Orientation SF erp
## <int> <int> <fct> <fct> <dbl>
## 1 1 1 upright low 1.08
## 2 1 1 upright medium -4.13
## 3 1 1 upright high -1.81
## 4 1 1 inverted low -3.90
## 5 1 1 inverted medium -5.15
## 6 1 1 inverted high -1.11
## 7 1 2 upright low -4.35
## 8 1 2 upright medium -4.94
## 9 1 2 upright high -4.67
## 10 1 2 inverted low -3.86
## # … with 5,390 more rows# true values for population parameters (used for simulating data)
(true_value <- simulation$population)## upr_low upr_med upr_hig inv_low inv_med inv_hig
## -1.750000 -3.250000 -4.000000 -4.583333 -5.083333 -5.333333# Note: the contrasts in this dataset was already set with contr.sdif()
contrasts(df_lmm$Orientation)## 2-1
## upright -0.5
## inverted 0.5contrasts(df_lmm$SF)## 2-1 3-2
## low -0.6666667 -0.3333333
## medium 0.3333333 -0.3333333
## high 0.3333333 0.66666672 Hypothesis and the contrast matrix
Suppose we are interested in the differences between upright and inverted conditions nested the three different spatial frequency levels, i.e., inv_low - upr_low, inv_med - upr_med, and inv_hig - upr_hig. These hypothese are constructed with library(hypr) as followings:
hypr_lmm <- hypr(inv_low - upr_low ~ 0, # 1
inv_med - upr_med ~ 0, # 2
inv_hig - upr_hig ~ 0, # 3
(inv_med + upr_med)/2 - (inv_low + upr_low)/2 ~ 0, # 4
(inv_hig + upr_hig)/2 - (inv_med + upr_med)/2 ~ 0, # 5
levels = names(simulation$population))
hypr_lmm## hypr object containing 5 null hypotheses:
## H0.1: 0 = inv_low - upr_low
## H0.2: 0 = inv_med - upr_med
## H0.3: 0 = inv_hig - upr_hig
## H0.4: 0 = 1/2*inv_med + 1/2*upr_med - 1/2*inv_low - 1/2*upr_low
## H0.5: 0 = 1/2*inv_hig + 1/2*upr_hig - 1/2*inv_med - 1/2*upr_med
##
## Hypothesis matrix (transposed):
## [,1] [,2] [,3] [,4] [,5]
## upr_low -1 0 0 -1/2 0
## upr_med 0 -1 0 1/2 -1/2
## upr_hig 0 0 -1 0 1/2
## inv_low 1 0 0 -1/2 0
## inv_med 0 1 0 1/2 -1/2
## inv_hig 0 0 1 0 1/2
##
## Contrast matrix:
## [,1] [,2] [,3] [,4] [,5]
## upr_low -1/2 0 0 -2/3 -1/3
## upr_med 0 -1/2 0 1/3 -1/3
## upr_hig 0 0 -1/2 1/3 2/3
## inv_low 1/2 0 0 -2/3 -1/3
## inv_med 0 1/2 0 1/3 -1/3
## inv_hig 0 0 1/2 1/3 2/3Then a new variable is created by combining Orientation and SF (spatial frequency):
df_lmm <- df_lmm %>%
mutate(I_ = paste(substr(Orientation, 1, 3), substr(SF, 1, 3), sep = "_"),
I_ = factor(I_, levels = names(simulation$population)))
# apply the contrast to I_
contrasts(df_lmm$I_) <- contr.hypothesis(hypr_lmm)Convert the contrast matrix to design matrix:
(df_var <- as_tibble(model.matrix(~ 1 + I_, df_lmm)))## # A tibble: 5,400 x 6
## `(Intercept)` I_1 I_2 I_3 I_4 I_5
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 -0.5 0 0 -0.667 -0.333
## 2 1 0 -0.5 0 0.333 -0.333
## 3 1 0 0 -0.5 0.333 0.667
## 4 1 0.5 0 0 -0.667 -0.333
## 5 1 0 0.5 0 0.333 -0.333
## 6 1 0 0 0.5 0.333 0.667
## 7 1 -0.5 0 0 -0.667 -0.333
## 8 1 0 -0.5 0 0.333 -0.333
## 9 1 0 0 -0.5 0.333 0.667
## 10 1 0.5 0 0 -0.667 -0.333
## # … with 5,390 more rows# combine df_var and df_lmm
df_lmm2 <- cbind(df_lmm, df_var)3 Fit linear mixed models
3.1 The maximal model with I_
if (!file.exists("2020-11-14-priori-contrast-lmm-maxi.rds")){
lmm_fit <- lmer(erp ~ I_ +
(1 + I_ | subject) +
(1 + I_ | stimulus),
data = df_lmm,
control = lmerControl(optimizer = "optimx", # calc.derivs = FALSE,
optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
save("lmm_fit", file = "2020-11-14-priori-contrast-lmm-maxi.rds")
} else {
load("2020-11-14-priori-contrast-lmm-maxi.rds")
}
summary(lmm_fit, corr = FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: erp ~ I_ + (1 + I_ | subject) + (1 + I_ | stimulus)
## Data: df_lmm
## Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb",
## starttests = FALSE, kkt = FALSE))
##
## REML criterion at convergence: 23108.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7943 -0.6498 0.0183 0.6703 3.3004
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 5.63000 2.3728
## I_1 0.24207 0.4920 -0.80
## I_2 0.08074 0.2841 0.07 -0.06
## I_3 0.42266 0.6501 0.78 -0.88 -0.40
## I_4 0.13773 0.3711 0.17 -0.55 0.06 0.37
## I_5 0.01974 0.1405 0.23 -0.01 -0.72 0.42 -0.59
## stimulus (Intercept) 1.31182 1.1453
## I_1 0.24454 0.4945 0.17
## I_2 0.20913 0.4573 0.24 0.54
## I_3 0.19300 0.4393 0.45 0.46 0.19
## I_4 0.06782 0.2604 0.13 0.36 0.30 -0.08
## I_5 0.15969 0.3996 0.21 0.33 0.37 0.86 0.07
## Residual 3.90761 1.9768
## Number of obs: 5400, groups: subject, 30; stimulus, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -4.3970 0.4818 -9.127
## I_1 -2.7476 0.1578 -17.411
## I_2 -1.6917 0.1354 -12.490
## I_3 -1.2640 0.1709 -7.396
## I_4 -0.9945 0.1058 -9.400
## I_5 -0.6048 0.1016 -5.953
## optimizer (optimx) convergence code: 1
## boundary (singular) fit: see ?isSingularWe can see that the singualr fit warning shows up in the above results. Thus, we have to remove some random effects from the maximal model. The whole process of removing random effects is following Bates et al., (2015).
3.2 The maximal model with design matrix
Before applying the steps to remove random effects, I re-fit the maximal model but with the design matrix (with each contrast as one preditive variable, e.g., I_1, I_2 etc) instead of the single variable I_.
if (!file.exists("2020-11-14-priori-contrast-lmm-max.rds")){
lmm_max <- lmer(erp ~ I_1 + I_2 + I_3 + I_4 + I_5 +
(1 + I_1 + I_2 + I_3 + I_4 + I_5 | subject) +
(1 + I_1 + I_2 + I_3 + I_4 + I_5 | stimulus),
data = df_lmm2,
control = lmerControl(optimizer = "optimx", # calc.derivs = FALSE,
optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
save("lmm_max", file = "2020-11-14-priori-contrast-lmm-max.rds")
} else {
load("2020-11-14-priori-contrast-lmm-max.rds")
}
summary(lmm_max, corr = FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: erp ~ I_1 + I_2 + I_3 + I_4 + I_5 + (1 + I_1 + I_2 + I_3 + I_4 +
## I_5 | subject) + (1 + I_1 + I_2 + I_3 + I_4 + I_5 | stimulus)
## Data: df_lmm2
## Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb",
## starttests = FALSE, kkt = FALSE))
##
## REML criterion at convergence: 23108.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7943 -0.6498 0.0183 0.6703 3.3004
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 5.63000 2.3728
## I_1 0.24207 0.4920 -0.80
## I_2 0.08074 0.2841 0.07 -0.06
## I_3 0.42266 0.6501 0.78 -0.88 -0.40
## I_4 0.13773 0.3711 0.17 -0.55 0.06 0.37
## I_5 0.01974 0.1405 0.23 -0.01 -0.72 0.42 -0.59
## stimulus (Intercept) 1.31182 1.1453
## I_1 0.24454 0.4945 0.17
## I_2 0.20913 0.4573 0.24 0.54
## I_3 0.19300 0.4393 0.45 0.46 0.19
## I_4 0.06782 0.2604 0.13 0.36 0.30 -0.08
## I_5 0.15969 0.3996 0.21 0.33 0.37 0.86 0.07
## Residual 3.90761 1.9768
## Number of obs: 5400, groups: subject, 30; stimulus, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -4.3970 0.4818 -9.127
## I_1 -2.7476 0.1578 -17.411
## I_2 -1.6917 0.1354 -12.490
## I_3 -1.2640 0.1709 -7.396
## I_4 -0.9945 0.1058 -9.400
## I_5 -0.6048 0.1016 -5.953
## optimizer (optimx) convergence code: 1
## boundary (singular) fit: see ?isSingularThe results of lmm_fit (fitted with I_) and lmm_max (fitted with I_1, I_2, etc) are the same. The main point here is that when “different” formula (e.g., lmm_fit and lmm_max) are used to fit the same model, same results could be obtained. Thus, either formula can be used. However, it is more convenient/easier to remove random effects from the formula used in lmm_max.
3.3 The zero-correlation-parameter model
At first, correlations between random effects are removed, making the zero-correlation-parameter (zcp) model (lmm_zcp):
if (!file.exists("2020-11-14-priori-contrast-lmm-zcp.rds")){
lmm_zcp <- lmer(erp ~ I_1 + I_2 + I_3 + I_4 + I_5 +
(1 + I_1 + I_2 + I_3 + I_4 + I_5 || subject) +
(1 + I_1 + I_2 + I_3 + I_4 + I_5 || stimulus),
data = df_lmm2,
control = lmerControl(optimizer = "optimx", # calc.derivs = FALSE,
optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
save("lmm_zcp", file = "2020-11-14-priori-contrast-lmm-zcp.rds")
} else {
load("2020-11-14-priori-contrast-lmm-zcp.rds")
}
summary(lmm_zcp, corr = FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: erp ~ I_1 + I_2 + I_3 + I_4 + I_5 + ((1 | subject) + (0 + I_1 |
## subject) + (0 + I_2 | subject) + (0 + I_3 | subject) + (0 +
## I_4 | subject) + (0 + I_5 | subject)) + ((1 | stimulus) +
## (0 + I_1 | stimulus) + (0 + I_2 | stimulus) + (0 + I_3 |
## stimulus) + (0 + I_4 | stimulus) + (0 + I_5 | stimulus))
## Data: df_lmm2
## Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb",
## starttests = FALSE, kkt = FALSE))
##
## REML criterion at convergence: 23161
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7912 -0.6389 0.0140 0.6653 3.3590
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 5.63073 2.3729
## subject.1 I_1 0.21642 0.4652
## subject.2 I_2 0.05023 0.2241
## subject.3 I_3 0.39015 0.6246
## subject.4 I_4 0.10820 0.3289
## subject.5 I_5 0.00000 0.0000
## stimulus (Intercept) 1.31099 1.1450
## stimulus.1 I_1 0.24396 0.4939
## stimulus.2 I_2 0.20841 0.4565
## stimulus.3 I_3 0.19188 0.4380
## stimulus.4 I_4 0.07075 0.2660
## stimulus.5 I_5 0.16393 0.4049
## Residual 3.91881 1.9796
## Number of obs: 5400, groups: subject, 30; stimulus, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -4.39704 0.48178 -9.127
## I_1 -2.74761 0.15510 -17.716
## I_2 -1.69168 0.13164 -12.851
## I_3 -1.26403 0.16766 -7.539
## I_4 -0.99446 0.10158 -9.790
## I_5 -0.60481 0.09909 -6.104
## optimizer (optimx) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular3.4 The reduced model
rePCA() is used to check the random effects.
summary(rePCA(lmm_zcp))## $subject
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6]
## Standard deviation 1.1987 0.3155 0.23500 0.16616 0.11321 0
## Proportion of Variance 0.8804 0.0610 0.03384 0.01692 0.00785 0
## Cumulative Proportion 0.8804 0.9414 0.97523 0.99215 1.00000 1
##
## $stimulus
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6]
## Standard deviation 0.5784 0.2495 0.23061 0.22128 0.20453 0.13437
## Proportion of Variance 0.5986 0.1114 0.09517 0.08762 0.07486 0.03231
## Cumulative Proportion 0.5986 0.7100 0.80521 0.89284 0.96769 1.00000Two of by-subject random effects (I_5 and I_2) explained less than .1% variances and they were removed from the zcp model, making the reduced model (lmm_rdc).
if (!file.exists("2020-11-14-priori-contrast-lmm-rdc.rds")){
lmm_rdc <- lmer(erp ~ I_1 + I_2 + I_3 + I_4 + I_5 +
(1 + I_1 + I_3 + I_4 || subject) + # I_2 + + I_5
(1 + I_1 + I_2 + I_3 + I_4 + I_5 || stimulus),
data = df_lmm2,
control = lmerControl(optimizer = "optimx", # calc.derivs = FALSE,
optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
save("lmm_rdc", file = "2020-11-14-priori-contrast-lmm-rdc.rds")
} else {
load("2020-11-14-priori-contrast-lmm-rdc.rds")
}
summary(lmm_rdc, corr = FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: erp ~ I_1 + I_2 + I_3 + I_4 + I_5 + ((1 | subject) + (0 + I_1 |
## subject) + (0 + I_3 | subject) + (0 + I_4 | subject)) + ((1 |
## stimulus) + (0 + I_1 | stimulus) + (0 + I_2 | stimulus) +
## (0 + I_3 | stimulus) + (0 + I_4 | stimulus) + (0 + I_5 | stimulus))
## Data: df_lmm2
## Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb",
## starttests = FALSE, kkt = FALSE))
##
## REML criterion at convergence: 23161.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7893 -0.6409 0.0125 0.6666 3.3335
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 5.63046 2.3729
## subject.1 I_1 0.21614 0.4649
## subject.2 I_3 0.38986 0.6244
## subject.3 I_4 0.10809 0.3288
## stimulus (Intercept) 1.31097 1.1450
## stimulus.1 I_1 0.24368 0.4936
## stimulus.2 I_2 0.20813 0.4562
## stimulus.3 I_3 0.19160 0.4377
## stimulus.4 I_4 0.07065 0.2658
## stimulus.5 I_5 0.16383 0.4048
## Residual 3.92308 1.9807
## Number of obs: 5400, groups: subject, 30; stimulus, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -4.3970 0.4818 -9.127
## I_1 -2.7476 0.1551 -17.719
## I_2 -1.6917 0.1251 -13.520
## I_3 -1.2640 0.1676 -7.541
## I_4 -0.9945 0.1016 -9.791
## I_5 -0.6048 0.0991 -6.1033.5 The extended model
Then the correlations are added back to the zcp model, making the extended model (lmm_etd).
if (!file.exists("2020-11-14-priori-contrast-lmm-etd.rds")){
lmm_etd <- lmer(erp ~ I_1 + I_2 + I_3 + I_4 + I_5 +
(1 + I_1 + I_3 + I_4 | subject) + # I_2 + + I_5
(1 + I_1 + I_2 + I_3 + I_4 + I_5 | stimulus),
data = df_lmm2,
control = lmerControl(optimizer = "optimx", # calc.derivs = FALSE,
optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
save("lmm_etd", file = "2020-11-14-priori-contrast-lmm-etd.rds")
} else {
load("2020-11-14-priori-contrast-lmm-etd.rds")
}
summary(lmm_etd, corr = FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: erp ~ I_1 + I_2 + I_3 + I_4 + I_5 + (1 + I_1 + I_3 + I_4 | subject) +
## (1 + I_1 + I_2 + I_3 + I_4 + I_5 | stimulus)
## Data: df_lmm2
## Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb",
## starttests = FALSE, kkt = FALSE))
##
## REML criterion at convergence: 23114.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7817 -0.6540 0.0174 0.6684 3.2965
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 5.62949 2.3727
## I_1 0.22918 0.4787 -0.82
## I_3 0.39459 0.6282 0.81 -0.99
## I_4 0.10827 0.3290 0.24 -0.62 0.54
## stimulus (Intercept) 1.31099 1.1450
## I_1 0.24390 0.4939 0.17
## I_2 0.20814 0.4562 0.25 0.54
## I_3 0.19168 0.4378 0.45 0.46 0.19
## I_4 0.06741 0.2596 0.14 0.36 0.30 -0.08
## I_5 0.15908 0.3988 0.21 0.33 0.37 0.86 0.08
## Residual 3.92190 1.9804
## Number of obs: 5400, groups: subject, 30; stimulus, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -4.39704 0.48174 -9.127
## I_1 -2.74761 0.15648 -17.559
## I_2 -1.69168 0.12511 -13.521
## I_3 -1.26403 0.16810 -7.519
## I_4 -0.99446 0.10106 -9.840
## I_5 -0.60481 0.09829 -6.154
## optimizer (optimx) convergence code: 1
## Model failed to converge with max|grad| = 0.0144438 (tol = 0.002, component 1)summary(rePCA(lmm_etd))## $subject
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 1.2435 0.25368 0.10706 0.001734
## Proportion of Variance 0.9533 0.03967 0.00707 0.000000
## Cumulative Proportion 0.9533 0.99293 1.00000 1.000000
##
## $stimulus
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6]
## Standard deviation 0.5977 0.3202 0.23152 0.16791 0.11889 0.02060
## Proportion of Variance 0.6424 0.1843 0.09638 0.05069 0.02541 0.00076
## Cumulative Proportion 0.6424 0.8267 0.92313 0.97382 0.99924 1.00000One by-subject random-effect (I_1) and one by-stimulus random-effect (I_4) were removed from lmm_etd, making the updated extended model (lmm_etd1).
if (!file.exists("2020-11-14-priori-contrast-lmm-etd1.rds")){
lmm_etd1 <- lmer(erp ~ I_1 + I_2 + I_3 + I_4 + I_5 +
(1 + I_3 + I_4 | subject) + # I_2 + I_5 + I_1 +
(1 + I_1 + I_2 + I_3 + I_5 | stimulus), # + I_4
data = df_lmm2,
control = lmerControl(optimizer = "optimx", # calc.derivs = FALSE,
optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
save("lmm_etd1", file = "2020-11-14-priori-contrast-lmm-etd1.rds")
} else {
load("2020-11-14-priori-contrast-lmm-etd1.rds")
}
summary(lmm_etd1, corr = FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: erp ~ I_1 + I_2 + I_3 + I_4 + I_5 + (1 + I_3 + I_4 | subject) +
## (1 + I_1 + I_2 + I_3 + I_5 | stimulus)
## Data: df_lmm2
## Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb",
## starttests = FALSE, kkt = FALSE))
##
## REML criterion at convergence: 23140.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7911 -0.6454 0.0076 0.6715 3.3044
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 5.6305 2.3729
## I_3 0.3879 0.6228 0.82
## I_4 0.1073 0.3276 0.24 0.55
## stimulus (Intercept) 1.3108 1.1449
## I_1 0.2417 0.4916 0.18
## I_2 0.2061 0.4540 0.25 0.54
## I_3 0.1896 0.4354 0.45 0.46 0.20
## I_5 0.1831 0.4279 0.24 0.42 0.44 0.78
## Residual 3.9530 1.9882
## Number of obs: 5400, groups: subject, 30; stimulus, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -4.39704 0.48178 -9.127
## I_1 -2.74761 0.12977 -21.173
## I_2 -1.69168 0.12512 -13.520
## I_3 -1.26403 0.16743 -7.549
## I_4 -0.99446 0.08928 -11.139
## I_5 -0.60481 0.10244 -5.9043.6 The optimal model
lmm_etd1 was compared to lmm_rdc:
anova(lmm_etd1, lmm_rdc, refit = FALSE)## Data: df_lmm2
## Models:
## lmm_rdc: erp ~ I_1 + I_2 + I_3 + I_4 + I_5 + ((1 | subject) + (0 + I_1 | subject) + (0 + I_3 | subject) + (0 + I_4 | subject)) + ((1 | stimulus) + (0 + I_1 | stimulus) + (0 + I_2 | stimulus) + (0 + I_3 | stimulus) + (0 + I_4 | stimulus) + (0 + I_5 | stimulus))
## lmm_etd1: erp ~ I_1 + I_2 + I_3 + I_4 + I_5 + (1 + I_3 + I_4 | subject) + (1 + I_1 + I_2 + I_3 + I_5 | stimulus)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lmm_rdc 17 23196 23308 -11581 23162
## lmm_etd1 28 23197 23381 -11570 23141 20.886 11 0.03457 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1It seems that lmm_rdc was better than lmm_etd1. Thus, lmm_rdc is taken as the optimal model.
lmm_opt <- lmm_rdc
summary(lmm_opt)## Linear mixed model fit by REML ['lmerMod']
## Formula: erp ~ I_1 + I_2 + I_3 + I_4 + I_5 + ((1 | subject) + (0 + I_1 |
## subject) + (0 + I_3 | subject) + (0 + I_4 | subject)) + ((1 |
## stimulus) + (0 + I_1 | stimulus) + (0 + I_2 | stimulus) +
## (0 + I_3 | stimulus) + (0 + I_4 | stimulus) + (0 + I_5 | stimulus))
## Data: df_lmm2
## Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb",
## starttests = FALSE, kkt = FALSE))
##
## REML criterion at convergence: 23161.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7893 -0.6409 0.0125 0.6666 3.3335
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 5.63046 2.3729
## subject.1 I_1 0.21614 0.4649
## subject.2 I_3 0.38986 0.6244
## subject.3 I_4 0.10809 0.3288
## stimulus (Intercept) 1.31097 1.1450
## stimulus.1 I_1 0.24368 0.4936
## stimulus.2 I_2 0.20813 0.4562
## stimulus.3 I_3 0.19160 0.4377
## stimulus.4 I_4 0.07065 0.2658
## stimulus.5 I_5 0.16383 0.4048
## Residual 3.92308 1.9807
## Number of obs: 5400, groups: subject, 30; stimulus, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -4.3970 0.4818 -9.127
## I_1 -2.7476 0.1551 -17.719
## I_2 -1.6917 0.1251 -13.520
## I_3 -1.2640 0.1676 -7.541
## I_4 -0.9945 0.1016 -9.791
## I_5 -0.6048 0.0991 -6.103
##
## Correlation of Fixed Effects:
## (Intr) I_1 I_2 I_3 I_4
## I_1 0.000
## I_2 0.000 0.000
## I_3 0.000 0.000 0.000
## I_4 0.000 0.000 0.000 0.000
## I_5 0.000 0.000 0.000 0.000 -0.2173.7 Comparisons between estimates and “the true values”
| coefficient | calculation | true value |
estimate |
|---|---|---|---|
(Intercept) |
(inv_low + upr_low + inv_med + upr_med + inv_hig + upr_hig)/6 |
-4 | -4.4 |
I_1 |
inv_low - upr_low |
-2.83 | -2.75 |
I_2 |
inv_med - upr_med |
-1.83 | -1.69 |
I_3 |
inv_hig - upr_hig |
-1.33 | -1.26 |
I_4 |
(inv_med + upr_med)/2 - (inv_low + upr_low)/2 |
-1 | -0.99 |
I_5 |
(inv_hig + upr_hig)/2 - (inv_med + upr_med)/2 |
-0.5 | -0.6 |
The estimates of lmm_opt are simliar to the “true values”, which were used to simulate the data. For the hypothesis of interests, we can check the estimates of I_1, I_2, and I_3.
4 Summary
In this post, I applied custom contrasts in linear mixed models with a \(2\times3\) design. It is similar to the previous two posts for fitting the fixed effects. The main difference is that, to remove random effects in linear mixed models, the design matrix needs to be converted from the contrast matrix.
Finally, after completing these three posts, I feel more comfortable to apply priori contrasts in linear models (and linear mixed models). It will be great if these posts are also helpful for you.
References
Bates, D., Kliegl, R., Vasishth, S., & Baayen, R. H. (2015). Parsimonious Mixed Models. Retrieved from http://arxiv.org/abs/1506.04967