# SCM Repository

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

# Diff of /pkg/R/lmer.R

revision 979, Tue Oct 11 00:02:45 2005 UTC revision 1165, Sat Jan 14 23:41:55 2006 UTC
# Line 2  Line 2
2
3  ## Some utilities  ## Some utilities
4
## Return the index into the packed lower triangle
Lind <- function(i,j) {
if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))
((i - 1) * i)/2 + j
}

5  ## Return the pairs of expressions separated by vertical bars  ## Return the pairs of expressions separated by vertical bars
6  findbars <- function(term)  findbars <- function(term)
7  {  {
# Line 23  Line 17
17  ## that are separated by vertical bars  ## that are separated by vertical bars
18  nobars <- function(term)  nobars <- function(term)
19  {  {
20      # FIXME: is the is.name in the condition redundant?      if (!('|' %in% all.names(term))) return(term)
#   A name won't satisfy the first condition.
if (!('|' %in% all.names(term)) || is.name(term)) return(term)
21      if (is.call(term) && term[[1]] == as.name('|')) return(NULL)      if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
22      if (length(term) == 2) {      if (length(term) == 2) {
23          nb <- nobars(term[[2]])          nb <- nobars(term[[2]])
# Line 60  Line 52
52
53  ## Expand an expression with colons to the sum of the lhs  ## Expand an expression with colons to the sum of the lhs
54  ## and the current expression  ## and the current expression

55  colExpand <- function(term)  colExpand <- function(term)
56  {  {
57      if (is.name(term) || is.numeric(term)) return(term)      if (is.name(term) || is.numeric(term)) return(term)
# Line 116  Line 107
107           analyticHessian = as.logical(analyticHessian))           analyticHessian = as.logical(analyticHessian))
108  }  }
109
110    setMethod("coef", signature(object = "lmer"),
111              function(object, ...)
112          {
113              fef <- data.frame(rbind(object@fixed), check.names = FALSE)
114              ref <- as(ranef(object), "list")
115              names(ref) <- names(object@flist)
116              val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
117              for (i in seq(a = val)) {
118                  refi <- ref[[i]]
119                  row.names(val[[i]]) <- row.names(refi)
120                  if (!all(names(refi) %in% names(fef)))
121                      stop("unable to align random and fixed effects")
122                  val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
123              }
124              new("lmer.coef", val)
125          })
126
127    ## 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"),  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(), start,                     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 , ...)
## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
182        {        {
183            ## match and check parameters            ## match and check parameters
184            if (length(formula) < 3) stop("formula must be a two-sided formula")            if (length(formula) < 3) stop("formula must be a two-sided formula")
185            cv <- do.call("lmerControl", control)            cv <- do.call("lmerControl", control)
186            ## evaluate glm.fit, a generalized linear fit of fixed effects only
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()            mf <- match.call()
191            m <- match(c("family", "data", "subset", "weights",            m <- match(c("family", "data", "subset", "weights",
192                         "na.action", "offset"), names(mf), 0)                         "na.action", "offset"), names(mf), 0)
193            mf <- mf[c(1, m)]            mf <- fe <- mf[c(1, m)]
194            frame.form <- subbars(formula) # substitute `+' for `|'            frame.form <- subbars(formula) # substitute `+' for `|'
195            fixed.form <- nobars(formula)  # remove any terms with `|'            fixed.form <- nobars(formula)  # remove any terms with `|'
196            if (inherits(fixed.form, "name")) # RHS is empty - use a constant            if (inherits(fixed.form, "name")) # RHS is empty - use a constant
197                fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))                fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
198            environment(fixed.form) <- environment(frame.form) <- environment(formula)            environment(fixed.form) <- environment(frame.form) <- environment(formula)
199            mf\$formula <- fixed.form
200            mf\$x <- mf\$model <- mf\$y <- TRUE            ## evaluate a model frame for fixed and random effects
201            mf[[1]] <- as.name("glm")            mf\$formula <- frame.form
202            glm.fit <- eval(mf, parent.frame())            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            x <- glm.fit\$x
215            y <- as.double(glm.fit\$y)            y <- as.double(glm.fit\$y)
216            family <- glm.fit\$family            family <- glm.fit\$family
217
218            ## check for a linear mixed model            ## check for a linear mixed model
219            lmm <- family\$family == "gaussian" && family\$link == "identity"            lmm <- family\$family == "gaussian" && family\$link == "identity"
220            if (lmm) { # linear mixed model            if (lmm) { # linear mixed model
# Line 165  Line 236
236                                '\nUsing method = "PQL".\n')                                '\nUsing method = "PQL".\n')
237                }                }
238            }            }
239              ## create factor list for the random effects
## evaluate a model frame for fixed and random effects
mf\$formula <- frame.form
mf\$x <- mf\$model <- mf\$y <- mf\$family <- NULL
mf\$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
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, function(x)            fl <- lapply(bars,
243                       list(model.matrix(eval(substitute(~ T, list(T = x[[2]]))),                         function(x)
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                  fl <- fl[ord]
252            ## Create the model matrices and a mixed-effects representation (mer)            }
253            mmats <- c(lapply(random, "[[", 1),            ## create list of transposed model matrices for random effects
254                       .fixed = list(cbind(glm.fit\$x, .response = glm.fit\$y)))            Ztl <- lapply(bars, function(x)
255            mer <- .Call("lmer_create", lapply(random, "[[", 2),                          t(model.matrix(eval(substitute(~ expr,
256                         mmats, method, PACKAGE = "Matrix")                                                         list(expr = x[[2]]))),
257            if (lmm) {                    ## linear mixed model                                         frm)))
258                if (missing(start)) .Call("lmer_initial", mer, PACKAGE="Matrix")            ## Create the mixed-effects representation (mer) object
259                else .Call("lmer_set_initial", mer, start, PACKAGE = "Matrix")            mer <- .Call("mer_create", fl,
260                .Call("lmer_ECMEsteps", mer, cv\$niterEM, cv\$EMverbose, PACKAGE = "Matrix")                         .Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"),
261                           x, y, method, sapply(Ztl, nrow),
262                           c(lapply(Ztl, rownames), list(.fixed = colnames(x))),
263                           !(family\$family %in% c("binomial", "poisson")),
264                           match.call(), family,
265                           PACKAGE = "Matrix")
266              if (lmm) {
267                  .Call("mer_ECMEsteps", mer, cv\$niterEM, cv\$EMverbose, PACKAGE = "Matrix")
268                LMEoptimize(mer) <- cv                LMEoptimize(mer) <- cv
269                fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")                return(mer)
return(new("lmer",
mer,
assign = attr(x, "assign"),
call = match.call(),
family = family, fitted = fits,
fixed = fixef(mer),
frame = if (model) frm else data.frame(),
logLik = logLik(mer),
residuals = unname(model.response(frm) - fits),
terms = glm.fit\$terms))
270            }            }
271
272            ## The rest of the function applies to generalized linear mixed models            ## The rest of the function applies to generalized linear mixed models
# Line 219  Line 276
276            wtssqr <- wts * wts            wtssqr <- wts * wts
277            offset <- glm.fit\$offset            offset <- glm.fit\$offset
278            if (is.null(offset)) offset <- numeric(length(eta))            if (is.null(offset)) offset <- numeric(length(eta))
mu <- numeric(length(eta))

dev.resids <- quote(family\$dev.resids(y, mu, wtssqr))
279            linkinv <- quote(family\$linkinv(eta))            linkinv <- quote(family\$linkinv(eta))
280            mu.eta <- quote(family\$mu.eta(eta))            mu.eta <- quote(family\$mu.eta(eta))
281              mu <- family\$linkinv(eta)
282            variance <- quote(family\$variance(mu))            variance <- quote(family\$variance(mu))
283              dev.resids <- quote(family\$dev.resids(y, mu, wtssqr))
284            LMEopt <- get("LMEoptimize<-")            LMEopt <- get("LMEoptimize<-")
285            doLMEopt <- quote(LMEopt(x = mer, value = cv))            doLMEopt <- quote(LMEopt(x = mer, value = cv))
286
287            GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")            GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
288            .Call("glmer_PQL", GSpt, PACKAGE = "Matrix")  # obtain PQL estimates            .Call("glmer_PQL", GSpt, PACKAGE = "Matrix")  # obtain PQL estimates

289            fixInd <- seq(ncol(x))            fixInd <- seq(ncol(x))
290            ## pars[fixInd] == beta, pars[-fixInd] == theta            ## pars[fixInd] == beta, pars[-fixInd] == theta
291            PQLpars <- c(fixef(mer),            PQLpars <- c(fixef(mer),
292                         .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))                         .Call("mer_coef", mer, 2, PACKAGE = "Matrix"))
293            ## set flag to skip fixed-effects in subsequent calls            .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")
mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
294            ## indicator of constrained parameters            ## indicator of constrained parameters
295            const <- c(rep(FALSE, length(fixInd)),            const <- c(rep(FALSE, length(fixInd)),
296                       unlist(lapply(mer@nc[seq(along = random)],                       unlist(lapply(mer@nc[seq(along = fl)],
297                                     function(k) 1:((k*(k+1))/2) <= k)                                     function(k) 1:((k*(k+1))/2) <= k)
298                              ))                              ))
299            devAGQ <- function(pars, n)            devLaplace <- function(pars)
300                .Call("glmer_devAGQ", pars, GSpt, n, PACKAGE = "Matrix")                .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)            deviance <- devAGQ(PQLpars, 1)
310
311  ### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs  ### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs
312  ### AGQ for nc > 1 first.  ### AGQ for nc > 1 first.
313            fxd <- PQLpars[fixInd]            fxd <- PQLpars[fixInd]
# Line 265  Line 327
327                }                }
328                obj <- function(pars)                obj <- function(pars)
329                    .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")                    .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
if (exists("nlminb", mode = "function")) {
330                    optimRes <-                    optimRes <-
331                        nlminb(PQLpars, obj,                        nlminb(PQLpars, obj,
332                               lower = ifelse(const, 5e-10, -Inf),                               lower = ifelse(const, 5e-10, -Inf),
# Line 275  Line 336
336                    if (optimRes\$convergence != 0)                    if (optimRes\$convergence != 0)
337                        warning("nlminb failed to converge")                        warning("nlminb failed to converge")
338                    deviance <- optimRes\$objective                    deviance <- optimRes\$objective
339                } else {                if (gVerb)
optimRes <-
optim(PQLpars, obj, method = "L-BFGS-B",
lower = ifelse(const, 5e-10, -Inf),
control = list(trace = getOption("verbose"),
maxit = cv\$msMaxIter))
optpars <- optimRes\$par
if (optimRes\$convergence != 0)
warning("optim failed to converge")
deviance <- optimRes\$value
}
if (gVerb) {
340                    cat(paste("convergence message", optimRes\$message, "\n"))                    cat(paste("convergence message", optimRes\$message, "\n"))
}
341                fxd[] <- optpars[fixInd]  ## preserve the names                fxd[] <- optpars[fixInd]  ## preserve the names
342                .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")                .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
343            }            }
# Line 302  Line 351
351                call = match.call(), family = family,                call = match.call(), family = family,
352                logLik = loglik, fixed = fxd)                logLik = loglik, fixed = fxd)
353        })        })
## end{  "lmer . formula " }
354
355
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"),  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
381                   constr <- unlist(lapply(nc[1:(length(nc) - 2)],                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
function(k) 1:((k*(k+1))/2) <= k))
382                   fn <- function(pars)                   fn <- function(pars)
383                       deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))                       deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
384                   gr <-                   gr <- if (value\$analyticGradient)
if (value\$analyticGradient)
385                           function(pars) {                           function(pars) {
386                               if (!isTRUE(all.equal(pars,                               if (!isTRUE(all.equal(pars,
387                                                     .Call("lmer_coef", x,                                                 .Call("mer_coef", x,
388                                                           2, PACKAGE = "Matrix"))))                                                           2, PACKAGE = "Matrix"))))
389                                   .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")                               .Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")
390                               .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")                           .Call("mer_gradient", x, 2, PACKAGE = "Matrix")
391                           }                           }
392                   ## else NULL                   else NULL
393                   optimRes <-                   optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"),
if (exists("nlminb", mode = "function"))
nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
394                                  fn, gr,                                  fn, gr,
395                                  lower = ifelse(constr, 5e-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")
optim(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
fn, gr, method = "L-BFGS-B",
lower = ifelse(constr, 5e-10, -Inf),
control = list(maxit = value\$msMaxIter,
trace = as.integer(value\$msVerbose)))
.Call("lmer_coefGets", x, optimRes\$par, 2, PACKAGE = "Matrix")
399                   if (optimRes\$convergence != 0) {                   if (optimRes\$convergence != 0) {
400                       warning(paste("optim or nlminb returned message",                       warning(paste("nlminb returned message",
401                                     optimRes\$message,"\n"))                                     optimRes\$message,"\n"))
402                   }                   }
403                   return(x)                   return(x)
404               })               })
405
406  setMethod("ranef", signature(object = "lmer"),  setMethod("deviance", signature(object = "mer"),
407            function(object, accumulate = FALSE, ...) {            function(object, ...) {
408                val <- new("lmer.ranef",                .Call("mer_factor", object, PACKAGE = "Matrix")
409                           lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),                object@deviance[[ifelse(object@method == "REML", "REML", "ML")]]
data.frame, check.names = FALSE),
varFac = object@bVar,
stdErr = .Call("lmer_sigma", object,
object@method == "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
410            })            })
411
412  setMethod("fixef", signature(object = "mer"),  setMethod("mcmcsamp", signature(object = "mer"),
413            function(object, ...)            function(object, n = 1, verbose = FALSE, saveb = FALSE,
414                .Call("lmer_fixef", object, PACKAGE = "Matrix"))                     trans = TRUE, ...)
415          {
416              ans <- t(.Call("mer_MCMCsamp", object, saveb, n,
417                             trans, PACKAGE = "Matrix"))
418              attr(ans, "mcpar") <- as.integer(c(1, n, 1))
419              class(ans) <- "mcmc"
420              glmer <- FALSE
421              gnms <- names(object@flist)
422              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("fixef", signature(object = "lmer"),  setMethod("simulate", signature(object = "mer"),
448            function(object, ...) object@fixed)            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  setMethod("VarCorr", signature(x = "lmer"),            family <- object@family
462  ##FIXME - change this for reasonable defaults of useScale according to            if (family\$family != "gaussian" ||
463            ##the family slot.                family\$link != "identity")
464            function(x, REML = TRUE, useScale = TRUE, ...) {                stop("simulation of generalized linear mixed models not yet implemented")
465                val <- .Call("lmer_variances", x, PACKAGE = "Matrix")            ## similate the linear predictors
466                for (i in seq(along = val)) {            lpred <- .Call("mer_simulate", object, nsim, PACKAGE = "Matrix")
467                    dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])            sc <- 1
468                    val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")            if (object@useScale)
469                }                sc <- .Call("mer_sigma", object, object@method == "REML",
470                new("VarCorr",                            PACKAGE = "Matrix")
471                    scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),            ## add fixed-effects contribution and per-observation noise term
472                    reSumry = val,            lpred <- as.data.frame(lpred + drop(object@X %*% fixef(object)) +
473                    useScale = useScale)                                   rnorm(prod(dim(lpred)), sd = sc))
474              ## save the seed
475              attr(lpred, "seed") <- RNGstate
476              lpred
477            })            })
478
setMethod("gradient", signature(x = "lmer"),
function(x, unconst, ...)
.Call("lmer_gradient", x, unconst, PACKAGE = "Matrix"))

setMethod("summary", signature(object = "lmer"),
function(object, ...)
new("summary.lmer", object,
showCorrelation = TRUE,
useScale = !((object@family)\$family %in% c("binomial", "poisson"))))

setMethod("show", signature(object = "lmer"),
function(object)
show(new("summary.lmer", object,
showCorrelation = FALSE,
useScale = !((object@family)\$family %in% c("binomial", "poisson")))))
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 <- object@method == "REML"                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 447  Line 559
559                                 row.names = ""))                                 row.names = ""))
560                }                }
561                cat("Random effects:\n")                cat("Random effects:\n")
562                show(VarCorr(object, useScale = useScale))                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, FALSE, 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 473  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              function(object, REML = object@method == "REML",
608                       useScale = object@useScale,...) {
609                  sc <- if (object@useScale) {
610                      .Call("mer_sigma", object, REML, PACKAGE = "Matrix")
611                  } else { 1 }
612                  rr <- as(sc^2 * tcrossprod(solve(object@RXX)), "dpoMatrix")
613                  rr@factors\$correlation <- as(rr, "correlation")
614                  rr
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 <- abs(nc[1]) - 1
n <- nc[2]
rep(n-p, p)
626        })        })
627
628  setMethod("logLik", signature(object="mer"),  setMethod("logLik", signature(object="mer"),
629            function(object, REML = object@method == "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") <- abs(nc[1]) +                attr(val, "df") <- abs(devc[2]) +
634                    length(.Call("lmer_coef", object, 0, PACKAGE = "Matrix"))                    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="lmer"),  setMethod("VarCorr", signature(x = "mer"),
641            function(object, ...) object@logLik)            function(x, REML = x@method == "REML", useScale = x@useScale, ...)
642          {
643              sc <- 1
644              if (useScale)
645                  sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix")
646              sc2 <- sc * sc
647              ans <- lapply(x@Omega, function(el) {
648                  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 566  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 580  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 595  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("deviance", "mer",
function(object, REML = NULL, ...) {
.Call("lmer_factor", object, PACKAGE = "Matrix")
if (is.null(REML))
REML <- object@method == "REML"
object@deviance[[ifelse(REML, "REML", "ML")]]
})

setMethod("deviance", "lmer",
function(object, ...) -2 * c(object@logLik))

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@method == "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              )
setMethod("resid", signature(object = "lmer"),
function(object, ...) do.call("residuals", c(list(object), list(...))))
755
756  setMethod("coef", signature(object = "lmer"),  setMethod("summary", signature(object = "mer"),
757            function(object, ...)            function(object, ...)
758        {            .NotYetImplemented()
759            fef <- data.frame(rbind(object@fixed), 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)
})

setMethod("plot", signature(x = "lmer.coef"),
function(x, y, ...)
{
## require("lattice", quietly = TRUE) -- now via Imports
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, ...))
})

setMethod("plot", signature(x = "lmer.ranef"),
function(x, y, ...)
{
## require("lattice", quietly = TRUE) -- now via Imports
lapply(x, function(x) {
cn <- lapply(colnames(x), as.name)
switch(min(ncol(x), 3),
qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
xyplot(eval(substitute(y ~ x,
list(y = cn[[1]],
x = cn[[2]]))), x, ...),
splom(~ x, ...))
})
})

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)

setMethod("show", signature(object="VarCorr"),
function(object)
{
digits <- max(3, getOption("digits") - 2)
useScale <- length(object@useScale) > 0 && object@useScale[1]
sc <- ifelse(useScale, object@scale,  1.)
reStdDev <- c(lapply(object@reSumry,
function(x, sc)
sc*x@stdDev,
sc = sc), list(Residual = sc))
reLens <- unlist(c(lapply(reStdDev, length)))
reMat <- array('', c(sum(reLens), 4),
list(rep('', sum(reLens)),
c("Groups", "Name", "Variance", "Std.Dev.")))
reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
reMat[,2] <- c(unlist(lapply(reStdDev, names)), "")
reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
reMat[,4] <- format(unlist(reStdDev), digits = digits)
if (any(reLens > 1)) {
maxlen <- max(reLens)
corr <-
do.call("rbind",
lapply(object@reSumry,
function(x, maxlen) {
cc <- format(round(x, 3), nsmall = 3)
cc[!lower.tri(cc)] <- ""
nr <- dim(cc)[1]
if (nr >= maxlen) return(cc)
cbind(cc, matrix("", nr, maxlen-nr))
}, maxlen))
colnames(corr) <- c("Corr", rep("", maxlen - 1))
reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr))))
}
if (!useScale) reMat <- reMat[-nrow(reMat),]
print(reMat, quote = FALSE)
})

setMethod("mcmcsamp", signature(object = "lmer"),
function(object, n = 1, verbose = FALSE, saveb = FALSE,
trans = TRUE, ...)
{
if (object@family\$family == "gaussian" &&
object@family\$link == "identity") {
glmer <- FALSE
ans <- .Call("lmer_MCMCsamp", object, saveb, n, trans,
PACKAGE = "Matrix")
} else {
glmer <- TRUE
if (trans)
warning("trans option not currently allowed for generalized models")
trans <- FALSE
## Check arguments
if (length(object@Omega) > 1 || object@nc[1] > 1)
stop("mcmcsamp currently defined for glmm models with only one variance component")
cv <- Matrix:::lmerControl()
if (verbose) cv\$msVerbose <- 1
family <- object@family
frm <- object@frame

## recreate model matrices
fixed.form <- Matrix:::nobars(object@call\$formula)
if (inherits(fixed.form, "name")) # RHS is empty - use a constant
fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE,
y = TRUE)
x <- glm.fit\$x
y <- as.double(glm.fit\$y)
bars <- Matrix:::findbars(object@call\$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]])))
if (any(names(random) != names(object@flist)))
random <- random[names(object@flist)]
mmats <- c(lapply(random, "[[", 1),
.fixed = list(cbind(glm.fit\$x, .response = glm.fit\$y)))
mer <- as(object, "mer")

## establish the GS object and the ans matrix
eta <- glm.fit\$linear.predictors # perhaps later change this to object@fitted?
wts <- glm.fit\$prior.weights
wtssqr <- wts * wts
offset <- glm.fit\$offset
if (is.null(offset)) offset <- numeric(length(eta))
off <- numeric(length(eta))
mu <- numeric(length(eta))
dev.resids <- quote(family\$dev.resids(y, mu, wtssqr))
linkinv <- quote(family\$linkinv(eta))
mu.eta <- quote(family\$mu.eta(eta))
variance <- quote(family\$variance(mu))
LMEopt <- getAnywhere("LMEoptimize<-")
doLMEopt <- quote(LMEopt(x = mer, value = cv))
GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
fixed <- object@fixed
varc <- .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")
b <- .Call("lmer_ranef", mer, PACKAGE = "Matrix")
ans <- .Call("glmer_MCMCsamp", GSpt, b, fixed, varc, saveb, n,
PACKAGE = "Matrix")
.Call("glmer_finalize", GSpt, PACKAGE = "Matrix");
}
gnms <- names(object@flist)
cnms <- object@cnames
ff <- fixef(object)
colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
unlist(lapply(seq(along = gnms),
function(i)
abbrvNms(gnms[i],cnms[[i]]))))
if (trans) {
## parameter type: 0 => fixed effect, 1 => variance,
##                 2 => covariance
ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1,
unlist(lapply(seq(along = gnms),
function(i)
{
k <- length(cnms[[i]])
rep(1:2, c(k, (k*(k-1))/2))
})))
colnms[ptyp == 1] <-
paste("log(", colnms[ptyp == 1], ")", sep = "")
colnms[ptyp == 2] <-
paste("atanh(", colnms[ptyp == 2], ")", sep = "")
}
colnames(ans) <- colnms
ans
})

rWishart <- function(n, df, invScal)
.Call("Matrix_rWishart", n, df, invScal)

760
761  setMethod("model.matrix", signature(object = "lmer"),  setMethod("update", signature(object = "mer"),
762            function(object, ...)            function(object, ...)
763        {            .NotYetImplemented()
764            frm <- object@frame            )
fixed.form <- Matrix:::nobars(object@call\$formula)
if (inherits(fixed.form, "name")) # RHS is empty - use a constant
fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
glm.fit <- glm(eval(fixed.form), object@family, frm, x = TRUE, y = TRUE)
fxd <- unname(drop(glm.fit\$x %*% fixef(object)))

## Create the random effects model matrices
bars <- Matrix:::findbars(object@call\$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]])))
## re-order the random effects pairs if necessary
if (any(names(random) != names(object@flist)))
random <- random[names(object@flist)]
c(lapply(random, "[[", 1),
.fixed = list(cbind(glm.fit\$x, .response = glm.fit\$y)))
})

setMethod("simulate", signature(object = "lmer"),
function(object, nsim = 1, seed = NULL, ...)
{
if(!exists(".Random.seed", envir = .GlobalEnv))
runif(1)               # initialize the RNG if necessary
if(is.null(seed))
RNGstate <- .Random.seed
else {
R.seed <- .Random.seed
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}

family <- object@family
if (family\$family != "gaussian" ||
family\$link != "identity")
stop("simulation of generalized linear mixed models not yet implemented")

## pieces we will need later
scale <- .Call("lmer_sigma", object, object@method == "REML",
PACKAGE = "Matrix")
mmats <- model.matrix(object)
ff <- fixef(object)

###_FIXME: If the factor levels have been permuted, has the
###        permutation been applied in the stored frame?  Otherwise we
###        need to check this.

## similate the linear predictors
lpred <- .Call("lmer_simulate", as(object, "mer"), nsim,
unname(drop(mmats[[length(mmats)]][,seq(a = ff)] %*% ff)),
mmats, TRUE, PACKAGE = "Matrix")
## add per-observation noise term
lpred <- as.data.frame(lpred + rnorm(prod(dim(lpred)), sd = scale))
765
## save the seed
attr(lpred, "seed") <- RNGstate
lpred
})
766
767  simulate2 <- function(object, n = 1, ...)  simss <- function(fm0, fma, nsim)
768  {  {
769      family <- object@family      ysim <- simulate(fm0, nsim)
770      if (family\$family != "gaussian" ||      cv <- list(analyticGradient = FALSE, msMaxIter = 200:200,
771          family\$link != "identity")                 msVerbose = 0:0)
772          stop("simulation of generalized linear mixed models not implemented yet")      sapply(ysim, function(yy) {
773            .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")
774      ## create the mean from the fixed effects          LMEoptimize(fm0) <- cv
775      frm <- object@frame          .Call("mer_update_y", fma, yy, PACKAGE = "Matrix")
776      fixed.form <- Matrix:::nobars(object@call\$formula)          LMEoptimize(fma) <- cv
777      if (inherits(fixed.form, "name")) # RHS is empty - use a constant          exp(c(H0 = fm0@devComp[["logryy2"]],
778          fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))                Ha = fma@devComp[["logryy2"]]))
glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE, y = TRUE)
lpred <- matrix(glm.fit\$x %*% fixef(object), nr = nrow(frm), nc = n)

## Create the random effects model matrices
bars <- Matrix:::findbars(object@call\$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]])))
## re-order the random effects pairs if necessary
flist <- object@flist
if (any(names(random) != names(flist)))
random <- random[names(flist)]
mmats <- lapply(random, "[[", 1)

## simulate the random effects
scale <- .Call("lmer_sigma", object, object@method == "REML",
PACKAGE = "Matrix")
Omega <- object@Omega
re <- lapply(seq(along = Omega),
function(i) {
om <- Omega[[i]]
nr <- nrow(om)
nlev <- length(levels(flist[[i]]))
scale * array(solve(chol(new("dpoMatrix", Dim = dim(om),
uplo = "U", x = c(om))),
matrix(rnorm(nr * n * nlev),
nr = nr))@x, c(nr, n, nlev))
779                   })                   })
## apply the random effects
for (j in seq(along = Omega)) {
for (i in 1:nrow(lpred))
lpred[i,] <- lpred[i,] + mmats[[j]][i,] %*% re[[j]][, , as.integer(flist[[j]])[i]]
}
## add per-observation noise term
lpred <- lpred + rnorm(prod(dim(lpred)), sd = scale)
attr(lpred, "re") <- re
lpred
}

refdist <- function(fm1, fm2, n, ...)
{
cv <- lmerControl()
obs <- deviance(fm2) - deviance(fm1)
newy <- simulate(fm2, n)
mm1 <- model.matrix(fm1)
mm2 <- model.matrix(fm2)
ref <- numeric(n)
mer1 <- as(fm1, "mer")
mer2 <- as(fm2, "mer")
for (j in 1:n) {
.Call("lmer_update_y", mer2, newy[[j]], mm2, PACKAGE = "Matrix")
LMEoptimize(mer2) <- cv
.Call("lmer_update_y", mer1, newy[[j]], mm1, PACKAGE = "Matrix")
LMEoptimize(mer1) <- cv
ref[j] <- deviance(mer2) - deviance(mer1)
}
attr(ref, "observed") <- obs
ref
780  }  }

mer2 <-
function(formula, data, family,
method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
control = list(), start,
subset, weights, na.action, offset,
model = TRUE, x = FALSE, y = FALSE , ...)
## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
{
## match and check parameters
if (length(formula) < 3) stop("formula must be a two-sided formula")
## cv <- do.call("Matrix:::lmerControl", control)
cv <- control
cv\$analyticGradient <- FALSE
cv\$msMaxIter <- as.integer(200)
if (is.null(cv\$msVerbose)) cv\$msVerbose <- as.integer(1)
## FIXME: Need to evaluate the model frame first and then fit the glm
##  in that frame.  Otherwise missing values in the grouping factors
##  cause problems.
## evaluate glm.fit, a generalized linear fit of fixed effects only
mf <- match.call()
m <- match(c("family", "data", "subset", "weights",
"na.action", "offset"), names(mf), 0)
mf <- mf[c(1, m)]
frame.form <- subbars(formula) # substitute `+' for `|'
fixed.form <- nobars(formula)  # remove any terms with `|'
if (inherits(fixed.form, "name")) # RHS is empty - use a constant
fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
environment(fixed.form) <- environment(frame.form) <- environment(formula)
mf\$formula <- fixed.form
mf\$x <- mf\$model <- mf\$y <- TRUE
mf[[1]] <- as.name("glm")
glm.fit <- eval(mf, parent.frame())
x <- glm.fit\$x
y <- as.double(glm.fit\$y)
family <- glm.fit\$family
## check for a linear mixed model
lmm <- family\$family == "gaussian" && family\$link == "identity"
if (lmm) { # linear mixed model
method <- match.arg(method)
if (method %in% c("PQL", "Laplace", "AGQ")) {
warning(paste('Argument method = "', method,
'" is not meaningful for a linear mixed model.\n',
'Using method = "REML".\n', sep = ''))
method <- "REML"
}
} else { # generalized linear mixed model
if (missing(method)) method <- "PQL"
else {
method <- match.arg(method)
if (method == "ML") method <- "PQL"
if (method == "REML")
warning('Argument method = "REML" is not meaningful ',
'for a generalized linear mixed model.',
'\nUsing method = "PQL".\n')
}
}

## evaluate a model frame for fixed and random effects
mf\$formula <- frame.form
mf\$x <- mf\$model <- mf\$y <- mf\$family <- NULL
mf\$drop.unused.levels <- TRUE
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(t(model.matrix(eval(substitute(~ expr,
list(expr = x[[2]]))),
frm)),
eval(substitute(as.factor(fac)[,drop = TRUE],
list(fac = x[[3]])), frm)))
names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))

## order factor list by decreasing number of levels
nlev <- sapply(random, function(x) length(levels(x[[2]])))
if (any(diff(nlev) > 0)) {
random <- random[rev(order(nlev))]
}

## Create the model matrices and a mixed-effects representation (mer)
mer <- .Call("mer2_create", random, x, y, method, PACKAGE = "Matrix")
LMEoptimize(mer) <- cv
mer
}

## Extract the permutation
setAs("mer2", "pMatrix", function(from) .Call("mer2_pMatrix", from))

## Extract the L matrix
setAs("mer2", "dtCMatrix", function(from) .Call("mer2_dtCMatrix", from))

## Optimization for mer2 objects
setReplaceMethod("LMEoptimize", signature(x="mer2", value="list"),
function(x, value)
{
if (value\$msMaxIter < 1) return(x)
nc <- x@nc
constr <- unlist(lapply(nc[1:(length(nc) - 2)],
function(k) 1:((k*(k+1))/2) <= k))
fn <- function(pars)
deviance(.Call("mer2_coefGets", x, pars, 2, PACKAGE = "Matrix"))
gr <- NULL  ## No gradient yet
##                     if (value\$analyticGradient)
##                         function(pars) {
##                             if (!isTRUE(all.equal(pars,
##                                                   .Call("lmer_coef", x,
##                                                         2, PACKAGE = "Matrix"))))
##                                 .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
##                             .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
##                         }
## else NULL
optimRes <- nlminb(.Call("mer2_coef", x, 2, PACKAGE = "Matrix"),
fn, gr,
lower = ifelse(constr, 5e-10, -Inf),
control = list(iter.max = value\$msMaxIter,
trace = as.integer(value\$msVerbose)))
.Call("mer2_coefGets", x, optimRes\$par, 2, PACKAGE = "Matrix")
if (optimRes\$convergence != 0) {
warning(paste("nlminb returned message",
optimRes\$message,"\n"))
}
return(x)
})

setMethod("deviance", "mer2",
function(object, REML = NULL, ...) {
.Call("mer2_factor", object, PACKAGE = "Matrix")
if (is.null(REML))
REML <- object@method == "REML"
object@deviance[[ifelse(REML, "REML", "ML")]]
})

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

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: