SCM Repository

[matrix] Diff of /pkg/R/lmer.R
 [matrix] / pkg / R / lmer.R

Diff of /pkg/R/lmer.R

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                      controlvalsmsMaxIter <- 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