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 831, Mon Aug 8 08:24:10 2005 UTC revision 832, Mon Aug 8 08:33:55 2005 UTC
# Line 51  Line 51 
51          return(term)          return(term)
52      }      }
53      stopifnot(length(term) == 3)      stopifnot(length(term) == 3)
54      if (is.call(term) && term[[1]] == as.name('|')) term[[1]] <- as.name('+')      if (is.call(term) && term[[1]] == as.name('|'))
55            term[[1]] <- as.name('+')
56      term[[2]] <- subbars(term[[2]])      term[[2]] <- subbars(term[[2]])
57      term[[3]] <- subbars(term[[3]])      term[[3]] <- subbars(term[[3]])
58      term      term
# Line 71  Line 72 
72    
73  ## Control parameters for lmer  ## Control parameters for lmer
74  lmerControl <-  lmerControl <-
75    function(maxIter = 200,    function(maxIter = 200, # used in ../src/lmer.c only
76               tolerance = sqrt(.Machine$double.eps),# dito
77             msMaxIter = 200,             msMaxIter = 200,
78             tolerance = sqrt((.Machine$double.eps)),             ## msTol = sqrt(.Machine$double.eps),
79             niterEM = 15,             ## FIXME:  should be able to pass tolerances to nlminb()
            msTol = sqrt(.Machine$double.eps),  
80             msVerbose = getOption("verbose"),             msVerbose = getOption("verbose"),
81             PQLmaxIt = 30,             niterEM = 15,
82             EMverbose = getOption("verbose"),             EMverbose = getOption("verbose"),
83               PQLmaxIt = 30,# unused _FIXME_
84             analyticGradient = TRUE,             analyticGradient = TRUE,
85             analyticHessian = FALSE)             analyticHessian = FALSE # unused _FIXME_
86               )
87  {  {
88      list(maxIter = as.integer(maxIter),      list(maxIter = as.integer(maxIter),
          msMaxIter = as.integer(msMaxIter),  
89           tolerance = as.double(tolerance),           tolerance = as.double(tolerance),
90             msMaxIter = as.integer(msMaxIter),
91             ## msTol = as.double(msTol),
92             msVerbose = as.integer(msVerbose),# "integer" on purpose
93           niterEM = as.integer(niterEM),           niterEM = as.integer(niterEM),
          msTol = as.double(msTol),  
          msVerbose = as.logical(msVerbose),  
          PQLmaxIt = as.integer(PQLmaxIt),  
94           EMverbose = as.logical(EMverbose),           EMverbose = as.logical(EMverbose),
95             PQLmaxIt = as.integer(PQLmaxIt),
96           analyticGradient = as.logical(analyticGradient),           analyticGradient = as.logical(analyticGradient),
97           analyticHessian = as.logical(analyticHessian))           analyticHessian = as.logical(analyticHessian))
98  }  }
# Line 100  Line 103 
103                     control = list(),                     control = list(),
104                     subset, weights, na.action, offset,                     subset, weights, na.action, offset,
105                     model = TRUE, x = FALSE, y = FALSE, ...)                     model = TRUE, x = FALSE, y = FALSE, ...)
106        {                                 ## match and check parameters            ## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
107          {
108              ## match and check parameters
109            if (length(formula) < 3) stop("formula must be a two-sided formula")            if (length(formula) < 3) stop("formula must be a two-sided formula")
110            cv <- do.call("lmerControl", control)            cv <- do.call("lmerControl", control)
           if (lmm <- missing(family)) { # linear mixed model  
               method <- match.arg(method)  
               if (method %in% c("PQL", "Laplace", "AGQ")) {  
                   warning(paste('Argument method = "', method,  
                                 '" is not meaningful for a linear mixed model.\n',  
                                 'Using method = "REML".\n', sep = ''))  
                   method <- "REML"  
               }  
           } else {                      # generalized linear mixed model  
               method <- if (missing(method)) "PQL" else match.arg(method)  
               if (method == "ML") method <- "PQL"  
               if (method == "REML")  
                   warning(paste('Argument method = "REML" is not meaningful',  
                                 'for a generalized linear mixed model.',  
                                 '\nUsing method = "PQL".\n'))  
           }  
   
111            ## evaluate glm.fit, a generalized linear fit of fixed effects only            ## evaluate glm.fit, a generalized linear fit of fixed effects only
112            mf <- match.call()            mf <- match.call()
113            m <- match(c("family", "data", "subset", "weights",            m <- match(c("family", "data", "subset", "weights",
# Line 136  Line 124 
124            glm.fit <- eval(mf, parent.frame())            glm.fit <- eval(mf, parent.frame())
125            x <- glm.fit$x            x <- glm.fit$x
126            y <- as.double(glm.fit$y)            y <- as.double(glm.fit$y)
127              lmm <- missing(family) # => linear mixed model
128            family <- glm.fit$family            family <- glm.fit$family
129    
130              if (lmm) { # linear mixed model
131                  method <- match.arg(method)
132                  if (method %in% c("PQL", "Laplace", "AGQ")) {
133                      warning(paste('Argument method = "', method,
134                                    '" is not meaningful for a linear mixed model.\n',
135                                    'Using method = "REML".\n', sep = ''))
136                      method <- "REML"
137                  }
138              } else { # generalized linear mixed model
139                  if (missing(method)) method <- "PQL"
140                  else {
141                      method <- match.arg(method)
142    
143                      ## <FIXME>: only
144                      ## if( !(family == "gaussian" && link == "identity") ) {
145                      if (method == "ML") method <- "PQL"
146                      if (method == "REML")
147                          warning('Argument method = "REML" is not meaningful ',
148                                  'for a generalized linear mixed model.',
149                                  '\nUsing method = "PQL".\n')
150                      ## } </FIXME>
151                  }
152              }
153    
154            ## evaluate a model frame for fixed and random effects            ## evaluate a model frame for fixed and random effects
155            mf$formula <- frame.form            mf$formula <- frame.form
156            mf$x <- mf$model <- mf$y <- mf$family <- NULL            mf$x <- mf$model <- mf$y <- mf$family <- NULL
# Line 148  Line 161 
161            ## grouping factors and model matrices for random effects            ## grouping factors and model matrices for random effects
162            bars <- findbars(formula[[3]])            bars <- findbars(formula[[3]])
163            random <-            random <-
164                lapply(bars,                lapply(bars, function(x)
165                       function(x) list(model.matrix(eval(substitute(~term,                       list(model.matrix(eval(substitute(~ T, list(T = x[[2]]))),
                                                                    list(term=x[[2]]))),  
166                                                     frm),                                                     frm),
167                                        eval(substitute(as.factor(fac)[,drop = TRUE],                                        eval(substitute(as.factor(fac)[,drop = TRUE],
168                                                        list(fac = x[[3]])), frm)))                                                        list(fac = x[[3]])), frm)))
# Line 267  Line 279 
279    
280            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
281            loglik[] <- -deviance/2            loglik[] <- -deviance/2
282            new("lmer", mer, frame = frm, terms = glm.fit$terms,            new("lmer", mer,
283                  frame = if (model) frm else data.frame(),
284                  terms = glm.fit$terms,
285                assign = attr(glm.fit$x, "assign"),                assign = attr(glm.fit$x, "assign"),
286                call = match.call(), family = family,                call = match.call(), family = family,
287                logLik = loglik, fixed = fxd)                logLik = loglik, fixed = fxd)
288        })        })
289    ## end{  "lmer . formula " }
290    
291  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
292                   function(x, value)                   function(x, value)
# Line 282  Line 297 
297                                           function(k) 1:((k*(k+1))/2) <= k))                                           function(k) 1:((k*(k+1))/2) <= k))
298                   fn <- function(pars)                   fn <- function(pars)
299                       deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))                       deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
                  gr <- NULL  
                  if (value$analyticGradient)  
300                       gr <-                       gr <-
301                         if (value$analyticGradient)
302                           function(pars) {                           function(pars) {
303                               if (!isTRUE(all.equal(pars,                               if (!isTRUE(all.equal(pars,
304                                                     .Call("lmer_coef", x,                                                     .Call("lmer_coef", x,
# Line 292  Line 306 
306                                   .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")                                   .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
307                               .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")                               .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
308                           }                           }
309                     ## else NULL
310                   optimRes <-                   optimRes <-
311                       if (exists("nlminb", mode = "function"))                       if (exists("nlminb", mode = "function"))
312                           nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),                           nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
# Line 737  Line 752 
752  setMethod("plot", signature(x = "lmer.coef"),  setMethod("plot", signature(x = "lmer.coef"),
753            function(x, y, ...)            function(x, y, ...)
754        {        {
755            if (require("lattice", quietly = TRUE)) {            ## require("lattice", quietly = TRUE) -- now via Imports
756                varying <- unique(do.call("c",                varying <- unique(do.call("c",
757                                          lapply(x, function(el)                                          lapply(x, function(el)
758                                                 names(el)[sapply(el,                                                 names(el)[sapply(el,
# Line 752  Line 767 
767                                              list(y = as.name(varying[1]),                                              list(y = as.name(varying[1]),
768                                                   x = as.name(varying[2])))), gf, ...),                                                   x = as.name(varying[2])))), gf, ...),
769                       splom(~ gf | .grp, ...))                       splom(~ gf | .grp, ...))
           }  
770        })        })
771    
772  setMethod("plot", signature(x = "lmer.ranef"),  setMethod("plot", signature(x = "lmer.ranef"),
773            function(x, y, ...)            function(x, y, ...)
774        {        {
775            if (require("lattice", quietly = TRUE))            ## require("lattice", quietly = TRUE) -- now via Imports
776                lapply(x, function(x) {                lapply(x, function(x) {
777                    cn <- lapply(colnames(x), as.name)                    cn <- lapply(colnames(x), as.name)
778                    switch(min(ncol(x), 3),                    switch(min(ncol(x), 3),
779                           qqmath(eval(substitute(~ x,                       qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
                                                 list(x = cn[[1]]))),  
                                 x, ...),  
780                           xyplot(eval(substitute(y ~ x,                           xyplot(eval(substitute(y ~ x,
781                                                  list(y = cn[[1]],                                                  list(y = cn[[1]],
782                                                       x = cn[[2]]))),                                                   x = cn[[2]]))), x, ...),
                                 x, ...),  
783                           splom(~ x, ...))                           splom(~ x, ...))
784                })                })
785        })        })

Legend:
Removed from v.831  
changed lines
  Added in v.832

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