This is an R Markdown document showing an example for bootmlm. The goal is to do a case bootstrapping of the fixed effects in an ‘lmer model’.
library(mlmRev) # contains data## Loading required package: lme4## Loading required package: Matrixlibrary(lme4)
data<-mlmRev::Exam From the slides
require(lme4)
schoolaverage<-lmer(normexam ~ standLRT + schavg + (1 + standLRT | school), data=Exam, REML=F)
summary(schoolaverage)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: normexam ~ standLRT + schavg + (1 + standLRT | school)
## Data: Exam
##
## AIC BIC logLik deviance df.resid
## 9324.4 9368.6 -4655.2 9310.4 4052
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8279 -0.6313 0.0337 0.6844 3.4370
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## school (Intercept) 0.07448 0.2729
## standLRT 0.01489 0.1220 0.38
## Residual 0.55362 0.7441
## Number of obs: 4059, groups: school, 65
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.001239 0.036675 -0.034
## standLRT 0.552392 0.020178 27.375
## schavg 0.294771 0.105590 2.792
##
## Correlation of Fixed Effects:
## (Intr) stnLRT
## standLRT 0.268
## schavg 0.089 -0.087Helper function to get the fixed effects.
mySumm <- function(x) {
c(getME(x, "beta"))
}Install the bootmlm package first from github, as described here
require(bootmlm)## Loading required package: bootmlmset.seed(1234)
boot_data <- bootstrap_mer(schoolaverage, mySumm, type = "case", nsim = 1000)## There are no names but corresponds to the fixed effects in the ‘schoolaverage’ model
require(boot)## Loading required package: bootrequire(skimr)## Loading required package: skimrskim(as.data.frame(boot_data$t)) | Name | as.data.frame(boot_data$t… |
| Number of rows | 1000 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| V1 | 0 | 1 | 0.00 | 0.04 | -0.13 | -0.03 | 0.00 | 0.02 | 0.13 | ▁▃▇▃▁ |
| V2 | 0 | 1 | 0.55 | 0.02 | 0.48 | 0.54 | 0.55 | 0.57 | 0.62 | ▁▃▇▅▁ |
| V3 | 0 | 1 | 0.30 | 0.13 | -0.11 | 0.22 | 0.30 | 0.38 | 0.83 | ▁▆▇▂▁ |
There are no missings. Somewhat wider than the example on the slides but still does not overlap with 0.
boots<-as.data.frame(boot_data$t)
quantile(boots$V3, prob = c(0.025, 0.975))## 2.5% 97.5%
## 0.0444086 0.5426102The papers underpinning ‘bootmlm’ suggests using ‘residual’ bootstraps. There are 7 bootstraps with convergence issues. In our case this confidence interval is a bit narrower and thus more supportive of an effect than the case bootstrap.
set.seed(12345)
boot_resid <- bootstrap_mer(schoolaverage, mySumm, type = "residual", nsim = 1000)## 1 warning(s) in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, lbound = lower) : Model failed to converge with max|grad| = 0.00215902 (tol = 0.002, component 1)
## 1 warning(s) in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, lbound = lower) : Model failed to converge with max|grad| = 0.00222617 (tol = 0.002, component 1)
## 1 warning(s) in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, lbound = lower) : Model failed to converge with max|grad| = 0.00239867 (tol = 0.002, component 1)
## 1 warning(s) in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, lbound = lower) : Model failed to converge with max|grad| = 0.00261775 (tol = 0.002, component 1)
## 1 warning(s) in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, lbound = lower) : Model failed to converge with max|grad| = 0.0028295 (tol = 0.002, component 1)
## 1 warning(s) in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, lbound = lower) : Model failed to converge with max|grad| = 0.00374849 (tol = 0.002, component 1)
## 1 warning(s) in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, lbound = lower) : Model failed to converge with max|grad| = 0.00675521 (tol = 0.002, component 1)skim(as.data.frame(boot_resid$t)) | Name | as.data.frame(boot_resid$… |
| Number of rows | 1000 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| V1 | 0 | 1 | 0.00 | 0.04 | -0.14 | -0.03 | 0.00 | 0.03 | 0.11 | ▁▂▇▅▁ |
| V2 | 0 | 1 | 0.55 | 0.02 | 0.48 | 0.53 | 0.55 | 0.56 | 0.62 | ▁▃▇▃▁ |
| V3 | 0 | 1 | 0.29 | 0.11 | -0.10 | 0.22 | 0.29 | 0.37 | 0.61 | ▁▂▇▆▁ |
boots_resid<-as.data.frame(boot_resid$t)
quantile(boots_resid$V3, prob = c(0.025, 0.975))## 2.5% 97.5%
## 0.0713145 0.5017247Note this was run on my home desktop which is R 4.0.2. Note that I have used a 1,000 bootstraps but you might want to use more.
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] skimr_2.1.3 boot_1.3-25 bootmlm_0.0.1.1000 mlmRev_1.0-8
## [5] lme4_1.1-26 Matrix_1.2-18
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.7 highr_0.9 compiler_4.0.2
## [4] pillar_1.6.2 nloptr_1.2.2.2 base64enc_0.1-3
## [7] tools_4.0.2 digest_0.6.27 statmod_1.4.35
## [10] jsonlite_1.7.2 tibble_3.1.3 evaluate_0.14
## [13] lifecycle_1.0.0 nlme_3.1-149 lattice_0.20-41
## [16] pkgconfig_2.0.3 rlang_0.4.11 DBI_1.1.0
## [19] yaml_2.2.1 xfun_0.25 repr_1.1.3
## [22] stringr_1.4.0 dplyr_1.0.7 knitr_1.33
## [25] generics_0.1.0 vctrs_0.3.8 tidyselect_1.1.1
## [28] grid_4.0.2 glue_1.4.2 R6_2.5.0
## [31] fansi_0.5.0 rmarkdown_2.7 minqa_1.2.4
## [34] tidyr_1.1.3 blob_1.2.1 purrr_0.3.4
## [37] magrittr_2.0.1 htmltools_0.5.1.1 ellipsis_0.3.2
## [40] MASS_7.3-53 splines_4.0.2 assertthat_0.2.1
## [43] numDeriv_2016.8-1.1 utf8_1.2.2 stringi_1.7.3
## [46] crayon_1.4.1