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 688, Sat Apr 2 13:18:16 2005 UTC revision 689, Sat Apr 2 14:41:46 2005 UTC
# Line 1  Line 1 
1    ## Methods for lmer and for the objects that it produces
2    
3    ## Some utilities
4    
5  contr.SAS <- function(n, contrasts = TRUE)  contr.SAS <- function(n, contrasts = TRUE)
6  {  ## Eliminate this function after R-2.1.0 is released
7      contr.treatment(n, if (is.numeric(n) && length(n) == 1) n else length(n), contrasts)      contr.treatment(n,
8                        if (is.numeric(n) && length(n) == 1) n else length(n),
9                        contrasts)
10    
11    Lind <- function(i,j) {
12        if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))
13        ((i - 1) * i)/2 + j
14    }
15    
16    Dhalf <- function(from) {
17        D <- from@D
18        nf <- length(D)
19        Gp <- from@Gp
20        res <- array(0, rep(Gp[nf+1],2))
21        for (i in 1:nf) {
22            DD <- D[[i]]
23            dd <- dim(DD)
24            for (k in 1:dd[3]) {
25                mm <- array(DD[ , , k], dd[1:2])
26                base <- Gp[i] + (k - 1)*dd[1]
27                res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)
28            }
29        }
30        res
31  }  }
32    
33  lmerControl <-                            # Control parameters for lmer  lmerControl <-                            # Control parameters for lmer
# Line 11  Line 38 
38             msTol = sqrt(.Machine$double.eps),             msTol = sqrt(.Machine$double.eps),
39             msVerbose = getOption("verbose"),             msVerbose = getOption("verbose"),
40             PQLmaxIt = 20,             PQLmaxIt = 20,
            .relStep = (.Machine$double.eps)^(1/3),  
41             EMverbose = getOption("verbose"),             EMverbose = getOption("verbose"),
42             analyticGradient = TRUE,             analyticGradient = TRUE,
43             analyticHessian=FALSE)             analyticHessian=FALSE)
# Line 23  Line 49 
49           msTol = msTol,           msTol = msTol,
50           msVerbose = msVerbose,           msVerbose = msVerbose,
51           PQLmaxIt = PQLmaxIt,           PQLmaxIt = PQLmaxIt,
          .relStep = .relStep,  
52           EMverbose=EMverbose,           EMverbose=EMverbose,
53           analyticHessian=analyticHessian,           analyticHessian=analyticHessian,
54           analyticGradient=analyticGradient)           analyticGradient=analyticGradient)
55  }  }
56    
57  setMethod("lmer", signature(formula = "formula"),  setMethod("lmer", signature(formula = "formula", family = "missing"),
58            function(formula, data,            function(formula, data, family,
59                     method = c("REML", "ML"),                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
60                     control = list(),                     control = list(),
61                     subset, weights, na.action, offset,                     subset, weights, na.action, offset,
62                     model = TRUE, x = FALSE, y = FALSE, ...)                     model = TRUE, x = FALSE, y = FALSE, ...)
# Line 76  Line 101 
101                       .fixed = list(cbind(Xmat, .response = model.response(frm))))                       .fixed = list(cbind(Xmat, .response = model.response(frm))))
102            ## FIXME: Use Xfrm and Xmat to get the terms and assign            ## FIXME: Use Xfrm and Xmat to get the terms and assign
103            ## slots, pass them to lmer_create, then destroy them            ## slots, pass them to lmer_create, then destroy them
104            obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")            obj <- .Call("lmer_create", lapply(random, "[[", 2),
105            obj@terms <- attr(model.frame(fixed.form, data), "terms")                         mmats, PACKAGE = "Matrix")
106            obj@assign <- attr(Xmat, "assign")            slot(obj, "terms") <- attr(model.frame(fixed.form, data), "terms")
107            obj@call <- match.call()            slot(obj, "assign") <- attr(Xmat, "assign")
108            obj@REML <- REML            slot(obj, "call") <- match.call()
109              slot(obj, "REML") <- REML
110            rm(Xmat)            rm(Xmat)
111            .Call("lmer_initial", obj, PACKAGE="Matrix")            .Call("lmer_initial", obj, PACKAGE="Matrix")
112            .Call("lmer_ECMEsteps", obj,            .Call("lmer_ECMEsteps", obj,
# Line 89  Line 115 
115                  controlvals$EMverbose,                  controlvals$EMverbose,
116                  PACKAGE = "Matrix")                  PACKAGE = "Matrix")
117            LMEoptimize(obj) <- controlvals            LMEoptimize(obj) <- controlvals
118            #fitted <- .Call("lmer_fitted", obj, mmats, TRUE, PACKAGE = "Matrix")            slot(obj, "residuals") <-
119            #residuals <- mmats$.Xy[,".response"] - fitted                unname(model.response(frm) -
120                         (slot(obj, "fitted") <-
121                          .Call("lmer_fitted", obj, mmats, TRUE, PACKAGE = "Matrix")))
122            obj            obj
123        })        })
124    
# Line 126  Line 154 
154               })               })
155    
156  setMethod("ranef", signature(object = "lmer"),  setMethod("ranef", signature(object = "lmer"),
157            function(object, ...) {            function(object, accumulate = FALSE, ...) {
158                .Call("lmer_ranef", object, PACKAGE = "Matrix")                val <- new("lmer.ranef",
159                             lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),
160                                    data.frame, check.names = FALSE),
161                             varFac = object@bVar,
162                             stdErr = .Call("lmer_sigma", object,
163                             object@REML, PACKAGE = "Matrix"))
164                  if (!accumulate || length(val@varFac) == 1) return(val)
165                  ## check for nested factors
166                  L <- object@L
167                  if (any(sapply(seq(a = val), function(i) length(L[[Lind(i,i)]]@i))))
168                      error("Require nested grouping factors to accumulate random effects")
169                  val
170            })            })
171    
172  setMethod("fixef", signature(object = "lmer"),  setMethod("fixef", signature(object = "lmer"),
# Line 243  Line 282 
282                invisible(object)                invisible(object)
283            })            })
284    
285    setMethod("lmer", signature(formula = "formula"),
286              function(formula, family, data,
287                       method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
288                       control = list(),
289                       subset, weights, na.action, offset,
290                       model = TRUE, x = FALSE, y = FALSE, ...)
291          {
292              gVerb <- getOption("verbose")
293                                            # match and check parameters
294              controlvals <- do.call("lmerControl", control)
295              controlvals$REML <- FALSE
296              if (length(formula) < 3) stop("formula must be a two-sided formula")
297    
298              ## initial glm fit
299              mf <- match.call()
300              m <- match(c("family", "data", "subset", "weights",
301                           "na.action", "offset"),
302                         names(mf), 0)
303              mf <- mf[c(1, m)]
304              mf[[1]] <- as.name("glm")
305              fixed.form <- nobars(formula)
306              if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula
307              environment(fixed.form) <- environment(formula)
308              mf$formula <- fixed.form
309              mf$x <- mf$model <- mf$y <- TRUE
310              glm.fit <- eval(mf, parent.frame())
311              family <- glm.fit$family
312              ## Note: offset is on the linear scale
313              offset <- glm.fit$offset
314              if (is.null(offset)) offset <- 0
315              weights <- sqrt(abs(glm.fit$prior.weights))
316              ## initial 'fitted' values on linear scale
317              etaold <- eta <- glm.fit$linear.predictors
318    
319              ## evaluation of model frame
320              mf$x <- mf$model <- mf$y <- mf$family <- NULL
321              mf$drop.unused.levels <- TRUE
322              this.form <- subbars(formula)
323              environment(this.form) <- environment(formula)
324              mf$formula <- this.form
325              mf[[1]] <- as.name("model.frame")
326              frm <- eval(mf, parent.frame())
327    
328              ## grouping factors and model matrices for random effects
329              bars <- findbars(formula[[3]])
330              random <-
331                  lapply(bars,
332                         function(x) list(model.matrix(eval(substitute(~term,
333                                                                       list(term=x[[2]]))),
334                                                       frm),
335                                          eval(substitute(as.factor(fac)[,drop = TRUE],
336                                                          list(fac = x[[3]])), frm)))
337              names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
338    
339              ## order factor list by decreasing number of levels
340              nlev <- sapply(random, function(x) length(levels(x[[2]])))
341              if (any(diff(nlev) > 0)) {
342                  random <- random[rev(order(nlev))]
343              }
344              mmats <- c(lapply(random, "[[", 1),
345                         .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
346              ## FIXME: Use Xfrm and Xmat to get the terms and assign
347              ## slots, pass these to lmer_create, then destroy Xfrm, Xmat, etc.
348              obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")
349              obj@terms <- attr(glm.fit$model, "terms")
350              obj@assign <- attr(glm.fit$x, "assign")
351              obj@call <- match.call()
352              obj@REML <- FALSE
353              rm(glm.fit)
354              .Call("lmer_initial", obj, PACKAGE="Matrix")
355              mmats.unadjusted <- mmats
356              mmats[[1]][1,1] <- mmats[[1]][1,1]
357              conv <- FALSE
358              firstIter <- TRUE
359              msMaxIter.orig <- controlvals$msMaxIter
360              responseIndex <- ncol(mmats$.fixed)
361    
362              for (iter in seq(length = controlvals$PQLmaxIt))
363              {
364                  mu <- family$linkinv(eta)
365                  dmu.deta <- family$mu.eta(eta)
366                  ## weights (note: weights is already square-rooted)
367                  w <- weights * dmu.deta / sqrt(family$variance(mu))
368                  ## adjusted response (should be comparable to X \beta, not including offset
369                  z <- eta - offset + (mmats.unadjusted$.fixed[, responseIndex] - mu) / dmu.deta
370                  .Call("nlme_weight_matrix_list",
371                        mmats.unadjusted, w, z, mmats, PACKAGE="Matrix")
372                  .Call("lmer_update_mm", obj, mmats, PACKAGE="Matrix")
373                  if (firstIter) {
374                      .Call("lmer_initial", obj, PACKAGE="Matrix")
375                      if (gVerb) cat(" PQL iterations convergence criterion\n")
376                  }
377                  .Call("lmer_ECMEsteps", obj,
378                        controlvals$niterEM,
379                        FALSE,
380                        controlvals$EMverbose,
381                        PACKAGE = "Matrix")
382                  LMEoptimize(obj) <- controlvals
383                  eta[] <- offset + ## FIXME: should the offset be here ?
384                      .Call("lmer_fitted", obj,
385                            mmats.unadjusted, TRUE, PACKAGE = "Matrix")
386                  crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))
387                  if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))
388                  ## use this to determine convergence
389                  if (crit < controlvals$tolerance) {
390                      conv <- TRUE
391                      break
392                  }
393                  etaold[] <- eta
394    
395                  ## Changing number of iterations on second and
396                  ## subsequent iterations.
397                  if (firstIter)
398                  {
399                      controlvals$niterEM <- 2
400                      controlvals$msMaxIter <- 10
401                      firstIter <- FALSE
402                  }
403              }
404              if (!conv) warning("IRLS iterations for glmm did not converge")
405              obj
406          })
407    
408  ## calculates degrees of freedom for fixed effects Wald tests  ## calculates degrees of freedom for fixed effects Wald tests
409  ## This is a placeholder.  The answers are generally wrong.  It will  ## This is a placeholder.  The answers are generally wrong.  It will
410  ## 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 262  Line 424 
424                val <- -deviance(object, REML = REML)/2                val <- -deviance(object, REML = REML)/2
425                nc <- object@nc[-seq(a = object@Omega)]                nc <- object@nc[-seq(a = object@Omega)]
426                attr(val, "nall") <- attr(val, "nobs") <- nc[2]                attr(val, "nall") <- attr(val, "nobs") <- nc[2]
427                attr(val, "df") <- nc[1] + length(coef(object))                attr(val, "df") <- nc[1] + length(ccoef(object))
428                attr(val, "REML") <- REML                attr(val, "REML") <- REML
429                class(val) <- "logLik"                class(val) <- "logLik"
430                val                val
# Line 384  Line 546 
546            ci            ci
547        })        })
548    
549  setReplaceMethod("coef", signature(object = "lmer", value = "numeric"),  setMethod("param", signature(object = "lmer"),
                  function(object, unconst = FALSE, ..., value)  
                  .Call("lmer_coefGets", object, as.double(value),  
                        unconst, PACKAGE = "Matrix"))  
   
 setMethod("coef", signature(object = "lmer"),  
550            function(object, unconst = FALSE, ...) {            function(object, unconst = FALSE, ...) {
551                .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")                .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
552            })            })
# Line 430  Line 587 
587                rr                rr
588            })            })
589    
 Lind <- function(i,j) {  
     if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))  
     ((i - 1) * i)/2 + j  
 }  
   
 Dhalf <- function(from) {  
     D <- from@D  
     nf <- length(D)  
     Gp <- from@Gp  
     res <- array(0, rep(Gp[nf+1],2))  
     for (i in 1:nf) {  
         DD <- D[[i]]  
         dd <- dim(DD)  
         for (k in 1:dd[3]) {  
             mm <- array(DD[ , , k], dd[1:2])  
             base <- Gp[i] + (k - 1)*dd[1]  
             res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)  
         }  
     }  
     res  
 }  
   
590  ## Extract the L matrix  ## Extract the L matrix
591  setAs("lmer", "dtTMatrix",  setAs("lmer", "dtTMatrix",
592        function(from)        function(from)
# Line 509  Line 644 
644        new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx,        new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx,
645            uplo = "U")            uplo = "U")
646    })    })
647    
648    setMethod("fitted", signature(object = "lmer"),
649              function(object, ...) object@fitted)
650    
651    setMethod("residuals", signature(object = "lmer"),
652              function(object, ...) object@residuals)
653    
654    setMethod("resid", signature(object = "lmer"),
655              function(object, ...) do.call("residuals", c(list(object), list(...))))
656    
657    setMethod("coef", signature(object = "lmer"),
658              function(object, ...)
659          {
660              fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
661              ref <- as(ranef(object), "list")
662              names(ref) <- names(object@flist)
663              val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
664              for (i in seq(a = val)) {
665                  refi <- ref[[i]]
666                  row.names(val[[i]]) <- row.names(refi)
667                  if (!all(names(refi) %in% names(fef)))
668                      stop("unable to align random and fixed effects")
669                  val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
670              }
671              new("lmer.coef", val)
672          })
673    
674    setMethod("plot", signature(x = "lmer.coef"),
675              function(x, y, ...)
676          {
677              varying <- unique(do.call("c",
678                                        lapply(x, function(el)
679                                               names(el)[sapply(el,
680                                                                function(col)
681                                                                any(col != col[1]))])))
682              gf <- do.call("rbind", lapply(x, "[", j = varying))
683              gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
684              switch(min(length(varying), 3),
685                     qqmath(eval(substitute(~ x | .grp,
686                                            list(x = as.name(varying[1])))), gf, ...),
687                     xyplot(eval(substitute(y ~ x | .grp,
688                                            list(y = as.name(varying[1]),
689                                                 x = as.name(varying[2])))), gf, ...),
690                     splom(~ gf | .grp, ...))
691          })
692    
693    setMethod("plot", signature(x = "lmer.ranef"),
694              function(x, y, ...)
695          {
696              lapply(x, function(x) {
697                  cn <- lapply(colnames(x), as.name)
698                  switch(min(ncol(x), 3),
699                         qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
700                         xyplot(eval(substitute(y ~ x, list(y = cn[[1]], x = cn[[2]]))),
701                                x, ...),
702                         splom(~ x, ...))
703              })
704          })
705    
706    setMethod("with", signature(data = "lmer"),
707              function(data, expr, ...)
708                  eval(substitute(substitute(expr, list(. = quote(data))),
709                                  append(data@flist, eval(data@call$data)),
710                                  enclos = parent.frame())
711              })

Legend:
Removed from v.688  
changed lines
  Added in v.689

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