# lme4/Matrix: Call to .Call("mer_update_y"...) and LMEoptimize gives unexpected side effect...

2 messages
Open this post in threaded view
|
Report Content as Inappropriate

## lme4/Matrix: Call to .Call("mer_update_y"...) and LMEoptimize gives unexpected side effect...

 Dear all I want to compute Monte Carlo p-values in lmer-models based on sampled data sets. To speed up calculations, I've tried to use internal functions from the Matrix package (as suggested ealier on the list by Doug Bates). So I did:  fm2 <- lmer(resistance ~ ET + position + (1|Grp), Semiconductor,method='ML')  simdata<-simulate(fm2,nsim=1)  ynew <- simdata[,1]  mer <- fm2  .Call("mer_update_y", mer, ynew, PACKAGE = "Matrix")  mer1u <- LMEoptimize(mer, lmerControl(mer)) What puzzles me is that this call alters my original model fm2 as some kind of side effect. In fact, after the call fm2 is the same as mer1u. Is this side effect intentional and is it possible to avoid? A detail is that "LMEoptmize" and "LMEoptimize<-" are not exported from the namespace in Matrix, so I simply copied the LMEoptimize function and made it an ordinary function as shown below. Thanks in advance Søren LMEoptimize <- function(x, value)              {                  if (value\$msMaxIter < 1) return(x)                  nc <- x@nc                  constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))                  fn <- function(pars)                      deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))                  gr <- if (value\$gradient)                      function(pars) {                          if (!isTRUE(all.equal(pars,                                                .Call("mer_coef", x,                                                      2, PACKAGE = "Matrix"))))                              .Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")                          .Call("mer_gradient", x, 2, PACKAGE = "Matrix")                      }                  else NULL                  optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"),                                     fn, gr,                                     lower = ifelse(constr, 5e-10, -Inf),                                     control = list(iter.max = value\$msMaxIter,                                     trace = as.integer(value\$msVerbose)))                  .Call("mer_coefGets", x, optimRes\$par, 2, PACKAGE = "Matrix")                  if (optimRes\$convergence != 0) {                      warning(paste("nlminb returned message",                                    optimRes\$message,"\n"))                  }                  return(x)              } lmerControl <-   function(maxIter = 200, # used in ../src/lmer.c only            tolerance = sqrt(.Machine\$double.eps),# ditto            msMaxIter = 200,            ## msTol = sqrt(.Machine\$double.eps),            ## FIXME:  should be able to pass tolerances to nlminb()            msVerbose = getOption("verbose"),            niterEM = 15,            EMverbose = getOption("verbose"),            PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead            gradient = TRUE,            Hessian = FALSE # unused _FIXME_            ) {     list(maxIter = as.integer(maxIter),          tolerance = as.double(tolerance),          msMaxIter = as.integer(msMaxIter),          ## msTol = as.double(msTol),          msVerbose = as.integer(msVerbose),# "integer" on purpose          niterEM = as.integer(niterEM),          EMverbose = as.logical(EMverbose),          PQLmaxIt = as.integer(PQLmaxIt),          gradient = as.logical(gradient),          Hessian = as.logical(Hessian)) } ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide! http://www.R-project.org/posting-guide.html