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 887, Wed Aug 31 20:31:47 2005 UTC revision 888, Thu Sep 1 13:46:21 2005 UTC
# Line 73  Line 73 
73  ## Control parameters for lmer  ## Control parameters for lmer
74  lmerControl <-  lmerControl <-
75    function(maxIter = 200, # used in ../src/lmer.c only    function(maxIter = 200, # used in ../src/lmer.c only
76             tolerance = sqrt(.Machine$double.eps),# dito             tolerance = sqrt(.Machine$double.eps),# ditto
77             msMaxIter = 200,             msMaxIter = 200,
78             ## msTol = sqrt(.Machine$double.eps),             ## msTol = sqrt(.Machine$double.eps),
79             ## FIXME:  should be able to pass tolerances to nlminb()             ## FIXME:  should be able to pass tolerances to nlminb()
# Line 474  Line 474 
474                invisible(object)                invisible(object)
475            })            })
476    
   
   
   
477  ## calculates degrees of freedom for fixed effects Wald tests  ## calculates degrees of freedom for fixed effects Wald tests
478  ## This is a placeholder.  The answers are generally wrong.  It will  ## This is a placeholder.  The answers are generally wrong.  It will
479  ## be very tricky to decide what a 'right' answer should be with  ## be very tricky to decide what a 'right' answer should be with
# Line 929  Line 926 
926  rWishart <- function(n, df, invScal)  rWishart <- function(n, df, invScal)
927    .Call("Matrix_rWishart", n, df, invScal)    .Call("Matrix_rWishart", n, df, invScal)
928    
929  setMethod("simulate", signature(object = "lmer"),  
930            function(object, n = 1, ...)  setMethod("model.matrix", signature(object = "lmer"),
931              function(object, ...)
932        {        {
           family <- object@family  
           if (family$family != "gaussian" ||  
               family$link != "identity")  
               stop("simulation of generalized linear mixed models not implemented yet")  
           ## create the mean from the fixed effects  
933            frm <- object@frame            frm <- object@frame
934            fixed.form <- Matrix:::nobars(object@call$formula)            fixed.form <- Matrix:::nobars(object@call$formula)
935            if (inherits(fixed.form, "name")) # RHS is empty - use a constant            if (inherits(fixed.form, "name")) # RHS is empty - use a constant
936                fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))                fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
937            glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE, y = TRUE)            glm.fit <- glm(eval(fixed.form), object@family, frm, x = TRUE, y = TRUE)
938            fxd <- unname(drop(glm.fit$x %*% fixef(object)))            fxd <- unname(drop(glm.fit$x %*% fixef(object)))
939    
940              ## Create the random effects model matrices
941            bars <- Matrix:::findbars(object@call$formula[[3]])            bars <- Matrix:::findbars(object@call$formula[[3]])
942            random <-            random <-
943                lapply(bars,                lapply(bars,
# Line 955  Line 950 
950            ## re-order the random effects pairs if necessary            ## re-order the random effects pairs if necessary
951            if (any(names(random) != names(object@flist)))            if (any(names(random) != names(object@flist)))
952                random <- random[names(object@flist)]                random <- random[names(object@flist)]
953            mer <- as(object, "mer")            c(lapply(random, "[[", 1),
954                .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
955          })
956    
957    setMethod("simulate", signature(object = "lmer"),
958              function(object, nsim = 1,
959                       seed = runif(1, 0, .Machine$integer.max),
960                       ...)
961          {
962              runif(1) ## to initialize the RNG if necessary
963              RNGstate <- .Random.seed
964              set.seed(seed)
965    
966              family <- object@family
967              if (family$family != "gaussian" ||
968                  family$link != "identity")
969                  stop("simulation of generalized linear mixed models not yet implemented")
970    
971              ## pieces we will need later
972              scale <- .Call("lmer_sigma", object, object@method == "REML",
973                             PACKAGE = "Matrix")
974              mmats <- model.matrix(object)
975              ff <- fixef(object)
976    
977    ###_FIXME: If the factor levels have been permuted, has the
978    ###        permutation been applied in the stored frame?  Otherwise we
979    ###        need to check this.
980    
981            ## similate the linear predictors            ## similate the linear predictors
982            lpred <- .Call("lmer_simulate", mer, n, fxd,            lpred <- .Call("lmer_simulate", as(object, "mer"), nsim,
983                           lapply(random, "[[", 1), TRUE, PACKAGE = "Matrix")                           unname(drop(mmats[[length(mmats)]][,seq(a = ff)] %*% ff)),
984                             mmats, TRUE, PACKAGE = "Matrix")
985            ## add per-observation noise term            ## add per-observation noise term
986            lpred + rnorm(prod(dim(lpred)),            lpred <- as.data.frame(lpred + rnorm(prod(dim(lpred)), sd = scale))
987                          sd = .Call("lmer_sigma", mer,  
988                          object@method == "REML", PACKAGE = "Matrix"))            ## save the seed and restore the RNG state
989              attr(lpred, "seed") <- seed
990              assign(".Random.seed", RNGstate, envir = .GlobalEnv)
991              lpred
992        })        })
993    
994  simulate2 <- function(object, n = 1, ...)  simulate2 <- function(object, n = 1, ...)
# Line 1023  Line 1049 
1049    
1050  refdist <- function(fm1, fm2, n, ...)  refdist <- function(fm1, fm2, n, ...)
1051  {  {
1052      frm <- fm1@frame      cv <- lmerControl()
1053        obs <- deviance(fm2) - deviance(fm1)
1054      newy <- simulate(fm2, n)      newy <- simulate(fm2, n)
1055        mm1 <- model.matrix(fm1)
1056        mm2 <- model.matrix(fm2)
1057      ref <- numeric(n)      ref <- numeric(n)
1058      f1 <- eval(fm1@call$formula)      mer1 <- as(fm1, "mer")
1059      f2 <- eval(fm2@call$formula)      mer2 <- as(fm2, "mer")
     nm <- as.character(as.name(f1[[2]]))  
1060      for (j in 1:n) {      for (j in 1:n) {
1061          frm[[nm]] <- newy[, j]          .Call("lmer_update_y", mer2, newy[[j]], mm2, PACKAGE = "Matrix")
1062          ref[j] <- deviance(lmer(f2, frm)) - deviance(lmer(f1, frm))          LMEoptimize(mer2) <- cv
1063            .Call("lmer_update_y", mer1, newy[[j]], mm1, PACKAGE = "Matrix")
1064            LMEoptimize(mer1) <- cv
1065            ref[j] <- deviance(mer2) - deviance(mer1)
1066      }      }
1067      attr(ref, "observed") <- deviance(fm2) - deviance(fm1)      attr(ref, "observed") <- obs
1068      ref      ref
1069  }  }

Legend:
Removed from v.887  
changed lines
  Added in v.888

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