knitr::opts_chunk$set(echo = TRUE)
suppressMessages(library("lme4"))
library("RePsychLing")
suppressMessages(library("tidyverse"))
data("Machines", package = "MEMSS")
Machines$Measure <- factor(rep(1:3))
Machines <- Machines[do.call(order, Machines), ] # reorder by Worker, then Machine
d <- as.tibble(Machines) For each of the contrasts (treatment, sliding difference, sum) I compute the contrast at the level of observations. We rename Machine to a factor VC for variance components. Then, I use these derived scores as input for the LMM with the (0 + VC | Worker) specification in FE and RE parts of the forumula to estimate variance and correlation parameters for the contrasts. The results should roughly agree with the results using the original observations in LMMs (1 + VC | Worker) specification.
# Levels and contrasts
d_trt <- d %>%
spread(Machine, score) %>%
mutate(dBA = B-A, dCA = C-A) %>%
select(Worker, Measure, A, dBA, dCA) %>%
gather(key=VC, score, -c(Worker, Measure), factor_key=TRUE)
d_trt## # A tibble: 54 x 4
## Worker Measure VC score
## <fct> <fct> <fct> <dbl>
## 1 1 1 A 52
## 2 1 2 A 52.8
## 3 1 3 A 53.1
## 4 2 1 A 51.8
## 5 2 2 A 52.8
## 6 2 3 A 53.1
## 7 3 1 A 60
## 8 3 2 A 60.2
## 9 3 3 A 58.4
## 10 4 1 A 51.1
## # ... with 44 more rowsprint(summary(m1b <- lmer(score ~ 0 + VC + (0 + VC | Worker), d_trt, REML=FALSE )),corr=FALSE)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + VC + (0 + VC | Worker)
## Data: d_trt
##
## AIC BIC logLik deviance df.resid
## 265.6 285.4 -122.8 245.6 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.30120 -0.40158 0.00849 0.47205 1.81511
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker VCA 13.432 3.665
## VCdBA 28.610 5.349 0.48
## VCdCA 11.167 3.342 -0.39 0.31
## Residual 2.077 1.441
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## VCA 52.356 1.534 34.124
## VCdBA 7.967 2.210 3.605
## VCdCA 13.917 1.406 9.899contrasts(d$Machine) <- contr.treatment(3)
print(summary(m1 <- lmer(score ~ 1 + Machine + (1 + Machine | Worker), d, REML=FALSE )),corr=FALSE)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.8158 3.7170
## Machine2 28.6862 5.3559 0.49
## Machine3 11.2431 3.3531 -0.36 0.30
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.534 34.124
## Machine2 7.967 2.210 3.605
## Machine3 13.917 1.406 9.899d_sld <- d %>%
spread(Machine, score) %>%
mutate(GM = (A+B+C)/3, dBA = B-A, dCB = C-B) %>%
select(Worker, Measure, GM, dBA, dCB) %>%
gather(key=VC, score, -c(Worker, Measure))
d_sld## # A tibble: 54 x 4
## Worker Measure VC score
## <fct> <fct> <chr> <dbl>
## 1 1 1 GM 60.5
## 2 1 2 GM 60.9
## 3 1 3 GM 61.3
## 4 2 1 GM 57.7
## 5 2 2 GM 58.2
## 6 2 3 GM 58.1
## 7 3 1 GM 66.5
## 8 3 2 GM 65.5
## 9 3 3 GM 66.4
## 10 4 1 GM 59.5
## # ... with 44 more rowsprint(summary(m1b_sld <- lmer(score ~ 0 + VC + (0 + VC | Worker), d_sld, REML=FALSE )),corr=FALSE)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + VC + (0 + VC | Worker)
## Data: d_sld
##
## AIC BIC logLik deviance df.resid
## 255.4 275.2 -117.7 235.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.55157 -0.36112 -0.02514 0.50343 2.00506
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker VCdBA 28.781 5.365
## VCdCB 29.405 5.423 -0.81
## VCGM 22.476 4.741 0.83 -0.78
## Residual 1.565 1.251
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## VCdBA 7.967 2.210 3.605
## VCdCB 5.950 2.233 2.664
## VCGM 59.650 1.958 30.468contrasts(d$Machine) <- MASS::contr.sdif(3)
print(summary(m1_sld <- lmer(score ~ 1 + Machine + (1 + Machine | Worker), d, REML=FALSE )),corr=FALSE)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 22.8953 4.7849
## Machine2-1 28.6862 5.3559 0.82
## Machine3-2 29.3098 5.4139 -0.77 -0.81
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.650 1.958 30.468
## Machine2-1 7.967 2.210 3.605
## Machine3-2 5.950 2.233 2.664# Note different ordering of VPs/CPs in outputd_sum <- d %>%
spread(Machine, score) %>%
mutate(GM = (A+B+C)/3, dA = A-GM, dB = B-GM) %>%
select(Worker, Measure, GM, dA, dB) %>%
gather(key=VC, score, -c(Worker, Measure))
d_sum## # A tibble: 54 x 4
## Worker Measure VC score
## <fct> <fct> <chr> <dbl>
## 1 1 1 GM 60.5
## 2 1 2 GM 60.9
## 3 1 3 GM 61.3
## 4 2 1 GM 57.7
## 5 2 2 GM 58.2
## 6 2 3 GM 58.1
## 7 3 1 GM 66.5
## 8 3 2 GM 65.5
## 9 3 3 GM 66.4
## 10 4 1 GM 59.5
## # ... with 44 more rowsprint(summary(m1b_sum <- lmer(score ~ 0 + VC + (0 + VC | Worker), d_sum, REML=FALSE )),corr=FALSE)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + VC + (0 + VC | Worker)
## Data: d_sum
##
## AIC BIC logLik deviance df.resid
## 210.5 230.4 -95.3 190.5 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.16841 -0.47852 -0.02479 0.43231 2.55509
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker VCdA 5.6055 2.368
## VCdB 11.6277 3.410 -0.72
## VCGM 22.7815 4.773 -0.65 0.84
## Residual 0.6496 0.806
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## VCdA -7.2944 0.9851 -7.405
## VCdB 0.6722 1.4050 0.478
## VCGM 59.6500 1.9578 30.468contrasts(d$Machine) <- contr.sum(3)
print(summary(m1_sum <- lmer(score ~ 1 + Machine + (1 + Machine | Worker), d, REML=FALSE )),corr=FALSE)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 22.8953 4.7849
## Machine1 5.6165 2.3699 -0.65
## Machine2 11.6388 3.4116 0.84 -0.71
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 1.9578 30.468
## Machine1 -7.2944 0.9851 -7.405
## Machine2 0.6722 1.4050 0.478# Note different ordering of VPs/CPs in outputComputing contrast vectors at the level of the observations outside the LMM and using these derived measures as levels of a factor VC in an LMM with the (0 + VC | Worker) specification in FE and RE yields (almost) equal results compared to using the original observations in an LMM with (1 + Machine | Worker) specification. There are small difference in VPs and CPs and somewhat larger differences in the residual variance of the VC-factor LMMs. I expected them to be smaller for the latter due to the shrinkage vs. within-worker computation of differences This expectation was met for treatment and sliding-difference contrasts, but not for sum-based contrasts. The reason could be that in the sum contrast the difference scores use the grand mean (reducing noise) and one of the factor levels whereas in the other two contrasts the difference scores are between two factor levels.
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 forcats_0.3.0 stringr_1.3.1
## [4] dplyr_0.7.6 purrr_0.2.5 readr_1.1.1
## [7] tidyr_0.8.1 tibble_1.4.2 ggplot2_3.0.0
## [10] tidyverse_1.2.1 RePsychLing_0.0.4 lme4_1.1-18-1
## [13] Matrix_1.2-14
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.4 splines_3.5.1 haven_1.1.2
## [4] lattice_0.20-35 colorspace_1.3-2 htmltools_0.3.6
## [7] yaml_2.2.0 utf8_1.1.4 rlang_0.2.2
## [10] nloptr_1.0.4 pillar_1.3.0.9000 withr_2.1.2
## [13] glue_1.3.0 modelr_0.1.2 readxl_1.1.0
## [16] bindr_0.1.1 plyr_1.8.4 munsell_0.5.0
## [19] gtable_0.2.0 cellranger_1.1.0 rvest_0.3.2
## [22] evaluate_0.11 knitr_1.20 fansi_0.3.0
## [25] broom_0.5.0 Rcpp_0.12.18 backports_1.1.2
## [28] scales_1.0.0 jsonlite_1.5 hms_0.4.2
## [31] digest_0.6.17 stringi_1.2.4 grid_3.5.1
## [34] cli_1.0.1 tools_3.5.1 magrittr_1.5
## [37] lazyeval_0.2.1 crayon_1.3.4 pkgconfig_2.0.2
## [40] MASS_7.3-50 xml2_1.2.0 lubridate_1.7.4
## [43] rstudioapi_0.7 assertthat_0.2.0 minqa_1.2.4
## [46] rmarkdown_1.10.11 httr_1.3.1 htmldeps_0.1.1
## [49] R6_2.2.2 nlme_3.1-137 compiler_3.5.1