R help - Mar 2006 - 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 at 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))
}

You provided the code for "LMEoptimize" and "lmerControl",
but your
data.frame Semiconductor was "not found".  Therefore, I could not 
replicate your problem.  However, some external functions use their 
inputs for temporary storage or even to return values.  This was smart 
programming decades ago when programmer time was much cheaper than 
memory.  Today, I regard it as a user hostile feature of some code.

	  hope this helps.
	  spencer graves

 > sessionInfo()
R version 2.2.1, 2005-12-20, i386-pc-mingw32

attached base packages:
[1] "methods"   "stats"     "graphics" 
"grDevices" "utils"     "datasets"
[7] "base"

other attached packages:
      lme4   lattice    Matrix
"0.995-2" "0.12-11" "0.995-5"
 >

S?ren H?jsgaard wrote:
> 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 at 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))
> }
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html