# SCM Repository

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

# Diff of /pkg/R/lmer.R

branches/trunk-lme4/R/lmeRep.R revision 382, Sun Dec 12 17:30:29 2004 UTC pkg/R/lmer.R revision 1165, Sat Jan 14 23:41:55 2006 UTC
# Line 1  Line 1
1  setReplaceMethod("LMEoptimize", signature(x="lmeRep", value="list"),  # Methods for lmer and for the objects that it produces
2
3    ## Some utilities
4
5    ## Return the pairs of expressions separated by vertical bars
6    findbars <- function(term)
7    {
8        if (is.name(term) || is.numeric(term)) return(NULL)
9        if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
10        if (!is.call(term)) stop("term must be of class call")
11        if (term[[1]] == as.name('|')) return(term)
12        if (length(term) == 2) return(findbars(term[[2]]))
13        c(findbars(term[[2]]), findbars(term[[3]]))
14    }
15
16    ## Return the formula omitting the pairs of expressions
17    ## that are separated by vertical bars
18    nobars <- function(term)
19    {
20        if (!('|' %in% all.names(term))) return(term)
21        if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
22        if (length(term) == 2) {
23            nb <- nobars(term[[2]])
24            if (is.null(nb)) return(NULL)
25            term[[2]] <- nb
26            return(term)
27        }
28        nb2 <- nobars(term[[2]])
29        nb3 <- nobars(term[[3]])
30        if (is.null(nb2)) return(nb3)
31        if (is.null(nb3)) return(nb2)
32        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"),
91               niterEM = 15,
92               EMverbose = getOption("verbose"),
93               PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead
95               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              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"),
177              function(formula, data, family,
178                       method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
179                       control = list(), start,
180                       subset, weights, na.action, offset,
181                       model = TRUE, x = FALSE, y = FALSE , ...)
182          {
183              ## 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)
222                  if (method %in% c("PQL", "Laplace", "AGQ")) {
223                      warning(paste('Argument method = "', method,
224                                    '" is not meaningful for a linear mixed model.\n',
225                                    'Using method = "REML".\n', sep = ''))
226                      method <- "REML"
227                  }
228              } else { # generalized linear mixed model
229                  if (missing(method)) method <- "PQL"
230                  else {
231                      method <- match.arg(method)
232                      if (method == "ML") method <- "PQL"
233                      if (method == "REML")
234                          warning('Argument method = "REML" is not meaningful ',
235                                  'for a generalized linear mixed model.',
236                                  '\nUsing method = "PQL".\n')
237                  }
238              }
239              ## create factor list for the random effects
240              bars <- findbars(formula[[3]])
241              names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
242              fl <- lapply(bars,
243                           function(x)
244                           eval(substitute(as.factor(fac)[,drop = TRUE],
245                                           list(fac = x[[3]])), frm))
246              ## order factor list by decreasing number of levels
247              nlev <- sapply(fl, function(x) length(levels(x)))
248              if (any(diff(nlev) > 0)) {
249                  ord <- rev(order(nlev))
250                  bars <- bars[ord]
251                  fl <- fl[ord]
252              }
253              ## create list of transposed model matrices for random effects
254              Ztl <- lapply(bars, function(x)
255                            t(model.matrix(eval(substitute(~ expr,
256                                                           list(expr = x[[2]]))),
257                                           frm)))
258              ## Create the mixed-effects representation (mer) object
259              mer <- .Call("mer_create", fl,
260                           .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
269                  return(mer)
270              }
271
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))
280              mu.eta <- quote(family\$mu.eta(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
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)
st <- ccoef(x)         # starting values
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                       ccoef(x) <- pars                       deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
deviance(x, REML = value\$REML)
}
385                       function(pars) {                       function(pars) {
386                           ccoef(x) <- pars                           if (!isTRUE(all.equal(pars,
387                           grad <- lme4:::gradient(x, REML = value\$REML, unconst = TRUE)                                                 .Call("mer_coef", x,
389                           grad                               .Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")
390                       } else NULL                           .Call("mer_gradient", x, 2, PACKAGE = "Matrix")
391                   optimRes <- optim(st, fn, gr,                       }
392                                     method = "L-BFGS-B",                   else NULL
393                                     lower = ifelse(constr, 1e-10, -Inf),                   optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"),
394                                     control = list(maxit = value\$msMaxIter,                                      fn, gr,
395                                     trace = ifelse(value\$msVerbose,6,0)))                                      lower = ifelse(constr, 5e-10, -Inf),
396                                        control = list(iter.max = value\$msMaxIter,
397                                        trace = as.integer(value\$msVerbose)))
398                     .Call("mer_coefGets", x, optimRes\$par, 2, PACKAGE = "Matrix")
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                   }                   }
ccoef(x) = optimRes\$par
403                   return(x)                   return(x)
404               })               })
405
406  setMethod("deviance", signature(object = "lmeRep"),  setMethod("deviance", signature(object = "mer"),
function(object, REML = FALSE, ...) {
.Call("lmeRep_factor", object, PACKAGE = "Matrix")
object@deviance[ifelse(REML, 2, 1)]
})

setMethod("ranef", signature(object = "lmeRep"),
407            function(object, ...) {            function(object, ...) {
408                .Call("lmeRep_ranef", object, PACKAGE = "Matrix")                .Call("mer_factor", object, PACKAGE = "Matrix")
409                  object@deviance[[ifelse(object@method == "REML", "REML", "ML")]]
410            })            })
411
412  setMethod("fixef", signature(object = "lmeRep"),  setMethod("mcmcsamp", signature(object = "mer"),
413            function(object, ...) {            function(object, n = 1, verbose = FALSE, saveb = FALSE,
414                val = .Call("lmeRep_fixef", object, PACKAGE = "Matrix")                     trans = TRUE, ...)
415                names(val) = object@cnames[[".fixed"]]        {
416                val[-length(val)]            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("vcov", signature(object = "lmeRep"),  setMethod("simulate", signature(object = "mer"),
448            function(object, REML = TRUE, useScale = TRUE,...) {            function(object, nsim = 1, seed = NULL, ...)
449                ## force an "lmeRep_invert"        {
450                sc <- .Call("lmeRep_sigma", object, REML, PACKAGE = "Matrix")            if(!exists(".Random.seed", envir = .GlobalEnv))
451                rr <- object@RXX                runif(1)               # initialize the RNG if necessary
452                nms <- object@cnames[[".fixed"]]            if(is.null(seed))
453                dimnames(rr) <- list(nms, nms)                RNGstate <- .Random.seed
454                nr <- nrow(rr)            else {
455                rr <- rr[-nr, -nr, drop = FALSE]                R.seed <- .Random.seed
456                rr <- rr %*% t(rr)                set.seed(seed)
457                if (useScale) {                RNGstate <- structure(seed, kind = as.list(RNGkind()))
458                    rr = sc^2 * rr                on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
459                }                }
rr
})
460
461  setMethod("VarCorr", signature(x = "lmeRep"),            family <- object@family
462            function(x, REML = TRUE, useScale = TRUE, ...) {            if (family\$family != "gaussian" ||
463                val = .Call("lmeRep_variances", x, PACKAGE = "Matrix")                family\$link != "identity")
464                for (i in seq(along = val)) {                stop("simulation of generalized linear mixed models not yet implemented")
465                    dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])            ## similate the linear predictors
466                    val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")            lpred <- .Call("mer_simulate", object, nsim, PACKAGE = "Matrix")
467                }            sc <- 1
468                new("VarCorr",            if (object@useScale)
469                    scale = .Call("lmeRep_sigma", x, REML),                sc <- .Call("mer_sigma", object, object@method == "REML",
470                    reSumry = val,                            PACKAGE = "Matrix")
471                    useScale = useScale)            ## 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  setMethod("gradient", signature(x = "lmeRep"),            ## save the seed
475            function(x, REML, unconst, ...)            attr(lpred, "seed") <- RNGstate
476            .Call("lmeRep_gradient", x, REML, unconst))            lpred

setMethod("summary", "lmeRep",
function(object, REML = TRUE, useScale = TRUE, ...) {
fcoef <- fixef(object)
corF <- as(as(vcov(object, REML, useScale), "pdmatrix"),
"corrmatrix")
DF <- getFixDF(object)
coefs <- cbind(fcoef, corF@stdDev, DF)
nc <- object@nc
dimnames(coefs) <-
list(names(fcoef), c("Estimate", "Std. Error", "DF"))
new("summary.ssclme",
coefficients = as.matrix(coefs),
scale = .Call("lmeRep_sigma", object, REML),
denomDF = as.integer(DF),
REML = REML,
ngrps = unlist(lapply(object@levels, length)),
nobs = nc[length(nc)],
corFixed = corF,
VarCorr = VarCorr(object, REML, useScale),
useScale = useScale,
showCorrelation = FALSE)
477            })            })
478
479  setMethod("show", "lmeRep",
480    setMethod("show", "mer",
481            function(object) {            function(object) {
482                fcoef <- fixef(object)                vcShow <- function(varc, useScale)
483                corF <- as(as(vcov(object, REML = TRUE, useScale = TRUE), "pdmatrix"),                {
484                           "corrmatrix")                    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
518                  corF <- vcov(object)@factors\$correlation
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                new("summary.ssclme",                digits <- max(3, getOption("digits") - 2)
524                    coefficients = as.matrix(coefs),                REML <- object@method == "REML"
525                    scale = .Call("lmeRep_sigma", object, REML = TRUE),                llik <- logLik(object, REML)
526                    denomDF = as.integer(DF),                dev <- object@deviance
527                    REML = TRUE,                devc <- object@devComp
528                    ngrps = unlist(lapply(object@levels, length)),
529                    nobs = nc[length(nc)],                rdig <- 5
530                    corFixed = corF,                if (glz <- !(object@method %in% c("REML", "ML"))) {
531                    VarCorr = VarCorr(object, REML = TRUE, useScale = TRUE),                    cat(paste("Generalized linear mixed model fit using",
532                    useScale = TRUE,                              object@method, "\n"))
533                    showCorrelation = FALSE)                } else {
534                      cat("Linear mixed-effects model fit by ")
535                      cat(if(REML) "REML\n" else "maximum likelihood\n")
536                  }
537                  if (!is.null(object@call\$formula)) {
538                      cat("Formula:", deparse(object@call\$formula),"\n")
539                  }
540                  if (!is.null(object@call\$data)) {
541                      cat("   Data:", deparse(object@call\$data), "\n")
542                  }
543                  if (!is.null(object@call\$subset)) {
544                      cat(" Subset:",
545                          deparse(asOneSidedFormula(object@call\$subset)[[2]]),"\n")
546                  }
547                  if (glz) {
548                      cat(" Family: ", object@family\$family, "(",
550                      print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
551                                       logLik = c(llik),
552                                       deviance = -2*llik,
553                                       row.names = ""))
554                  } else {
555                      print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
556                                       logLik = c(llik),
557                                       MLdeviance = dev["ML"],
558                                       REMLdeviance = dev["REML"],
559                                       row.names = ""))
560                  }
561                  cat("Random effects:\n")
562                  vcShow(VarCorr(object), useScale)
563                  ngrps <- lapply(object@flist, function(x) length(levels(x)))
564                  cat(sprintf("# of obs: %d, groups: ", devc[1]))
565                  cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
566                  cat("\n")
567                  if (!useScale)
568                      cat("\nEstimated scale (compare to 1) ",
569                          .Call("mer_sigma", object, FALSE, PACKAGE = "Matrix"),
570                          "\n")
571                  if (nrow(coefs) > 0) {
572                      if (useScale) {
573                          stat <- coefs[,1]/coefs[,2]
574                          pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
575                          nms <- colnames(coefs)
576                          coefs <- cbind(coefs, stat, pval)
577                          colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")
578                      } else {
579                          coefs <- coefs[, 1:2, drop = FALSE]
580                          stat <- coefs[,1]/coefs[,2]
581                          pval <- 2*pnorm(abs(stat), lower = FALSE)
582                          nms <- colnames(coefs)
583                          coefs <- cbind(coefs, stat, pval)
584                          colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
585                      }
586                      cat("\nFixed effects:\n")
587                      printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
588                      rn <- rownames(coefs)
589                      if (!is.null(corF)) {
590                          p <- ncol(corF)
591                          if (p > 1) {
592                              cat("\nCorrelation of Fixed Effects:\n")
593                              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)] <- ""
599                              print(corF[-1, -p, drop=FALSE], quote = FALSE)
600                          }
601                      }
602                  }
603                  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
# Line 131  Line 619
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="lmeRep"),  setMethod("getFixDF", signature(object="mer"),
623            function(object, ...)            function(object, ...) {
624                  devc <- object@devComp
625                  rep(as.integer(devc[1]- devc[2]), devc[2])
626              })
627
628    setMethod("logLik", signature(object="mer"),
629              function(object, REML = object@method == "REML", ...) {
630                  val <- -deviance(object, REML = REML)/2
631                  devc <- as.integer(object@devComp[1:2])
632                  attr(val, "nall") <- attr(val, "nobs") <- devc[1]
633                  attr(val, "df") <- abs(devc[2]) +
634                      length(.Call("mer_coef", object, 0, PACKAGE = "Matrix"))
635                  attr(val, "REML") <- REML
636                  class(val) <- "logLik"
637                  val
638              })
639
640    setMethod("VarCorr", signature(x = "mer"),
641              function(x, REML = x@method == "REML", useScale = x@useScale, ...)
642        {        {
643            nc <- object@nc[-seq(along = object@Omega)]            sc <- 1
644            p <- nc[1] - 1            if (useScale)
645            n <- nc[2]                sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix")
646            rep(n-p, p)            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="lmeRep"),  setMethod("anova", signature(object = "mer"),
657            function(object, ...)            function(object, ...)
658        {        {
659              mCall <- match.call(expand.dots = TRUE)
660              dots <- list(...)
661              modp <- logical(0)
662              if (length(dots))
663                  modp <- sapply(dots, inherits, "mer") | sapply(dots, inherits, "lm")
664              if (any(modp)) {              # multiple models - form table
665                  opts <- dots[!modp]
666                  mods <- c(list(object), dots[modp])
667                  names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)
668                  mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE), attr, "df"))]
669                  calls <- lapply(mods, slot, "call")
670                  data <- lapply(calls, "[[", "data")
671                  if (any(data != data[[1]])) stop("all models must be fit to the same data object")
673                  subset <- lapply(calls, "[[", "subset")
674                  if (any(subset != subset[[1]])) stop("all models must use the same subset")
675                  if (!is.null(subset[[1]]))
677                          c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
678                  llks <- lapply(mods, logLik, REML = FALSE)
679                  Df <- sapply(llks, attr, "df")
680                  llk <- unlist(llks)
681                  chisq <- 2 * pmax(0, c(NA, diff(llk)))
682                  dfChisq <- c(NA, diff(Df))
683                  val <- data.frame(Df = Df,
684                                    AIC = sapply(llks, AIC),
685                                    BIC = sapply(llks, BIC),
686                                    logLik = llk,
687                                    "Chisq" = chisq,
688                                    "Chi Df" = dfChisq,
689                                    "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
690                                    check.names = FALSE)
691                  class(val) <- c("anova", class(val))
694                        paste(names(mods),
695                              unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
696                              sep = ": "))
697                  return(val)
698              } else {
699            foo <- object            foo <- object
700            foo@status["factored"] <- FALSE            foo@status["factored"] <- FALSE
701            .Call("lmeRep_factor", foo, PACKAGE="Matrix")                .Call("mer_factor", foo, PACKAGE="Matrix")
702            dfr <- getFixDF(foo)            dfr <- getFixDF(foo)
703            ss <- foo@RXX[ , ".response"]^2                ss <- foo@rXy^2
704            ssr <- ss[[".response"]]                ssr <- exp(foo@devComp["logryy2"])
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            # FIXME: This only gives single degree of freedom tests                asgn <- foo@assign
708            ms <- ss                terms <- foo@terms
709            df <- rep(1, length(ms))                nmeffects <- attr(terms, "term.labels")
710            f <- ms/(ssr/dfr)                if ("(Intercept)" %in% names(ss))
711            P <- pf(f, df, dfr, lower.tail = FALSE)                    nmeffects <- c("(Intercept)", nmeffects)
712            table <- data.frame(df, ss, ms, dfr, f, P)                ss <- unlist(lapply(split(ss, asgn), sum))
713                  df <- unlist(lapply(split(asgn,  asgn), length))
714                  #dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
715                  ms <- ss/df
716                  #f <- ms/(ssr/dfr)
717                  #P <- pf(f, df, dfr, lower.tail = FALSE)
718                  #table <- data.frame(df, ss, ms, dfr, f, P)
719                  table <- data.frame(df, ss, ms)
720            dimnames(table) <-            dimnames(table) <-
721                list(names(ss),                    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            if (any(match(names(ss), "(Intercept)", nomatch = 0)))                         c("Df", "Sum Sq", "Mean Sq"))
724                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")
727            table            table
728              }
729        })        })
730
731    setMethod("confint", signature(object = "mer"),
732              function(object, parm, level = 0.95, ...)
733              .NotYetImplemented()
734              )
735
736    setMethod("fitted", signature(object = "mer"),
737              function(object, ...)
738              .Call("mer_fitted", object, PACKAGE = "Matrix")
739              )
740
741    setMethod("formula", signature(x = "mer"),
742              function(x, ...)
743              x@call\$formula
744              )
745
746    setMethod("residuals", signature(object = "mer"),
747              function(object, ...)
748              .NotYetImplemented()
749              )
750
751    setMethod("resid", signature(object = "mer"),
752              function(object, ...)
753              .NotYetImplemented()
754              )
755
756    setMethod("summary", signature(object = "mer"),
757              function(object, ...)
758              .NotYetImplemented()
759              )
760
761    setMethod("update", signature(object = "mer"),
762              function(object, ...)
763              .NotYetImplemented()
764              )
765
766
767    simss <- function(fm0, fma, nsim)
768    {
769        ysim <- simulate(fm0, nsim)
770        cv <- list(analyticGradient = FALSE, msMaxIter = 200:200,
771                   msVerbose = 0:0)
772        sapply(ysim, function(yy) {
773            .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")
774            LMEoptimize(fm0) <- cv
775            .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    }

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