[R-sig-ME] simulate lme function value

[Ben Bolker]

bbonit at ... <bbonit at ...> writes:

> 
>  Dear list i have a question :
> 
> how can i extract simulated value from simulate function in nlme ? (if it
is possible)
> thank You so much 
> sorry for noise
> Bonitta Gianluca 

  I think you may be out of luck -- that is, I don't think that
simulate.lme stores the simulated values.  If I run the example
in ?simulate.lme and examine the resulting object, I get:

str(orthSim)
List of 2
 $ null:List of 2
  ..$ ML  : num [1:1000, 1:2] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:1000] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "info" "logLik"
  ..$ REML: num [1:1000, 1:2] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:1000] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "info" "logLik"
 $ alt :List of 2
  ..$ ML  : num [1:1000, 1:2] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:1000] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "info" "logLik"
  ..$ REML: num [1:1000, 1:2] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:1000] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "info" "logLik"
 - attr(*, "call")= language simulate.lme(object = 
  list(fixed = distance ~ age, data = Orthodont,
   random = ~1 |      Subject), nsim = 1000, 
   m2 = list(random = ~age | Subject))
 - attr(*, "seed")= int 61390140
 - attr(*, "df")= num 2
 - attr(*, "useGen")= logi FALSE
 - attr(*, "class")= chr "simulate.lme"

In words -- the ML and REML criteria are stored for the reduced
and full models for each simulation, but not the values themselves.
You *might* be able to hack the function: the critical lines are

base2 <- base + rnorm(N, sd = sig)
        for (j in 1:Q) {
            base2 <- base2 + ((array(rnorm(ngrp[j] * qvec[j]), 
                c(ngrp[j], qvec[j]), list(1:ngrp[j], NULL)) %*% 
                DeltaInv[[j]])[ind[[j]], , drop = FALSE] * condL1$Xy[, 
                csq1[j]:csq2[j], drop = FALSE]) %*% rep(1, qvec[j])
        }