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 434, Thu Jan 13 19:36:12 2005 UTC revision 435, Thu Jan 13 19:37:26 2005 UTC
# Line 1  Line 1 
1    lmerControl <-                            # Control parameters for lmer
2      function(maxIter = 50,
3               msMaxIter = 50,
4               tolerance = sqrt((.Machine$double.eps)),
5               niterEM = 20,
6               msTol = sqrt(.Machine$double.eps),
7               msVerbose = getOption("verbose"),
8               PQLmaxIt = 20,
9               .relStep = (.Machine$double.eps)^(1/3),
10               EMverbose = getOption("verbose"),
11               analyticGradient = TRUE,
12               analyticHessian=FALSE)
13    {
14        list(maxIter = maxIter,
15             msMaxIter = msMaxIter,
16             tolerance = tolerance,
17             niterEM = niterEM,
18             msTol = msTol,
19             msVerbose = msVerbose,
20             PQLmaxIt = PQLmaxIt,
21             .relStep = .relStep,
22             EMverbose=EMverbose,
23             analyticHessian=analyticHessian,
24             analyticGradient=analyticGradient)
25    }
26    
27    setMethod("lmer", signature(formula = "formula"),
28              function(formula, data,
29                       method = c("REML", "ML"),
30                       control = list(),
31                       subset, weights, na.action, offset,
32                       model = TRUE, x = FALSE, y = FALSE, ...)
33          {
34                                            # match and check parameters
35              method <- match.arg(method)
36              controlvals <- do.call("lmerControl", control)
37              controlvals$REML <- method == "REML"
38              if (length(formula) < 3) stop("formula must be a two-sided formula")
39                                            # create the model frame as frm
40              mf <- match.call()
41              m <- match(c("data", "subset", "weights", "na.action", "offset"),
42                         names(mf), 0)
43              mf <- mf[c(1, m)]
44              mf[[1]] <- as.name("model.frame")
45              frame.form <- subbars(formula)
46              environment(frame.form) <- environment(formula)
47              mf$formula <- frame.form
48              mf$drop.unused.levels <- TRUE
49              frm <- eval(mf, parent.frame())
50    
51              ## grouping factors and model matrices for random effects
52              bars <- findbars(formula[[3]])
53              random <-
54                  lapply(bars,
55                         function(x) list(model.matrix(eval(substitute(~term,
56                                                                       list(term=x[[2]]))),
57                                                       frm),
58                                          eval(substitute(as.factor(fac),
59                                                          list(fac = x[[3]])), frm)))
60              names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
61    
62              ## order factor list by decreasing number of levels
63              ford <- rev(order(sapply(random, function(x) length(levels(x[[2]])))))
64              if (any(ford != seq(a = random))) { # re-order both facs and random
65                  random <- random[ford]
66              }
67              mmats <- c(lapply(random, "[[", 1),
68                         .fixed = list(cbind(model.matrix(nobars(formula), frm),
69                         .response = model.response(frm))))
70              obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")
71              .Call("lmer_initial", obj, PACKAGE="Matrix")
72              .Call("lmer_ECMEsteps", obj,
73                    controlvals$niterEM,
74                    controlvals$REML,
75                    controlvals$EMverbose,
76                    PACKAGE = "Matrix")
77              LMEoptimize(obj) <- controlvals
78              obj@call <- match.call()
79              #fitted = .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")
80              #residuals = mmats$.Xy[,".response"] - fitted
81              #if (as.logical(x)[1]) x = mmats else x = list()
82              #rm(mmats)
83              obj
84          })
85    
86  setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),  setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),
87                   function(x, value)                   function(x, value)
88               {               {

Legend:
Removed from v.434  
changed lines
  Added in v.435

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