sample of exactRLRT

[IrisTeng]

exactRLRT(m0)$sample is a function instead of simulated test statistics. How should I fix it?

[fabian-s]

Sorry, can't reproduce this -- I get:

library(RLRsim)
library(lme4)
#> Loading required package: Matrix
m0 <- lmer(Reaction ~ I(Days-4.5) + (1|Subject), data = sleepstudy)
rlr <- exactRLRT(m0)
str(rlr$sample)
#>  num [1:10000] 4.101 0 0.256 0.7 0 ...
#>  - attr(*, "lambda")= num [1:10000] 0.913 0 0.195 0.336 0 ...

^{Created on 2019-04-23 by the reprex package (v0.2.1)}

Session info

devtools::session_info()
#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.5.3 (2019-03-11)
#>  os       Linux Mint 19.1             
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language en_GB                       
#>  collate  en_GB.UTF-8                 
#>  ctype    en_GB.UTF-8                 
#>  tz       Europe/Berlin               
#>  date     2019-04-23                  
#> 
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package     * version  date       lib source        
#>  assertthat    0.2.1    2019-03-21 [1] CRAN (R 3.5.3)
#>  backports     1.1.4    2019-04-10 [1] CRAN (R 3.5.3)
#>  boot          1.3-20   2017-07-30 [1] CRAN (R 3.5.3)
#>  callr         3.2.0    2019-03-15 [1] CRAN (R 3.5.3)
#>  cli           1.1.0    2019-03-19 [1] CRAN (R 3.5.3)
#>  crayon        1.3.4    2017-09-16 [1] CRAN (R 3.5.3)
#>  desc          1.2.0    2018-05-01 [1] CRAN (R 3.5.3)
#>  devtools      2.0.2    2019-04-08 [1] CRAN (R 3.5.3)
#>  digest        0.6.18   2018-10-10 [1] CRAN (R 3.5.3)
#>  evaluate      0.13     2019-02-12 [1] CRAN (R 3.5.3)
#>  fs            1.2.7    2019-03-19 [1] CRAN (R 3.5.3)
#>  glue          1.3.1    2019-03-12 [1] CRAN (R 3.5.3)
#>  highr         0.8      2019-03-20 [1] CRAN (R 3.5.3)
#>  htmltools     0.3.6    2017-04-28 [1] CRAN (R 3.5.3)
#>  knitr         1.22     2019-03-08 [1] CRAN (R 3.5.3)
#>  lattice       0.20-38  2018-11-04 [4] CRAN (R 3.5.1)
#>  lme4        * 1.1-21   2019-03-05 [1] CRAN (R 3.5.3)
#>  magrittr      1.5      2014-11-22 [1] CRAN (R 3.5.3)
#>  MASS          7.3-51.1 2018-11-01 [4] CRAN (R 3.5.1)
#>  Matrix      * 1.2-17   2019-03-22 [4] CRAN (R 3.5.3)
#>  memoise       1.1.0    2017-04-21 [1] CRAN (R 3.5.3)
#>  mgcv          1.8-28   2019-03-21 [4] CRAN (R 3.5.3)
#>  minqa         1.2.4    2014-10-09 [1] CRAN (R 3.5.3)
#>  nlme          3.1-139  2019-04-09 [1] CRAN (R 3.5.3)
#>  nloptr        1.2.1    2018-10-03 [1] CRAN (R 3.5.3)
#>  pkgbuild      1.0.3    2019-03-20 [1] CRAN (R 3.5.3)
#>  pkgload       1.0.2    2018-10-29 [1] CRAN (R 3.5.3)
#>  prettyunits   1.0.2    2015-07-13 [1] CRAN (R 3.5.3)
#>  processx      3.3.0    2019-03-10 [1] CRAN (R 3.5.3)
#>  ps            1.3.0    2018-12-21 [1] CRAN (R 3.5.3)
#>  R6            2.4.0    2019-02-14 [1] CRAN (R 3.5.3)
#>  Rcpp          1.0.1    2019-03-17 [1] CRAN (R 3.5.3)
#>  remotes       2.0.4    2019-04-10 [1] CRAN (R 3.5.3)
#>  rlang         0.3.4    2019-04-07 [1] CRAN (R 3.5.3)
#>  RLRsim      * 3.1-3    2016-11-04 [1] CRAN (R 3.5.3)
#>  rmarkdown     1.12     2019-03-14 [1] CRAN (R 3.5.3)
#>  rprojroot     1.3-2    2018-01-03 [1] CRAN (R 3.5.3)
#>  sessioninfo   1.1.1    2018-11-05 [1] CRAN (R 3.5.3)
#>  stringi       1.4.3    2019-03-12 [1] CRAN (R 3.5.3)
#>  stringr       1.4.0    2019-02-10 [1] CRAN (R 3.5.3)
#>  testthat      2.0.1    2018-10-13 [1] CRAN (R 3.5.3)
#>  usethis       1.5.0    2019-04-07 [1] CRAN (R 3.5.3)
#>  withr         2.1.2    2018-03-15 [1] CRAN (R 3.5.3)
#>  xfun          0.6      2019-04-02 [1] CRAN (R 3.5.3)
#>  yaml          2.2.0    2018-07-25 [1] CRAN (R 3.5.3)
#> 
#> [1] /home/lmmista-wap218/R/x86_64-pc-linux-gnu-library/3.5
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/library

[fabian-s]

Please include a reproducible example of the problem.

[IrisTeng]

yes!

library(RLRsim)
library(nlme)
library(lme4)
form <- as.formula(c("~env+(env|ID)"))
beta <- c(1, -1)
names(beta) <- c("(Intercept)", "env")
N_ind <- 10
N_obs <- 100
simdat <- data.frame(ID = factor(rep(1:N_ind, each = N_obs)),
env = rnorm(N_ind * N_obs, 0, 1))
V_err0 <- 1
vcov <- matrix(c(0, 0, 0, 0), 2, 2)
theta <- c(0, 0, 0)
names(theta) <- c("ID.(Intercept)", "ID.env.(Intercept)", "ID.env")
response <- simulate(form, newdata = simdat, family = gaussian,
newparams = list(theta = theta, beta = beta, sigma = sqrt(V_err0)))
simdat$resp <- as.vector(response[, 1])

m0 <- lmer(respenv+(1|ID), data = simdat)
mA <- update(m0, .. + (0 + env|ID))
mSlope <- update(mA, .~. - (1|ID))
rlrt_cr <- as.matrix(exactRLRT(mSlope, mA, m0, nsim = 1e5)$sample, ncol = 1)

Sometimes it works but sometimes exactRLRT(mSlope, mA, m0, nsim = 1e5)$sample is a function.

[fabian-s]

Thx for the example -- this is due to a shortcut I take in case when the observed test statistic is 0.
The p-value in that case is always 1, so there's no need to simulate anything.
If you do need a sample of the test statistic for such cases, I'm afraid you'll have to use the RLRsim function directly, sorry.

[IrisTeng]

Gotcha, thank you so much!