# SCM Repository

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

# Diff of /pkg/R/lmer.R

branches/trunk-lme4/R/lmer.R revision 754, Mon May 30 18:40:45 2005 UTC pkg/R/lmer.R revision 1165, Sat Jan 14 23:41:55 2006 UTC
# Line 1  Line 1
1  ## Methods for lmer and for the objects that it produces  # Methods for lmer and for the objects that it produces
2
3  ## Some utilities  ## Some utilities
4
5  Lind <- function(i,j) {  ## Return the pairs of expressions separated by vertical bars
6      if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))  findbars <- function(term)
7      ((i - 1) * i)/2 + j  {
8  }      if (is.name(term) || is.numeric(term)) return(NULL)
9        if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
10  Dhalf <- function(from) {      if (!is.call(term)) stop("term must be of class call")
11      D <- from@D      if (term[[1]] == as.name('|')) return(term)
12      nf <- length(D)      if (length(term) == 2) return(findbars(term[[2]]))
13      Gp <- from@Gp      c(findbars(term[[2]]), findbars(term[[3]]))
14      res <- array(0, rep(Gp[nf+1],2))  }
15      for (i in 1:nf) {
16          DD <- D[[i]]  ## Return the formula omitting the pairs of expressions
17          dd <- dim(DD)  ## that are separated by vertical bars
18          for (k in 1:dd[3]) {  nobars <- function(term)
19              mm <- array(DD[ , , k], dd[1:2])  {
20              base <- Gp[i] + (k - 1)*dd[1]      if (!('|' %in% all.names(term))) return(term)
21              res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)      if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
22          }      if (length(term) == 2) {
23      }          nb <- nobars(term[[2]])
24      res          if (is.null(nb)) return(NULL)
25  }          term[[2]] <- nb
26            return(term)
27  lmerControl <-                            # Control parameters for lmer      }
28    function(maxIter = 50,      nb2 <- nobars(term[[2]])
29             msMaxIter = 50,      nb3 <- nobars(term[[3]])
30             tolerance = sqrt((.Machine$double.eps)), if (is.null(nb2)) return(nb3) 31 niterEM = 15, if (is.null(nb3)) return(nb2) 32 msTol = sqrt(.Machine$double.eps),      term[[2]] <- nb2
33        term[[3]] <- nb3
34        term
35    }
36
37    ## Substitute the '+' function for the '|' function
38    subbars <- function(term)
39    {
40        if (is.name(term) || is.numeric(term)) return(term)
41        if (length(term) == 2) {
42            term[[2]] <- subbars(term[[2]])
43            return(term)
44        }
45        stopifnot(length(term) == 3)
46        if (is.call(term) && term[[1]] == as.name('|'))
47            term[[1]] <- as.name('+')
48        term[[2]] <- subbars(term[[2]])
49        term[[3]] <- subbars(term[[3]])
50        term
51    }
52
53    ## Expand an expression with colons to the sum of the lhs
54    ## and the current expression
55    colExpand <- function(term)
56    {
57        if (is.name(term) || is.numeric(term)) return(term)
58        if (length(term) == 2) {
59            term[[2]] <- colExpand(term[[2]])
60            return(term)
61        }
62        stopifnot(length(term) == 3)
63        if (is.call(term) && term[[1]] == as.name(':')) {
64            return(substitute(A+B, list(A = term, B = colExpand(term[[2]]))))
65        }
66        term[[2]] <- colExpand(term[[2]])
67        term[[3]] <- colExpand(term[[3]])
68        term
69    }
70
71    abbrvNms <- function(gnm, cnms)
72    {
73        ans <- paste(abbreviate(gnm), abbreviate(cnms), sep = '.')
74        if (length(cnms) > 1) {
75            anms <- lapply(cnms, abbreviate, minlength = 3)
76            nmmat <- outer(anms, anms, paste, sep = '.')
77            ans <- c(ans, paste(abbreviate(gnm, minlength = 3),
78                                nmmat[upper.tri(nmmat)], sep = '.'))
79        }
80        ans
81    }
82
83    ## Control parameters for lmer
84    lmerControl <-
85      function(maxIter = 200, # used in ../src/lmer.c only
86               tolerance = sqrt(.Machine$double.eps),# ditto 87 msMaxIter = 200, 88 ## msTol = sqrt(.Machine$double.eps),
89               ## FIXME:  should be able to pass tolerances to nlminb()
90             msVerbose = getOption("verbose"),             msVerbose = getOption("verbose"),
91             PQLmaxIt = 20,             niterEM = 15,
92             EMverbose = getOption("verbose"),             EMverbose = getOption("verbose"),
93               PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead
95             analyticHessian=FALSE)             analyticHessian = FALSE # unused _FIXME_
96               )
97    {
98        list(maxIter = as.integer(maxIter),
99             tolerance = as.double(tolerance),
100             msMaxIter = as.integer(msMaxIter),
101             ## msTol = as.double(msTol),
102             msVerbose = as.integer(msVerbose),# "integer" on purpose
103             niterEM = as.integer(niterEM),
104             EMverbose = as.logical(EMverbose),
105             PQLmaxIt = as.integer(PQLmaxIt),
107             analyticHessian = as.logical(analyticHessian))
108    }
109
110    setMethod("coef", signature(object = "lmer"),
111              function(object, ...)
112  {  {
113      list(maxIter = maxIter,            fef <- data.frame(rbind(object@fixed), check.names = FALSE)
114           msMaxIter = msMaxIter,            ref <- as(ranef(object), "list")
115           tolerance = tolerance,            names(ref) <- names(object@flist)
116           niterEM = niterEM,            val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
117           msTol = msTol,            for (i in seq(a = val)) {
118           msVerbose = msVerbose,                refi <- ref[[i]]
119           PQLmaxIt = PQLmaxIt,                row.names(val[[i]]) <- row.names(refi)
120           EMverbose=EMverbose,                if (!all(names(refi) %in% names(fef)))
121           analyticHessian=analyticHessian,                    stop("unable to align random and fixed effects")
123  }  }
124              new("lmer.coef", val)
125          })
126
127  setMethod("lmer", signature(formula = "formula", family = "missing"),  ## setMethod("plot", signature(x = "lmer.coef"),
128    ##           function(x, y, ...)
129    ##       {
130    ##           varying <- unique(do.call("c",
131    ##                                     lapply(x, function(el)
132    ##                                            names(el)[sapply(el,
133    ##                                                             function(col)
134    ##                                                             any(col != col[1]))])))
135    ##           gf <- do.call("rbind", lapply(x, "[", j = varying))
136    ##           gf$.grp <- factor(rep(names(x), sapply(x, nrow))) 137 ## switch(min(length(varying), 3), 138 ## qqmath(eval(substitute(~ x | .grp, 139 ## list(x = as.name(varying[1])))), gf, ...), 140 ## xyplot(eval(substitute(y ~ x | .grp, 141 ## list(y = as.name(varying[1]), 142 ## x = as.name(varying[2])))), gf, ...), 143 ## splom(~ gf | .grp, ...)) 144 ## }) 145 146 ## setMethod("plot", signature(x = "lmer.ranef"), 147 ## function(x, y, ...) 148 ## { 149 ## lapply(x, function(x) { 150 ## cn <- lapply(colnames(x), as.name) 151 ## switch(min(ncol(x), 3), 152 ## qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...), 153 ## xyplot(eval(substitute(y ~ x, 154 ## list(y = cn[[1]], 155 ## x = cn[[2]]))), x, ...), 156 ## splom(~ x, ...)) 157 ## }) 158 ## }) 159 160 setMethod("with", signature(data = "lmer"), 161 function(data, expr, ...) { 162 dat <- eval(data@call$data)
163                  if (!is.null(na.act <- attr(data@frame, "na.action")))
164                      dat <- dat[-na.act, ]
165                  lst <- c(list(. = data), data@flist, data@frame, dat)
166                  eval(substitute(expr), lst[unique(names(lst))])
167              })
168
169    setMethod("terms", signature(x = "lmer"),
170              function(x, ...) x@terms)
171
172    rWishart <- function(n, df, invScal)
173        .Call("Matrix_rWishart", n, df, invScal, PACKAGE = "Matrix")
174
175
176    setMethod("lmer", signature(formula = "formula"),
177            function(formula, data, family,            function(formula, data, family,
178                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
179                     control = list(),                     control = list(), start,
180                     subset, weights, na.action, offset,                     subset, weights, na.action, offset,
181                     model = TRUE, x = FALSE, y = FALSE, ...)                     model = TRUE, x = FALSE, y = FALSE, ...)
182        {        {
183                                          # match and check parameters            ## match and check parameters
184              if (length(formula) < 3) stop("formula must be a two-sided formula")
185              cv <- do.call("lmerControl", control)
186
187              ## Must evaluate the model frame first and then fit the glm using
188              ## that frame.  Otherwise missing values in the grouping factors
189              ## cause inconsistent numbers of observations.
190              mf <- match.call()
191              m <- match(c("family", "data", "subset", "weights",
192                           "na.action", "offset"), names(mf), 0)
193              mf <- fe <- mf[c(1, m)]
194              frame.form <- subbars(formula) # substitute +' for |'
195              fixed.form <- nobars(formula)  # remove any terms with |'
196              if (inherits(fixed.form, "name")) # RHS is empty - use a constant
197                  fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
198              environment(fixed.form) <- environment(frame.form) <- environment(formula)
199
200              ## evaluate a model frame for fixed and random effects
201              mf$formula <- frame.form 202 mf$family <- NULL
203              mf$drop.unused.levels <- TRUE 204 mf[[1]] <- as.name("model.frame") 205 frm <- eval(mf, parent.frame()) 206 207 ## fit a glm model to the fixed formula 208 fe$formula <- fixed.form
209              fe$subset <- NULL # subset has already been created in call to data.frame 210 fe$data <- frm
211              fe$x <- fe$model <- fe$y <- TRUE 212 fe[[1]] <- as.name("glm") 213 glm.fit <- eval(fe, parent.frame()) 214 x <- glm.fit$x
215              y <- as.double(glm.fit$y) 216 family <- glm.fit$family
217
218              ## check for a linear mixed model
219              lmm <- family$family == "gaussian" && family$link == "identity"
220              if (lmm) { # linear mixed model
221            method <- match.arg(method)            method <- match.arg(method)
222            if (method %in% c("PQL", "Laplace", "AGQ")) {            if (method %in% c("PQL", "Laplace", "AGQ")) {
223                warning(paste('Argument method = "', method,                warning(paste('Argument method = "', method,
# Line 63  Line 225
225                              'Using method = "REML".\n', sep = ''))                              'Using method = "REML".\n', sep = ''))
226                method <- "REML"                method <- "REML"
227            }            }
228            controlvals <- do.call("lmerControl", control)            } else { # generalized linear mixed model
229            controlvals$REML <- REML <- method == "REML" if (missing(method)) method <- "PQL" 230 else { 231 if (length(formula) < 3) stop("formula must be a two-sided formula") method <- match.arg(method) 232 if (method == "ML") method <- "PQL" 233 mf <- match.call() # create the model frame as frm if (method == "REML") 234 m <- match(c("data", "subset", "weights", "na.action", "offset"), warning('Argument method = "REML" is not meaningful ', 235 names(mf), 0) 'for a generalized linear mixed model.', 236 mf <- mf[c(1, m)] '\nUsing method = "PQL".\n') 237 mf[[1]] <- as.name("model.frame") } 238 frame.form <- subbars(formula) } 239 environment(frame.form) <- environment(formula) ## create factor list for the random effects mf$formula <- frame.form
mf$drop.unused.levels <- TRUE frm <- eval(mf, parent.frame()) ## grouping factors and model matrices for random effects 240 bars <- findbars(formula[[3]]) bars <- findbars(formula[[3]]) 241 random <- names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]]))) 242 lapply(bars, fl <- lapply(bars, 243 function(x) list(model.matrix(eval(substitute(~term, function(x) list(term=x[[2]]))), frm), 244 eval(substitute(as.factor(fac)[,drop = TRUE], eval(substitute(as.factor(fac)[,drop = TRUE], 245 list(fac = x[[3]])), frm))) list(fac = x[[3]])), frm)) names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]]))) 246 ## order factor list by decreasing number of levels ## order factor list by decreasing number of levels 247 nlev <- sapply(random, function(x) length(levels(x[[2]]))) nlev <- sapply(fl, function(x) length(levels(x))) 248 if (any(diff(nlev) > 0)) { if (any(diff(nlev) > 0)) { 249 random <- random[rev(order(nlev))] ord <- rev(order(nlev)) 250 } bars <- bars[ord] 251 fixed.form <- nobars(formula) fl <- fl[ord] 252 if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula } 253 Xmat <- model.matrix(fixed.form, frm) ## create list of transposed model matrices for random effects 254 mmats <- c(lapply(random, "[[", 1), Ztl <- lapply(bars, function(x) 255 .fixed = list(cbind(Xmat, .response = model.response(frm)))) t(model.matrix(eval(substitute(~ expr, 256 obj <- .Call("lmer_create", lapply(random, "[[", 2), list(expr = x[[2]]))), 257 mmats, PACKAGE = "Matrix") frm))) 258 slot(obj, "frame") <- frm ## Create the mixed-effects representation (mer) object 259 slot(obj, "terms") <- attr(model.frame(fixed.form, data), "terms") mer <- .Call("mer_create", fl, 260 slot(obj, "assign") <- attr(Xmat, "assign") .Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"), 261 slot(obj, "call") <- match.call() x, y, method, sapply(Ztl, nrow), 262 slot(obj, "REML") <- REML c(lapply(Ztl, rownames), list(.fixed = colnames(x))), 263 rm(Xmat) !(family$family %in% c("binomial", "poisson")),
264            .Call("lmer_initial", obj, PACKAGE="Matrix")                         match.call(), family,
.Call("lmer_ECMEsteps", obj,
controlvals$niterEM, controlvals$REML,
controlvals$EMverbose, 265 PACKAGE = "Matrix") PACKAGE = "Matrix") 266 LMEoptimize(obj) <- controlvals if (lmm) { 267 slot(obj, "residuals") <- .Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix") 268 unname(model.response(frm) - LMEoptimize(mer) <- cv 269 (slot(obj, "fitted") <- return(mer) 270 .Call("lmer_fitted", obj, mmats, TRUE, PACKAGE = "Matrix"))) } 271 obj 272 ## The rest of the function applies to generalized linear mixed models 273 gVerb <- getOption("verbose") 274 eta <- glm.fit$linear.predictors
275              wts <- glm.fit$prior.weights 276 wtssqr <- wts * wts 277 offset <- glm.fit$offset
278              if (is.null(offset)) offset <- numeric(length(eta))
279              linkinv <- quote(family$linkinv(eta)) 280 mu.eta <- quote(family$mu.eta(eta))
281              mu <- family$linkinv(eta) 282 variance <- quote(family$variance(mu))
283              dev.resids <- quote(family$dev.resids(y, mu, wtssqr)) 284 LMEopt <- get("LMEoptimize<-") 285 doLMEopt <- quote(LMEopt(x = mer, value = cv)) 286 287 GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix") 288 .Call("glmer_PQL", GSpt, PACKAGE = "Matrix") # obtain PQL estimates 289 fixInd <- seq(ncol(x)) 290 ## pars[fixInd] == beta, pars[-fixInd] == theta 291 PQLpars <- c(fixef(mer), 292 .Call("mer_coef", mer, 2, PACKAGE = "Matrix")) 293 .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix") 294 ## indicator of constrained parameters 295 const <- c(rep(FALSE, length(fixInd)), 296 unlist(lapply(mer@nc[seq(along = fl)], 297 function(k) 1:((k*(k+1))/2) <= k) 298 )) 299 devLaplace <- function(pars) 300 .Call("glmer_devLaplace", pars, GSpt, PACKAGE = "Matrix") 301 302 optimRes <- 303 nlminb(PQLpars, devLaplace, 304 lower = ifelse(const, 5e-10, -Inf), 305 control = list(trace = getOption("verbose"), 306 iter.max = cv$msMaxIter))
307              .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
308              return(mer)
309              deviance <- devAGQ(PQLpars, 1)
310
311    ### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs
312    ### AGQ for nc > 1 first.
313              fxd <- PQLpars[fixInd]
314              loglik <- logLik(mer)
315
316              if (method %in% c("Laplace", "AGQ")) {
317                  nAGQ <- 1
318                  if (method == "AGQ") {    # determine nAGQ at PQL estimates
319                      dev11 <- devAGQ(PQLpars, 11)
320                      ## FIXME: Should this be an absolute or a relative tolerance?
321                      devTol <- sqrt(.Machine$double.eps) * abs(dev11) 322 for (nAGQ in c(9, 7, 5, 3, 1)) 323 if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break 324 nAGQ <- nAGQ + 2 325 if (gVerb) 326 cat(paste("Using", nAGQ, "quadrature points per column\n")) 327 } 328 obj <- function(pars) 329 .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix") 330 optimRes <- 331 nlminb(PQLpars, obj, 332 lower = ifelse(const, 5e-10, -Inf), 333 control = list(trace = getOption("verbose"), 334 iter.max = cv$msMaxIter))
335                  optpars <- optimRes$par 336 if (optimRes$convergence != 0)
337                      warning("nlminb failed to converge")
338                  deviance <- optimRes$objective 339 if (gVerb) 340 cat(paste("convergence message", optimRes$message, "\n"))
341                  fxd[] <- optpars[fixInd]  ## preserve the names
342                  .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
343              }
344
345              .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
346              loglik[] <- -deviance/2
347              new("lmer", mer,
348                  frame = if (model) frm else data.frame(),
349                  terms = glm.fit$terms, 350 assign = attr(glm.fit$x, "assign"),
351                  call = match.call(), family = family,
352                  logLik = loglik, fixed = fxd)
353        })        })
354
355  setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),
356    ## Extract the permutation
357    setAs("mer", "pMatrix", function(from)
358          .Call("mer_pMatrix", from, PACKAGE = "Matrix"))
359
360    ## Extract the L matrix
361    setAs("mer", "dtCMatrix", function(from)
362          .Call("mer_dtCMatrix", from, PACKAGE = "Matrix"))
363
364    ## Extract the fixed effects
365    setMethod("fixef", signature(object = "mer"),
366              function(object, ...)
367              .Call("mer_fixef", object, PACKAGE = "Matrix"))
368
369    ## Extract the random effects
370    setMethod("ranef", signature(object = "mer"),
371              function(object, ...)
372                  .Call("mer_ranef", object, PACKAGE = "Matrix")
373              )
374
375    ## Optimization for mer objects
376    setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
377                   function(x, value)                   function(x, value)
378               {               {
379                   if (value$msMaxIter < 1) return(x) if (value$msMaxIter < 1) return(x)
380                   nc <- x@nc                   nc <- x@nc
nc <- nc[1:(length(nc) - 2)]
381                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
382                   fn <- function(pars)                   fn <- function(pars)
383                       deviance(.Call("lmer_coefGets", x, pars,                       deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
2, PACKAGE = "Matrix"),
REML = value$REML) 384 gr <- if (value$analyticGradient)                   gr <- if (value$analyticGradient) 385 function(pars) { function(pars) { 386 if (!isTRUE(all.equal(pars, if (!isTRUE(all.equal(pars, 387 .Call("lmer_coef", x, 2)))) .Call("mer_coef", x, 388 .Call("lmer_coefGets", x, pars, 2) 2, PACKAGE = "Matrix")))) 389 .Call("lmer_gradient", x, value$REML, 2)                               .Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")
390                       } else NULL                           .Call("mer_gradient", x, 2, PACKAGE = "Matrix")
391                   RV <- lapply(R.Version()[c("major", "minor")], as.numeric)                       }
392                   if ((RV$major > 1 && RV$minor >= 2.0) ||                   else NULL
393                       require("port", quietly = TRUE)) {                   optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"),
optimRes <- nlminb(.Call("lmer_coef", x, 2),
394                                          fn, gr,                                          fn, gr,
395                                          lower = ifelse(constr, 1e-10, -Inf),                                      lower = ifelse(constr, 5e-10, -Inf),
396                                          control = list(iter.max = value$msMaxIter, control = list(iter.max = value$msMaxIter,
397                                          trace = as.integer(value$msVerbose))) trace = as.integer(value$msVerbose)))
398                   } else {                   .Call("mer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix") optimRes <- optim(.Call("lmer_coef", x, 2), fn, gr, method = "L-BFGS-B", lower = ifelse(constr, 1e-10, -Inf), control = list(maxit = value$msMaxIter,
trace = as.integer(value$msVerbose))) } .Call("lmer_coefGets", x, optimRes$par, 2)
399                   if (optimRes$convergence != 0) { if (optimRes$convergence != 0) {
400                       warning(paste("optim returned message",optimRes$message,"\n")) warning(paste("nlminb returned message", 401 optimRes$message,"\n"))
402                   }                   }
403                   return(x)                   return(x)
404               })               })
405
406  setMethod("ranef", signature(object = "lmer"),  setMethod("deviance", signature(object = "mer"),
function(object, accumulate = FALSE, ...) {
val <- new("lmer.ranef",
lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),
data.frame, check.names = FALSE),
varFac = object@bVar,
stdErr = .Call("lmer_sigma", object,
object@REML, PACKAGE = "Matrix"))
if (!accumulate || length(val@varFac) == 1) return(val)
## check for nested factors
L <- object@L
if (any(sapply(seq(a = val), function(i) length(L[[Lind(i,i)]]@i))))
error("Require nested grouping factors to accumulate random effects")
val
})

setMethod("fixef", signature(object = "lmer"),
function(object, ...) {
val <- .Call("lmer_fixef", object, PACKAGE = "Matrix")
val[-length(val)]
})

setMethod("fixef", signature(object = "glmer"),
407            function(object, ...) {            function(object, ...) {
408                object@fixed                .Call("mer_factor", object, PACKAGE = "Matrix")
409                  object@deviance[[ifelse(object@method == "REML", "REML", "ML")]]
410            })            })
411
412  setMethod("VarCorr", signature(x = "lmer"),  setMethod("mcmcsamp", signature(object = "mer"),
413            function(x, REML = TRUE, useScale = TRUE, ...) {            function(object, n = 1, verbose = FALSE, saveb = FALSE,
414                val <- .Call("lmer_variances", x, PACKAGE = "Matrix")                     trans = TRUE, ...)
415                for (i in seq(along = val)) {        {
416                    dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])            ans <- t(.Call("mer_MCMCsamp", object, saveb, n,
417                    val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")                           trans, PACKAGE = "Matrix"))
418                }            attr(ans, "mcpar") <- as.integer(c(1, n, 1))
419                new("VarCorr",            class(ans) <- "mcmc"
420                    scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),            glmer <- FALSE
421                    reSumry = val,            gnms <- names(object@flist)
422                    useScale = useScale)            cnms <- object@cnames
423              ff <- fixef(object)
424              colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
425                          unlist(lapply(seq(along = gnms),
426                                        function(i)
427                                        abbrvNms(gnms[i],cnms[[i]]))))
428              if (trans) {
429                  ## parameter type: 0 => fixed effect, 1 => variance,
430                  ##                 2 => covariance
431                  ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1,
432                            unlist(lapply(seq(along = gnms),
433                                          function(i)
434                                      {
435                                          k <- length(cnms[[i]])
436                                          rep(1:2, c(k, (k*(k-1))/2))
437                                      })))
438                  colnms[ptyp == 1] <-
439                      paste("log(", colnms[ptyp == 1], ")", sep = "")
440                  colnms[ptyp == 2] <-
441                      paste("atanh(", colnms[ptyp == 2], ")", sep = "")
442              }
443              colnames(ans) <- colnms
444              ans
445          })
446
447    setMethod("simulate", signature(object = "mer"),
448              function(object, nsim = 1, seed = NULL, ...)
449          {
450              if(!exists(".Random.seed", envir = .GlobalEnv))
451                  runif(1)               # initialize the RNG if necessary
452              if(is.null(seed))
453                  RNGstate <- .Random.seed
454              else {
455                  R.seed <- .Random.seed
456                  set.seed(seed)
457                  RNGstate <- structure(seed, kind = as.list(RNGkind()))
458                  on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
459              }
460
461              family <- object@family
462              if (family$family != "gaussian" || 463 family$link != "identity")
464                  stop("simulation of generalized linear mixed models not yet implemented")
465              ## similate the linear predictors
466              lpred <- .Call("mer_simulate", object, nsim, PACKAGE = "Matrix")
467              sc <- 1
468              if (object@useScale)
469                  sc <- .Call("mer_sigma", object, object@method == "REML",
470                              PACKAGE = "Matrix")
471              ## add fixed-effects contribution and per-observation noise term
472              lpred <- as.data.frame(lpred + drop(object@X %*% fixef(object)) +
473                                     rnorm(prod(dim(lpred)), sd = sc))
474              ## save the seed
475              attr(lpred, "seed") <- RNGstate
476              lpred
477            })            })
478
function(x, REML, unconst, ...)
.Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix"))

setMethod("summary", signature(object = "lmer"),
function(object, ...)
new("summary.lmer", object, useScale = TRUE,
showCorrelation = TRUE,
method = if (object@REML) "REML" else "ML",
family = gaussian(),
logLik = logLik(object),
fixed = fixef(object)))

## FIXME: glmm-s with scale not handled
setMethod("summary", signature(object = "glmer"),
function(object, ...)
new("summary.lmer", object, useScale = FALSE,
showCorrelation = TRUE,
method = object@method,
family = object@family,
logLik = logLik(object),
fixed = fixef(object)))

setMethod("show", signature(object = "lmer"),
function(object)
show(new("summary.lmer", object, useScale = TRUE,
showCorrelation = FALSE,
method = if (object@REML) "REML" else "ML",
family = gaussian(),
logLik = logLik(object),
fixed = fixef(object)))
)

setMethod("show", signature(object = "glmer"),
function(object)
show(new("summary.lmer", object, useScale = FALSE,
showCorrelation = FALSE,
method = object@method,
family = object@family,
logLik = logLik(object),
fixed = fixef(object)))
)
479
480  setMethod("show", "summary.lmer",  setMethod("show", "mer",
481            function(object) {            function(object) {
482                fcoef <- object@fixed                vcShow <- function(varc, useScale)
483                  {
484                      digits <- max(3, getOption("digits") - 2)
485                      sc <- attr(varc, "sc")
486                      recorr <- lapply(varc, function(el) el@factors$correlation) 487 reStdDev <- c(lapply(recorr, slot, "sd"), list(Residual = sc)) 488 reLens <- unlist(c(lapply(reStdDev, length))) 489 reMat <- array('', c(sum(reLens), 4), 490 list(rep('', sum(reLens)), 491 c("Groups", "Name", "Variance", "Std.Dev."))) 492 reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens) 493 reMat[,2] <- c(unlist(lapply(reStdDev, names)), "") 494 reMat[,3] <- format(unlist(reStdDev)^2, digits = digits) 495 reMat[,4] <- format(unlist(reStdDev), digits = digits) 496 if (any(reLens > 1)) { 497 maxlen <- max(reLens) 498 corr <- 499 do.call("rbind", 500 lapply(recorr, 501 function(x, maxlen) { 502 x <- as(x, "matrix") 503 cc <- format(round(x, 3), nsmall = 3) 504 cc[!lower.tri(cc)] <- "" 505 nr <- dim(cc)[1] 506 if (nr >= maxlen) return(cc) 507 cbind(cc, matrix("", nr, maxlen-nr)) 508 }, maxlen)) 509 colnames(corr) <- c("Corr", rep("", maxlen - 1)) 510 reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr)))) 511 } 512 if (!useScale) reMat <- reMat[-nrow(reMat),] 513 print(reMat, quote = FALSE) 514 } 515 516 fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix") 517 useScale <- object@useScale useScale <- object@useScale 518 corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"), corF <- vcov(object)@factors$correlation
"corrmatrix")
519                DF <- getFixDF(object)                DF <- getFixDF(object)
520                coefs <- cbind(fcoef, corF@stdDev, DF)                coefs <- cbind(fcoef, corF@sd, DF)
nc <- object@nc
521                dimnames(coefs) <-                dimnames(coefs) <-
522                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))
523                              digits <- max(3, getOption("digits") - 2)                              digits <- max(3, getOption("digits") - 2)
524                REML <- length(object@REML) > 0 && object@REML[1]                REML <- object@method == "REML"
525                llik <- object@logLik                llik <- logLik(object, REML)
526                dev <- object@deviance                dev <- object@deviance
527                  devc <- object@devComp
528
529                rdig <- 5                rdig <- 5
530                if (glz <- !(object@method %in% c("REML", "ML"))) {                if (glz <- !(object@method %in% c("REML", "ML"))) {
# Line 267  Line 532
532                              object@method, "\n"))                              object@method, "\n"))
533                } else {                } else {
534                    cat("Linear mixed-effects model fit by ")                    cat("Linear mixed-effects model fit by ")
535                    cat(if(object@REML) "REML\n" else "maximum likelihood\n")                    cat(if(REML) "REML\n" else "maximum likelihood\n")
536                }                }
537                if (!is.null(object@call$formula)) { if (!is.null(object@call$formula)) {
538                    cat("Formula:", deparse(object@call$formula),"\n") cat("Formula:", deparse(object@call$formula),"\n")
# Line 294  Line 559
559                                 row.names = ""))                                 row.names = ""))
560                }                }
561                cat("Random effects:\n")                cat("Random effects:\n")
562                show(VarCorr(object))                vcShow(VarCorr(object), useScale)
563                ngrps <- lapply(object@flist, function(x) length(levels(x)))                ngrps <- lapply(object@flist, function(x) length(levels(x)))
564                cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))                cat(sprintf("# of obs: %d, groups: ", devc[1]))
565                cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))                cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
566                cat("\n")                cat("\n")
567                if (!useScale)                if (!useScale)
568                    cat("\nEstimated scale (compare to 1) ",                    cat("\nEstimated scale (compare to 1) ",
569                        .Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"),                        .Call("mer_sigma", object, FALSE, PACKAGE = "Matrix"),
570                        "\n")                        "\n")
571                if (nrow(coefs) > 0) {                if (nrow(coefs) > 0) {
572                    if (useScale) {                    if (useScale) {
# Line 320  Line 585
585                    }                    }
586                    cat("\nFixed effects:\n")                    cat("\nFixed effects:\n")
587                    printCoefmat(coefs, tst.ind = 4, zap.ind = 3)                    printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
588                        rn <- rownames(coefs)                        rn <- rownames(coefs)
dimnames(corF) <- list(
abbreviate(rn, minlen=11),
abbreviate(rn, minlen=6))
589                        if (!is.null(corF)) {                        if (!is.null(corF)) {
590                            p <- NCOL(corF)                        p <- ncol(corF)
591                            if (p > 1) {                            if (p > 1) {
592                                cat("\nCorrelation of Fixed Effects:\n")                                cat("\nCorrelation of Fixed Effects:\n")
593                                corF <- format(round(corF, 3), nsmall = 3)                            corF <- matrix(format(round(corF@x, 3), nsmall = 3),
594                                             nc = p)
595                              dimnames(corF) <- list(
596                                                     abbreviate(rn, minlen=11),
597                                                     abbreviate(rn, minlen=6))
598                                corF[!lower.tri(corF)] <- ""                                corF[!lower.tri(corF)] <- ""
599                                print(corF[-1, -p, drop=FALSE], quote = FALSE)                                print(corF[-1, -p, drop=FALSE], quote = FALSE)
600                            }                            }
601                        }                        }
602                    }                    }
}
603                invisible(object)                invisible(object)
604            })            })
605
606    setMethod("vcov", signature(object = "mer"),
607  setMethod("lmer", signature(formula = "formula"),            function(object, REML = object@method == "REML",
608            function(formula, family, data,                     useScale = object@useScale,...) {
609                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),                sc <- if (object@useScale) {
610                     control = list(),                    .Call("mer_sigma", object, REML, PACKAGE = "Matrix")
611                     subset, weights, na.action, offset,                } else { 1 }
612                     model = TRUE, x = FALSE, y = FALSE, ...)                rr <- as(sc^2 * tcrossprod(solve(object@RXX)), "dpoMatrix")
613        {                rr@factors$correlation <- as(rr, "correlation") 614 if (missing(method)) { rr method <- "PQL" } else { method <- 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')) } ## if (method %in% c("Laplace", "AGQ")) ## stop("Laplace and AGQ methods not yet implemented") if (method %in% c("AGQ")) stop("AGQ method not yet implemented") gVerb <- getOption("verbose") # match and check parameters controlvals <- do.call("lmerControl", control) controlvals$REML <- FALSE
if (length(formula) < 3) stop("formula must be a two-sided formula")
## initial glm fit
mf <- match.call()
m <- match(c("family", "data", "subset", "weights",
"na.action", "offset"),
names(mf), 0)
mf <- mf[c(1, m)]
mf[[1]] <- as.name("glm")
fixed.form <- nobars(formula)
if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula
environment(fixed.form) <- environment(formula)
mf$formula <- fixed.form mf$x <- mf$model <- mf$y <- TRUE
glm.fit <- eval(mf, parent.frame())
family <- glm.fit$family ## Note: offset is on the linear scale offset <- glm.fit$offset
if (is.null(offset)) offset <- 0
weights <- sqrt(abs(glm.fit$prior.weights)) ## initial 'fitted' values on linear scale etaold <- eta <- glm.fit$linear.predictors

## evaluation of model frame
mf$x <- mf$model <- mf$y <- mf$family <- NULL
mf$drop.unused.levels <- TRUE this.form <- subbars(formula) environment(this.form) <- environment(formula) mf$formula <- this.form
mf[[1]] <- as.name("model.frame")
frm <- eval(mf, parent.frame())

## grouping factors and model matrices for random effects
bars <- findbars(formula[[3]])
random <-
lapply(bars,
function(x) list(model.matrix(eval(substitute(~term,
list(term=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))]
}
mmats <- c(lapply(random, "[[", 1),
.fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
## FIXME: Use Xfrm and Xmat to get the terms and assign
## slots, pass these to lmer_create, then destroy Xfrm, Xmat, etc.
obj <- .Call("lmer_create", lapply(random, "[[", 2),
mmats, PACKAGE = "Matrix")
slot(obj, "frame") <- frm
slot(obj, "terms") <- attr(glm.fit$model, "terms") slot(obj, "assign") <- attr(glm.fit$x, "assign")
slot(obj, "call") <- match.call()
slot(obj, "REML") <- FALSE
rm(glm.fit)
.Call("lmer_initial", obj, PACKAGE="Matrix")
mmats[[1]][1,1] <- mmats[[1]][1,1]
conv <- FALSE
firstIter <- TRUE
msMaxIter.orig <- controlvals$msMaxIter responseIndex <- ncol(mmats$.fixed)

for (iter in seq(length = controlvals$PQLmaxIt)) { mu <- family$linkinv(eta)
dmu.deta <- family$mu.eta(eta) ## weights (note: weights is already square-rooted) w <- weights * dmu.deta / sqrt(family$variance(mu))
## adjusted response (should be comparable to X \beta, not including offset
z <- eta - offset + (mmats.unadjusted$.fixed[, responseIndex] - mu) / dmu.deta .Call("nlme_weight_matrix_list", mmats.unadjusted, w, z, mmats, PACKAGE="Matrix") .Call("lmer_update_mm", obj, mmats, PACKAGE="Matrix") if (firstIter) { .Call("lmer_initial", obj, PACKAGE="Matrix") if (gVerb) cat(" PQL iterations convergence criterion\n") } .Call("lmer_ECMEsteps", obj, controlvals$niterEM,
FALSE,
controlvals$EMverbose, PACKAGE = "Matrix") LMEoptimize(obj) <- controlvals eta[] <- offset + ## FIXME: should the offset be here ? .Call("lmer_fitted", obj, mmats.unadjusted, TRUE, PACKAGE = "Matrix") crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta))) if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit)) ## use this to determine convergence if (crit < controlvals$tolerance) {
conv <- TRUE
break
}
etaold[] <- eta

## Changing number of iterations on second and
## subsequent iterations.
if (firstIter)
{
controlvals$niterEM <- 2 controlvals$msMaxIter <- 10
firstIter <- FALSE
}
}
if (!conv) warning("IRLS iterations for glmm did not converge")
controlvals$msMaxIter <- msMaxIter.orig ### if (TRUE) ## Laplace ### { ## Need to optimize L(theta, beta) using Laplace approximation ## Things needed for that: ## ## 1. reduced ssclme object, offset, weighted model matrices ## 2. facs, reduced model matrices ## Of these, those in 2 will be fixed given theta and beta, ## and can be thought of arguments to the L(theta, beta) ## function. However, the ones in 1 will have the same ## structure. So the plan is to pre-allocate them and pass ## them in too so they can be used without creating/copying ## them more than once ## reduced ssclme reducedObj <- .Call("lmer_collapse", obj, PACKAGE = "Matrix") reducedMmats.unadjusted <- mmats.unadjusted reducedMmats.unadjusted$.fixed <-
reducedMmats.unadjusted$.fixed[, responseIndex, drop = FALSE] reducedMmats <- mmats reducedMmats$.fixed <-
reducedMmats$.fixed[, responseIndex, drop = FALSE] ## define function that calculates bhats given theta and beta bhat <- function(pars = NULL) # 1:(responseIndex-1) - beta, rest - theta { if (is.null(pars)) { off <- drop(mmats.unadjusted$.fixed %*%
c(fixef(obj), 0)) + offset
}
else
{
.Call("lmer_coefGets",
reducedObj,
as.double(pars[responseIndex:length(pars)]),
TRUE,
PACKAGE = "Matrix")
off <- drop(mmats.unadjusted$.fixed %*% c(pars[1:(responseIndex-1)], 0)) + offset } niter <- 20 conv <- FALSE eta <- offset + .Call("lmer_fitted", obj, mmats.unadjusted, TRUE, PACKAGE = "Matrix") etaold <- eta for (iter in seq(length = niter)) { mu <- family$linkinv(eta)
dmu.deta <- family$mu.eta(eta) w <- weights * dmu.deta / sqrt(family$variance(mu))
z <- eta - off + (reducedMmats.unadjusted$.fixed[, 1] - mu) / dmu.deta .Call("nlme_weight_matrix_list", reducedMmats.unadjusted, w, z, reducedMmats, PACKAGE="Matrix") .Call("lmer_update_mm", reducedObj, reducedMmats, PACKAGE="Matrix") eta[] <- off + .Call("lmer_fitted", reducedObj, reducedMmats.unadjusted, TRUE, PACKAGE = "Matrix") ##cat(paste("bhat Criterion:", max(abs(eta - etaold)) / ## (0.1 + max(abs(eta))), "\n")) ## use this to determine convergence if (max(abs(eta - etaold)) < (0.1 + max(abs(eta))) * controlvals$tolerance)
{
conv <- TRUE
break
}
etaold[] <- eta

}
if (!conv) warning("iterations for bhat did not converge")

## bhat doesn't really need to return anything, we
## just want the side-effect of modifying reducedObj
## In particular, we are interested in
## ranef(reducedObj) and reducedObj@bVar (?). But
## the mu-scale response will be useful for log-lik
## calculations later, so return them anyway

invisible(family$linkinv(eta)) } ## function that calculates log likelihood (the thing that ## needs to get evaluated at each Gauss-Hermite location) ## log scale ? worry about details later, get the pieces in place ## this is for the Laplace approximation only. GH is more ## complicated devLaplace <- function(pars = NULL) { ## FIXME: This actually returns half the deviance. ## gets correct values of bhat and bvars. As a side ## effect, mu now has fitted values mu <- bhat(pars = pars) ## GLM family log likelihood (upto constant ?)(log scale) ## FIXME: need to adjust for sigma^2 for appropriate models (not trivial!) ## Keep everything on (log) likelihood scale ## log lik from observations given fixed and random effects ## get deviance, then multiply by -1/2 (since deviance = -2 log lik) ans <- -sum(family$dev.resids(y = mmats.unadjusted$.fixed[, responseIndex], mu = mu, wt = weights^2))/2 ## if (is.null(getOption("laplaceinR"))) ## { ans <- ans + .Call("lmer_laplace_devComp", reducedObj, PACKAGE = "Matrix") ## } ## else ## { ## ranefs <- .Call("lmer_ranef", reducedObj, PACKAGE = "Matrix") ## ## ans <- ans + reducedObj@devComp[2]/2 # log-determinant of Omega ## Omega <- reducedObj@Omega ## for (i in seq(along = ranefs)) ## { ## ## contribution for random effects (get it working, ## ## optimize later) ## ## symmetrize RE variance ## Omega[[i]] <- Omega[[i]] + t(Omega[[i]]) ## diag(Omega[[i]]) <- diag(Omega[[i]]) / 2 ## ## want log of const det(Omega) exp(-1/2 b' ## ## Omega b )` i.e., const + log det(Omega) - .5 ## ## * (b' Omega b) ## ## FIXME: need to adjust for sigma^2 for appropriate ## ## models (easy). These are all the b'Omega b, ## ## summed as they eventually need to be. Think of ## ## this as sum(rowSums((ranefs[[i]] %*% Omega[[i]]) ## ## * ranefs[[i]])) ## ranef.loglik.det <- nrow(ranefs[[i]]) * ## determinant(Omega[[i]], logarithm = TRUE)$modulus/2

## #                      print(ranef.loglik.det)

##                       ranef.loglik.re <-
##                           -sum((ranefs[[i]] %*% Omega[[i]]) * ranefs[[i]])/2

## #                      print(ranef.loglik.re)

##                       ranef.loglik <- ranef.loglik.det + ranef.loglik.re

##                       log.jacobian <-
##                           sum(log(abs(apply(reducedObj@bVar[[i]],
##                                             3,

##                                             ## next line depends on
##                                             ## whether bVars are variances
##                                             ## or Cholesly factors

##                                             ## function(x) sum(diag(x)))
## ## Was this a bug?                          function(x) sum(diag( La.chol( x ) )))
##                                             function(x) prod(diag( La.chol( x ) )))
##                                       )))

## #                      print(log.jacobian)

##                       ## the constant terms from the r.e. and the final
##                       ## Laplacian integral cancel out both being:
##                       ## ranef.loglik.constant <- 0.5 * length(ranefs[[i]]) * log(2 * base::pi)

##                       ans <- ans + ranef.loglik + log.jacobian
##                   }
##               }
## ans is (up to some constant) log of the Laplacian
## approximation of the likelihood. Return it's negative
## to be minimized

##              cat("Parameters: ")
##              print(pars)

##              cat("Value: ")
##              print(as.double(-ans))

-ans
}

if (method == "Laplace")
{
RV <- lapply(R.Version()[c("major", "minor")], as.numeric)
if ((RV$major > 1 && RV$minor >= 2.0) ||
require("port", quietly = TRUE)) {
optimRes <-
nlminb(c(fixef(obj),
.Call("lmer_coef", obj, TRUE, PACKAGE = "Matrix")),
devLaplace,
control = list(trace = getOption("verbose"),
iter.max = controlvals$msMaxIter)) optpars <- optimRes$par
if (optimRes$convergence != 0) warning("nlminb failed to converge") } else { optimRes <- optim(c(fixef(obj), .Call("lmer_coef", obj, TRUE, PACKAGE = "Matrix")), devLaplace, method = "BFGS", hessian = TRUE, control = list(trace = getOption("verbose"), reltol = controlvals$msTol,
maxit = controlvals$msMaxIter)) optpars <- optimRes$par
if (optimRes$convergence != 0) warning("optim failed to converge") Hessian <- optimRes$hessian
}

##fixef(obj) <- optimRes$par[seq(length = responseIndex - 1)] if (getOption("verbose")) { cat(paste("optim convergence code", optimRes$convergence, "\n"))
cat("Fixed effects:\n")
print(fixef(obj))
print(optimRes$par[seq(length = responseIndex - 1)]) cat("(Unconstrained) variance coefficients:\n") print( .Call("lmer_coef", obj, TRUE, PACKAGE = "Matrix")) ##coef(obj, unconst = TRUE) <- ## optimRes$par[responseIndex:length(optimRes$par)] ##print(coef(obj, unconst = TRUE)) print( optimRes$par[responseIndex:length(optimRes$par)] ) } ## need to calculate likelihood. also need to store ## new estimates of fixed effects somewhere ## (probably cannot update standard errors) ###Rprof(NULL) } else { optpars <- c(fixef(obj), .Call("lmer_coef", obj, TRUE, PACKAGE = "Matrix")) Hessian <- new("matrix") } ## Before finishing, we need to call devLaplace with the ## optimum pars to get the final log likelihood (still need ## to make sure it's the actual likelihood and not a ## multiple). This would automatically call bhat() and hence ## have the 'correct' random effects in reducedObj. loglik <- -devLaplace(optpars) ##print(loglik) ff <- optpars[1:(responseIndex-1)] names(ff) <- names(fixef(obj)) if (!x) mmats <- list() ### } ##obj new("glmer", obj, family = family, glmmll = loglik, method = method, fixed = ff) 615 }) }) 616 617 ## calculates degrees of freedom for fixed effects Wald tests ## calculates degrees of freedom for fixed effects Wald tests 618 ## This is a placeholder. The answers are generally wrong. It will ## This is a placeholder. The answers are generally wrong. It will 619 ## be very tricky to decide what a 'right' answer should be with ## be very tricky to decide what a 'right' answer should be with 620 ## crossed random effects. ## crossed random effects. 621 622 setMethod("getFixDF", signature(object="lmer"), setMethod("getFixDF", signature(object="mer"), 623 function(object, ...) function(object, ...) { 624 { devc <- object@devComp 625 nc <- object@nc[-seq(along = object@Omega)] rep(as.integer(devc[1]- devc[2]), devc[2]) p <- nc[1] - 1 n <- nc[2] rep(n-p, p) 626 }) }) 627 628 setMethod("logLik", signature(object="lmer"), setMethod("logLik", signature(object="mer"), 629 function(object, REML = object@REML, ...) { function(object, REML = object@method == "REML", ...) { 630 val <- -deviance(object, REML = REML)/2 val <- -deviance(object, REML = REML)/2 631 nc <- object@nc[-seq(a = object@Omega)] devc <- as.integer(object@devComp[1:2]) 632 attr(val, "nall") <- attr(val, "nobs") <- nc[2] attr(val, "nall") <- attr(val, "nobs") <- devc[1] 633 attr(val, "df") <- attr(val, "df") <- abs(devc[2]) + 634 nc[1] + length(.Call("lmer_coef", object, 2)) length(.Call("mer_coef", object, 0, PACKAGE = "Matrix")) 635 attr(val, "REML") <- REML attr(val, "REML") <- REML 636 class(val) <- "logLik" class(val) <- "logLik" 637 val val 638 }) }) 639 640 setMethod("logLik", signature(object="glmer"), setMethod("VarCorr", signature(x = "mer"), 641 function(object, ...) { function(x, REML = x@method == "REML", useScale = x@useScale, ...) 642 val <- object@glmmll { 643 nc <- object@nc[-seq(a = object@Omega)] sc <- 1 644 attr(val, "nall") <- attr(val, "nobs") <- nc[2] if (useScale) 645 attr(val, "df") <- sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix") 646 nc[1] + length(.Call("lmer_coef", object, 2)) sc2 <- sc * sc 647 class(val) <- "logLik" ans <- lapply(x@Omega, function(el) { 648 val el <- as(sc2 * solve(el), "dpoMatrix") 649 el@factors$correlation <- as(el, "correlation")
650                  el
651              })
652              attr(ans, "sc") <- sc
653              ans
654            })            })
655
656  setMethod("anova", signature(object = "lmer"),  setMethod("anova", signature(object = "mer"),
657            function(object, ...)            function(object, ...)
658        {        {
659            mCall <- match.call(expand.dots = TRUE)            mCall <- match.call(expand.dots = TRUE)
660            dots <- list(...)            dots <- list(...)
661            modp <- logical(0)            modp <- logical(0)
662            if (length(dots))            if (length(dots))
663                modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm")                modp <- sapply(dots, inherits, "mer") | sapply(dots, inherits, "lm")
664            if (any(modp)) {              # multiple models - form table            if (any(modp)) {              # multiple models - form table
665                opts <- dots[!modp]                opts <- dots[!modp]
666                mods <- c(list(object), dots[modp])                mods <- c(list(object), dots[modp])
# Line 855  Line 698
698            } else {            } else {
699                foo <- object                foo <- object
700                foo@status["factored"] <- FALSE                foo@status["factored"] <- FALSE
701                .Call("lmer_factor", foo, PACKAGE="Matrix")                .Call("mer_factor", foo, PACKAGE="Matrix")
702                dfr <- getFixDF(foo)                dfr <- getFixDF(foo)
703                rcol <- ncol(foo@RXX)                ss <- foo@rXy^2
704                ss <- foo@RXX[ , rcol]^2                ssr <- exp(foo@devComp["logryy2"])
ssr <- ss[[rcol]]
705                ss <- ss[seq(along = dfr)]                ss <- ss[seq(along = dfr)]
706                names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]                names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
707                asgn <- foo@assign                asgn <- foo@assign
# Line 869  Line 711
711                    nmeffects <- c("(Intercept)", nmeffects)                    nmeffects <- c("(Intercept)", nmeffects)
712                ss <- unlist(lapply(split(ss, asgn), sum))                ss <- unlist(lapply(split(ss, asgn), sum))
713                df <- unlist(lapply(split(asgn,  asgn), length))                df <- unlist(lapply(split(asgn,  asgn), length))
714                dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))                #dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
715                ms <- ss/df                ms <- ss/df
716                f <- ms/(ssr/dfr)                #f <- ms/(ssr/dfr)
717                P <- pf(f, df, dfr, lower.tail = FALSE)                #P <- pf(f, df, dfr, lower.tail = FALSE)
718                table <- data.frame(df, ss, ms, dfr, f, P)                #table <- data.frame(df, ss, ms, dfr, f, P)
719                  table <- data.frame(df, ss, ms)
720                dimnames(table) <-                dimnames(table) <-
721                    list(nmeffects,                    list(nmeffects,
722                         c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))  #                       c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
723                           c("Df", "Sum Sq", "Mean Sq"))
724                if ("(Intercept)" %in% nmeffects) table <- table[-1,]                if ("(Intercept)" %in% nmeffects) table <- table[-1,]
725                attr(table, "heading") <- "Analysis of Variance Table"                attr(table, "heading") <- "Analysis of Variance Table"
726                class(table) <- c("anova", "data.frame")                class(table) <- c("anova", "data.frame")
# Line 884  Line 728
728            }            }
729        })        })
730
731  setMethod("update", signature(object = "lmer"),  setMethod("confint", signature(object = "mer"),
function(object, formula., ..., evaluate = TRUE)
{
call <- object@call
if (is.null(call))
stop("need an object with call component")
extras <- match.call(expand.dots = FALSE)$... if (!missing(formula.)) call$formula <- update.formula(formula(object), formula.)
if (length(extras) > 0) {
existing <- !is.na(match(names(extras), names(call)))
for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
if (any(!existing)) {
call <- c(as.list(call), extras[!existing])
call <- as.call(call)
}
}
if (evaluate)
eval(call, parent.frame())
else call
})

setMethod("confint", signature(object = "lmer"),
732            function (object, parm, level = 0.95, ...)            function (object, parm, level = 0.95, ...)
733        {            .NotYetImplemented()
cf <- fixef(object)
pnames <- names(cf)
if (missing(parm))
parm <- seq(along = pnames)
else if (is.character(parm))
parm <- match(parm, pnames, nomatch = 0)
a <- (1 - level)/2
a <- c(a, 1 - a)
pct <- paste(round(100 * a, 1), "%")
ci <- array(NA, dim = c(length(parm), 2),
dimnames = list(pnames[parm], pct))
ses <- sqrt(diag(vcov(object)))[parm]
ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
ci
})

setMethod("param", signature(object = "lmer"),
function(object, unconst = FALSE, ...) {
.Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
})

setMethod("deviance", "lmer",
function(object, REML = NULL, ...) {
.Call("lmer_factor", object, PACKAGE = "Matrix")
if (is.null(REML))
REML <- if (length(oR <- object@REML)) oR else FALSE
object@deviance[[ifelse(REML, "REML", "ML")]]
})

setMethod("deviance", "glmer",
function(object, ...) {
-2 * object@glmmll
})

setMethod("chol", signature(x = "lmer"),
function(x, pivot = FALSE, LINPACK = pivot) {
x@status["factored"] <- FALSE # force a decomposition
.Call("lmer_factor", x, PACKAGE = "Matrix")
})

setMethod("solve", signature(a = "lmer", b = "missing"),
function(a, b, ...)
.Call("lmer_invert", a, PACKAGE = "Matrix")
734            )            )
735
736  setMethod("formula", "lmer", function(x, ...) x@call$formula) setMethod("fitted", signature(object = "mer"), 737 function(object, ...) 738 setMethod("vcov", signature(object = "lmer"), .Call("mer_fitted", object, PACKAGE = "Matrix") 739 function(object, REML = object@REML, useScale = TRUE,...) { ) sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix") rr <- object@RXX nms <- object@cnames[[".fixed"]] dimnames(rr) <- list(nms, nms) nr <- nrow(rr) rr <- rr[-nr, -nr, drop = FALSE] rr <- rr %*% t(rr) if (useScale) { rr = sc^2 * rr } rr }) ## Extract the L matrix setAs("lmer", "dtTMatrix", function(from) { ## force a refactorization if the factors have been inverted if (from@status["inverted"]) from@status["factored"] <- FALSE .Call("lmer_factor", from, PACKAGE = "Matrix") L <- lapply(from@L, as, "dgTMatrix") nf <- length(from@D) Gp <- from@Gp nL <- Gp[nf + 1] Li <- integer(0) Lj <- integer(0) Lx <- double(0) for (i in 1:nf) { for (j in 1:i) { Lij <- L[[Lind(i, j)]] Li <- c(Li, Lij@i + Gp[i]) Lj <- c(Lj, Lij@j + Gp[j]) Lx <- c(Lx, Lij@x) } } new("dtTMatrix", Dim = as.integer(c(nL, nL)), i = Li, j = Lj, x = Lx, uplo = "L", diag = "U") }) 740 741 ## Extract the ZZX matrix setMethod("formula", signature(x = "mer"), 742 setAs("lmer", "dsTMatrix", function(x, ...) 743 function(from) x@call$formula
744    {            )
.Call("lmer_inflate", from, PACKAGE = "Matrix")
ZZpO <- lapply(from@ZZpO, as, "dgTMatrix")
ZZ <- lapply(from@ZtZ, as, "dgTMatrix")
nf <- length(ZZpO)
Gp <- from@Gp
nZ <- Gp[nf + 1]
Zi <- integer(0)
Zj <- integer(0)
Zx <- double(0)
for (i in 1:nf) {
ZZpOi <- ZZpO[[i]]
Zi <- c(Zi, ZZpOi@i + Gp[i])
Zj <- c(Zj, ZZpOi@j + Gp[i])
Zx <- c(Zx, ZZpOi@x)
if (i > 1) {
for (j in 1:(i-1)) {
ZZij <- ZZ[[Lind(i, j)]]
## off-diagonal blocks are transposed
Zi <- c(Zi, ZZij@j + Gp[j])
Zj <- c(Zj, ZZij@i + Gp[i])
Zx <- c(Zx, ZZij@x)
}
}
}
new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx,
uplo = "U")
})
745
746  setMethod("fitted", signature(object = "lmer"),  setMethod("residuals", signature(object = "mer"),
747            function(object, ...)            function(object, ...)
748            napredict(attr(object@frame, "na.action"), object@fitted))            .NotYetImplemented()
749              )
750
751  setMethod("residuals", signature(object = "lmer"),  setMethod("resid", signature(object = "mer"),
752            function(object, ...)            function(object, ...)
753            naresid(attr(object@frame, "na.action"), object@residuals))            .NotYetImplemented()
754              )
755
756  setMethod("resid", signature(object = "lmer"),  setMethod("summary", signature(object = "mer"),
757            function(object, ...) do.call("residuals", c(list(object), list(...))))            function(object, ...)
758              .NotYetImplemented()
759              )
760
761  setMethod("coef", signature(object = "lmer"),  setMethod("update", signature(object = "mer"),
762            function(object, ...)            function(object, ...)
763        {            .NotYetImplemented()
764            fef <- data.frame(rbind(fixef(object)), check.names = FALSE)            )
ref <- as(ranef(object), "list")
names(ref) <- names(object@flist)
val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
for (i in seq(a = val)) {
refi <- ref[[i]]
row.names(val[[i]]) <- row.names(refi)
if (!all(names(refi) %in% names(fef)))
stop("unable to align random and fixed effects")
val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
}
new("lmer.coef", val)
})
765
setMethod("plot", signature(x = "lmer.coef"),
function(x, y, ...)
{
varying <- unique(do.call("c",
lapply(x, function(el)
names(el)[sapply(el,
function(col)
any(col != col[1]))])))
gf <- do.call("rbind", lapply(x, "[", j = varying))
gf$.grp <- factor(rep(names(x), sapply(x, nrow))) switch(min(length(varying), 3), qqmath(eval(substitute(~ x | .grp, list(x = as.name(varying[1])))), gf, ...), xyplot(eval(substitute(y ~ x | .grp, list(y = as.name(varying[1]), x = as.name(varying[2])))), gf, ...), splom(~ gf | .grp, ...)) }) 766 767 setMethod("plot", signature(x = "lmer.ranef"), simss <- function(fm0, fma, nsim) function(x, y, ...) 768 { { 769 lapply(x, function(x) { ysim <- simulate(fm0, nsim) 770 cn <- lapply(colnames(x), as.name) cv <- list(analyticGradient = FALSE, msMaxIter = 200:200, 771 switch(min(ncol(x), 3), msVerbose = 0:0) 772 qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...), sapply(ysim, function(yy) { 773 xyplot(eval(substitute(y ~ x, list(y = cn[[1]], x = cn[[2]]))), .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix") 774 x, ...), LMEoptimize(fm0) <- cv 775 splom(~ x, ...)) .Call("mer_update_y", fma, yy, PACKAGE = "Matrix") 776 LMEoptimize(fma) <- cv 777 exp(c(H0 = fm0@devComp[["logryy2"]], 778 Ha = fma@devComp[["logryy2"]])) 779 }) }) 780 }) } setMethod("with", signature(data = "lmer"), function(data, expr, ...) { dat <- eval(data@call$data)
if (!is.null(na.act <- attr(data@frame, "na.action")))
dat <- dat[-na.act, ]
lst <- c(list(. = data), data@flist, data@frame, dat)
eval(substitute(expr), lst[unique(names(lst))])
})

setMethod("terms", signature(x = "lmer"),
function(x, ...) x@terms)

Legend:
 Removed from v.754 changed lines Added in v.1165