SCM

SCM Repository

[matrix] Diff of /pkg/R/lmer.R
ViewVC logotype

Diff of /pkg/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 864, Thu Aug 18 21:36:49 2005 UTC revision 878, Sun Aug 28 13:57:39 2005 UTC
# Line 929  Line 929 
929  rWishart <- function(n, df, invScal)  rWishart <- function(n, df, invScal)
930    .Call("Matrix_rWishart", n, df, invScal)    .Call("Matrix_rWishart", n, df, invScal)
931    
932    setMethod("simulate", signature(obj = "lmer"),
933              function(obj, nsamp = 1, verbose = FALSE, ...)
934          {
935              family <- obj@family
936              if (family$family != "gaussian" ||
937                  family$link != "identity")
938                  stop("simulation of generalized linear mixed models not implemented yet")
939              ## create the mean from the fixed effects
940              frm <- obj@frame
941              n <- nrow(frm)
942              fixed.form <- Matrix:::nobars(obj@call$formula)
943              if (inherits(fixed.form, "name")) # RHS is empty - use a constant
944                  fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
945              glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE, y = TRUE)
946              fxd <- drop(glm.fit$x %*% fixef(obj))
947              bars <- Matrix:::findbars(obj@call$formula[[3]])
948              random <-
949                  lapply(bars,
950                         function(x) list(model.matrix(eval(substitute(~term,
951                                                                       list(term=x[[2]]))),
952                                                       frm),
953                                          eval(substitute(as.factor(fac)[,drop = TRUE],
954                                                          list(fac = x[[3]])), frm)))
955              names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
956              ## re-order the random effects pairs if necessary
957              if (any(names(random) != names(obj@flist)))
958                  random <- random[names(obj@flist)]
959              rmmats <- c(lapply(random, "[[", 1))
960              vc <- VarCorr(obj)
961              ans <- fxd + matrix(rnorm(n * nsamp, sd = vc@scale), nr = n)
962              ans
963          })

Legend:
Removed from v.864  
changed lines
  Added in v.878

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge