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 984, Tue Oct 11 18:08:46 2005 UTC revision 986, Fri Oct 21 19:36:30 2005 UTC
# Line 1090  Line 1090 
1090      ref      ref
1091  }  }
1092    
 mer2 <-  
     function(formula, data, family,  
              method = c("REML", "ML", "PQL", "Laplace", "AGQ"),  
              control = list(), start,  
              subset, weights, na.action, offset,  
              model = TRUE, x = FALSE, y = FALSE , ...)  
           ## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(  
 {  
     ## match and check parameters  
     if (length(formula) < 3) stop("formula must be a two-sided formula")  
     ## cv <- do.call("Matrix:::lmerControl", control)  
     cv <- control  
     cv$analyticGradient <- FALSE  
     cv$msMaxIter <- as.integer(200)  
     if (is.null(cv$msVerbose)) cv$msVerbose <- as.integer(1)  
   
     ## Must evaluate the model frame first and then fit the glm using  
     ## that frame.  Otherwise missing values in the grouping factors  
     ## cause inconsistent numbers of observations.  
   
     ## evaluate glm.fit, a generalized linear fit of fixed effects only  
     mf <- match.call()  
     m <- match(c("family", "data", "subset", "weights",  
                  "na.action", "offset"), names(mf), 0)  
     mf <- fe <- mf[c(1, m)]  
     frame.form <- subbars(formula) # substitute `+' for `|'  
     fixed.form <- nobars(formula)  # remove any terms with `|'  
     if (inherits(fixed.form, "name")) # RHS is empty - use a constant  
         fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))  
     environment(fixed.form) <- environment(frame.form) <- environment(formula)  
   
     ## evaluate a model frame for fixed and random effects  
     mf$formula <- frame.form  
     mf$family <- NULL  
     mf$drop.unused.levels <- TRUE  
     mf[[1]] <- as.name("model.frame")  
     frm <- eval(mf, parent.frame())  
   
     ## fit a glm model to the fixed formula  
     fe$formula <- fixed.form  
     fe$data <- frm  
     fe$x <- fe$model <- fe$y <- TRUE  
     fe[[1]] <- as.name("glm")  
     glm.fit <- eval(fe, parent.frame())  
     x <- glm.fit$x  
     y <- as.double(glm.fit$y)  
     family <- glm.fit$family  
   
     ## check for a linear mixed model  
     lmm <- family$family == "gaussian" && family$link == "identity"  
     if (lmm) { # 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  
         if (missing(method)) method <- "PQL"  
         else {  
             method <- match.arg(method)  
             if (method == "ML") method <- "PQL"  
             if (method == "REML")  
                 warning('Argument method = "REML" is not meaningful ',  
                         'for a generalized linear mixed model.',  
                         '\nUsing method = "PQL".\n')  
         }  
     }  
   
     ## grouping factors and model matrices for random effects  
     bars <- findbars(formula[[3]])  
     random <-  
         lapply(bars, function(x)  
                list(t(model.matrix(eval(substitute(~ expr,  
                                                    list(expr = x[[2]]))),  
                                    frm)),  
                     eval(substitute(as.factor(fac)[,drop = TRUE],  
                                     list(fac = x[[3]])), frm)))  
     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))  
   
     ## order factor list by decreasing number of levels  
     nlev <- sapply(random, function(x) length(levels(x[[2]])))  
     if (any(diff(nlev) > 0)) {  
         random <- random[rev(order(nlev))]  
     }  
   
     ## Create the model matrices and a mixed-effects representation (mer)  
     mer <- .Call("mer2_create", random, x, y, method, PACKAGE = "Matrix")  
     LMEoptimize(mer) <- cv  
     mer  
 }  
   
 ## Extract the permutation  
 setAs("mer2", "pMatrix", function(from) .Call("mer2_pMatrix", from))  
   
 ## Extract the L matrix  
 setAs("mer2", "dtCMatrix", function(from) .Call("mer2_dtCMatrix", from))  
   
 ## Optimization for mer2 objects  
 setReplaceMethod("LMEoptimize", signature(x="mer2", value="list"),  
                  function(x, value)  
              {  
                  if (value$msMaxIter < 1) return(x)  
                  nc <- x@nc  
                  constr <- unlist(lapply(nc[1:(length(nc) - 1)],  
                                          function(k) 1:((k*(k+1))/2) <= k))  
                  fn <- function(pars)  
                      deviance(.Call("mer2_coefGets", x, pars, 2, PACKAGE = "Matrix"))  
                  gr <- NULL  ## No gradient yet  
 ##                     if (value$analyticGradient)  
 ##                         function(pars) {  
 ##                             if (!isTRUE(all.equal(pars,  
 ##                                                   .Call("lmer_coef", x,  
 ##                                                         2, PACKAGE = "Matrix"))))  
 ##                                 .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")  
 ##                             .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")  
 ##                         }  
                  ## else NULL  
                  optimRes <- nlminb(.Call("mer2_coef", x, 2, PACKAGE = "Matrix"),  
                                     fn, gr,  
                                     lower = ifelse(constr, 5e-10, -Inf),  
                                     control = list(iter.max = value$msMaxIter,  
                                     trace = as.integer(value$msVerbose)))  
                  .Call("mer2_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")  
                  if (optimRes$convergence != 0) {  
                      warning(paste("nlminb returned message",  
                                    optimRes$message,"\n"))  
                  }  
                  return(x)  
              })  
   
 setMethod("deviance", "mer2",  
           function(object, REML = NULL, ...) {  
               .Call("mer2_factor", object, PACKAGE = "Matrix")  
               if (is.null(REML))  
                   REML <- object@method == "REML"  
               object@deviance[[ifelse(REML, "REML", "ML")]]  
           })  

Legend:
Removed from v.984  
changed lines
  Added in v.986

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