SCM

SCM Repository

[matrix] Annotation of /branches/trunk-lme4/R/lmer.R
ViewVC logotype

Annotation of /branches/trunk-lme4/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2084 - (view) (download)

1 : bates 767 # Methods for lmer and for the objects that it produces
2 : bates 689
3 : bates 1303 ## To Do: Determine why the names of the components of the values of
4 :     ## the ranef and coef extractor methods are not printed.
5 :    
6 : bates 689 ## Some utilities
7 :    
8 : bates 775 ## Return the pairs of expressions separated by vertical bars
9 : bates 769 findbars <- function(term)
10 :     {
11 : bates 1284 if (is.name(term) || !is.language(term)) return(NULL)
12 : bates 769 if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
13 :     if (!is.call(term)) stop("term must be of class call")
14 :     if (term[[1]] == as.name('|')) return(term)
15 :     if (length(term) == 2) return(findbars(term[[2]]))
16 :     c(findbars(term[[2]]), findbars(term[[3]]))
17 :     }
18 :    
19 : bates 775 ## Return the formula omitting the pairs of expressions
20 :     ## that are separated by vertical bars
21 : bates 769 nobars <- function(term)
22 :     {
23 : bates 1150 if (!('|' %in% all.names(term))) return(term)
24 : bates 769 if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
25 :     if (length(term) == 2) {
26 : maechler 1222 nb <- nobars(term[[2]])
27 :     if (is.null(nb)) return(NULL)
28 :     term[[2]] <- nb
29 :     return(term)
30 : bates 769 }
31 :     nb2 <- nobars(term[[2]])
32 :     nb3 <- nobars(term[[3]])
33 :     if (is.null(nb2)) return(nb3)
34 :     if (is.null(nb3)) return(nb2)
35 :     term[[2]] <- nb2
36 :     term[[3]] <- nb3
37 :     term
38 :     }
39 :    
40 : bates 775 ## Substitute the '+' function for the '|' function
41 : bates 769 subbars <- function(term)
42 :     {
43 : bates 1284 if (is.name(term) || !is.language(term)) return(term)
44 : bates 769 if (length(term) == 2) {
45 : maechler 1222 term[[2]] <- subbars(term[[2]])
46 :     return(term)
47 : bates 769 }
48 : bates 1607 stopifnot(length(term) >= 3)
49 : maechler 832 if (is.call(term) && term[[1]] == as.name('|'))
50 : maechler 1222 term[[1]] <- as.name('+')
51 : bates 1607 for (j in 2:length(term)) term[[j]] <- subbars(term[[j]])
52 : bates 769 term
53 :     }
54 : bates 824
55 : bates 1856 ## Substitute any names from nlist in term with 1
56 :     subnms <- function(term, nlist)
57 :     {
58 :     if (!is.language(term)) return(term)
59 :     if (is.name(term)) {
60 :     if (any(unlist(lapply(nlist, get("=="), term)))) return(1)
61 :     return(term)
62 :     }
63 :     stopifnot(length(term) >= 2)
64 :     for (j in 2:length(term)) term[[j]] <- subnms(term[[j]], nlist)
65 :     term
66 :     }
67 :    
68 : bates 1286 ## Return the list of '/'-separated terms in an expression that
69 :     ## contains slashes
70 :     slashTerms <- function(x) {
71 :     if (!("/" %in% all.names(x))) return(x)
72 :     if (x[[1]] != as.name("/"))
73 :     stop("unparseable formula for grouping factor")
74 :     list(slashTerms(x[[2]]), slashTerms(x[[3]]))
75 : bates 979 }
76 :    
77 : bates 1286 ## from a list of length 2 return recursive interaction terms
78 :     makeInteraction <- function(x) {
79 :     if (length(x) < 2) return(x)
80 :     trm1 <- makeInteraction(x[[1]])
81 :     trm11 <- if(is.list(trm1)) trm1[[1]] else trm1
82 :     list(substitute(foo:bar, list(foo=x[[2]], bar = trm11)), trm1)
83 :     }
84 :    
85 : maechler 1343
86 :     factorNames2char <- function(nms, collapse = ", ") {
87 :     ## utility in messages / print etc:
88 :     nms <- sQuote(nms)
89 :     if(length(nms) == 1) paste("factor", nms)
90 :     else paste("factors", paste(nms, collapse = collapse))
91 :     }
92 :    
93 : bates 1286 ## expand any slashes in the grouping factors returned by findbars
94 :     expandSlash <- function(bb) {
95 : bates 1561 if (!is.list(bb)) return(expandSlash(list(bb)))
96 : bates 1286 ## I really do mean lapply(unlist(... - unlist returns a
97 :     ## flattened list in this case
98 :     unlist(lapply(bb, function(x) {
99 : bates 1561 if (length(x) > 2 && is.list(trms <- slashTerms(x[[3]])))
100 : bates 1286 return(lapply(unlist(makeInteraction(trms)),
101 :     function(trm) substitute(foo|bar,
102 :     list(foo = x[[2]],
103 :     bar = trm))))
104 :     x
105 :     }))
106 :     }
107 : maechler 1299
108 : bates 824 abbrvNms <- function(gnm, cnms)
109 :     {
110 :     ans <- paste(abbreviate(gnm), abbreviate(cnms), sep = '.')
111 :     if (length(cnms) > 1) {
112 : maechler 1222 anms <- lapply(cnms, abbreviate, minlength = 3)
113 :     nmmat <- outer(anms, anms, paste, sep = '.')
114 :     ans <- c(ans, paste(abbreviate(gnm, minlength = 3),
115 :     nmmat[upper.tri(nmmat)], sep = '.'))
116 : bates 824 }
117 :     ans
118 :     }
119 :    
120 : bates 775 ## Control parameters for lmer
121 :     lmerControl <-
122 : bates 1605 function(maxIter = 200, # used in ../src/lmer.c only
123 :     tolerance = sqrt(.Machine$double.eps),# ditto
124 :     msMaxIter = 200,
125 :     ## msTol = sqrt(.Machine$double.eps),
126 :     ## FIXME: should be able to pass tolerances to nlminb()
127 :     msVerbose = getOption("verbose"),
128 :     niterEM = 15,
129 :     EMverbose = getOption("verbose"),
130 :     PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead
131 : bates 1683 usePQL = FALSE,
132 : bates 1605 gradient = TRUE,
133 :     Hessian = FALSE # unused _FIXME_
134 :     )
135 : bates 435 {
136 : bates 775 list(maxIter = as.integer(maxIter),
137 : maechler 1222 tolerance = as.double(tolerance),
138 :     msMaxIter = as.integer(msMaxIter),
139 :     ## msTol = as.double(msTol),
140 :     msVerbose = as.integer(msVerbose),# "integer" on purpose
141 :     niterEM = as.integer(niterEM),
142 :     EMverbose = as.logical(EMverbose),
143 :     PQLmaxIt = as.integer(PQLmaxIt),
144 :     usePQL = as.logical(usePQL),
145 :     gradient = as.logical(gradient),
146 :     Hessian = as.logical(Hessian))
147 : bates 435 }
148 :    
149 : bates 1630 ## check for predefined families
150 :     mkFltype <- function(family)
151 :     {
152 :     fltype <- 0 # not a predefined type
153 :     if (family$family == "gaussian" && family$link == "identity") fltype <- -1
154 :     if (family$family == "binomial" && family$link == "logit") fltype <- 1
155 :     if (family$family == "binomial" && family$link == "probit") fltype <- 2
156 : bates 1637 if (family$family == "poisson" && family$link == "log") fltype <- 3
157 : bates 1630 as.integer(fltype)
158 :     }
159 :    
160 : bates 1166 rWishart <- function(n, df, invScal)
161 : bates 1504 .Call(lme4_rWishart, n, df, invScal)
162 : bates 1166
163 :     setMethod("coef", signature(object = "mer"),
164 : maechler 1222 function(object, ...)
165 : bates 1150 {
166 : bates 1240 if (length(list(...)))
167 :     warning(paste('arguments named "',
168 :     paste(names(list(...)), collapse = ", "),
169 : bates 1445 '" ignored', sep = ''))
170 : bates 1236 fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
171 :     ref <- ranef(object)
172 :     val <- lapply(ref, function(x) fef[rep(1, nrow(x)),,drop = FALSE])
173 :     for (i in seq(a = val)) {
174 :     refi <- ref[[i]]
175 :     row.names(val[[i]]) <- row.names(refi)
176 :     nmsi <- colnames(refi)
177 :     if (!all(nmsi %in% names(fef)))
178 :     stop("unable to align random and fixed effects")
179 :     for (nm in nmsi) val[[i]][[nm]] <- val[[i]][[nm]] + refi[,nm]
180 :     }
181 :     new("coef.lmer", val)
182 : bates 1166 })
183 : bates 1150
184 : bates 1234 setMethod("plot", signature(x = "coef.lmer"),
185 :     function(x, y, ...)
186 :     {
187 :     varying <- unique(do.call("c",
188 :     lapply(x, function(el)
189 :     names(el)[sapply(el,
190 :     function(col)
191 :     any(col != col[1]))])))
192 : bates 1811 gf <- do.call("rBind", lapply(x, "[", j = varying))
193 : bates 1234 gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
194 :     switch(min(length(varying), 3),
195 :     qqmath(eval(substitute(~ x | .grp,
196 :     list(x = as.name(varying[1])))), gf, ...),
197 :     xyplot(eval(substitute(y ~ x | .grp,
198 :     list(y = as.name(varying[1]),
199 :     x = as.name(varying[2])))), gf, ...),
200 :     splom(~ gf | .grp, ...))
201 :     })
202 : bates 1150
203 : bates 1234 setMethod("plot", signature(x = "ranef.lmer"),
204 : maechler 1222 function(x, y, ...)
205 : bates 1166 {
206 : maechler 1222 lapply(x, function(x) {
207 :     cn <- lapply(colnames(x), as.name)
208 :     switch(min(ncol(x), 3),
209 :     qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
210 :     xyplot(eval(substitute(y ~ x,
211 :     list(y = cn[[1]],
212 :     x = cn[[2]]))), x, ...),
213 :     splom(~ x, ...))
214 :     })
215 : bates 1166 })
216 : bates 1150
217 :     setMethod("with", signature(data = "lmer"),
218 : maechler 1222 function(data, expr, ...) {
219 :     dat <- eval(data@call$data)
220 :     if (!is.null(na.act <- attr(data@frame, "na.action")))
221 :     dat <- dat[-na.act, ]
222 :     lst <- c(list(. = data), data@flist, data@frame, dat)
223 :     eval(substitute(expr), lst[unique(names(lst))])
224 :     })
225 : bates 1150
226 :     setMethod("terms", signature(x = "lmer"),
227 : maechler 1222 function(x, ...) x@terms)
228 : bates 1150
229 : bates 1605 ## Utility functions used in lmer and variations on lmer
230 :    
231 :     ## Check that the 'start' argument matches the form of the Omega
232 :     ## slot in mer. If so, install start as the Omega slot.
233 : bates 1590 setOmega <- function(mer, start)
234 : bates 1586 {
235 : bates 1590 Om <- mer@Omega
236 :     if (!is.list(start) || length(start) != length(Om) ||
237 :     !all.equal(names(Om), names(start)))
238 :     stop(paste("start must be a list of length", length(Om),
239 :     "with names\n", paste(names(Om), collapse = ',')))
240 :     for (i in seq(along = start)) {
241 :     if (class(start[[i]]) != class(Om[[i]]) ||
242 :     !all.equal(dim(start[[i]]), dim(Om[[i]])))
243 :     stop(paste("start[[", i, "]] must be of class '", class(Om[[i]]),
244 :     "' and dimension ", paste(dim(Om[[i]]), collapse = ','),
245 :     sep = ''))
246 :     }
247 :     mer@Omega <- start
248 :     mer
249 : bates 1586 }
250 : bates 1150
251 : bates 1605 ## Create model frame, terms, weights and offset from an lmer formula
252 :     ##
253 : bates 1590 ## mc - matched call to parent function
254 :     ## formula - two-sided formula
255 :     ## contrasts - contrasts argument
256 : bates 1856 ## vnms - names of variables to be included in the model frame
257 :     lmerFrames <- function(mc, formula, contrasts, vnms = character(0))
258 : bates 1590 {
259 :     ## Must evaluate the model frame first and then fit the glm using
260 :     ## that frame. Otherwise missing values in the grouping factors
261 :     ## cause inconsistent numbers of observations.
262 :     mf <- mc
263 :     m <- match(c("data", "subset", "weights", "na.action", "offset"),
264 :     names(mf), 0)
265 :     mf <- mf[c(1, m)]
266 : bates 1778 frame.form <- subbars(formula) # substitute `+' for `|'
267 : bates 1856 if (length(vnms) > 0)
268 :     frame.form[[3]] <-
269 :     substitute(foo + bar,
270 :     list(foo = parse(text = paste(vnms, collapse = ' + '))[[1]],
271 :     bar = frame.form[[3]]))
272 : bates 1778 fixed.form <- nobars(formula) # remove any terms with `|'
273 : bates 1590 if (inherits(fixed.form, "name")) # RHS is empty - use a constant
274 :     fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
275 :     environment(fixed.form) <- environment(frame.form) <- environment(formula)
276 : bates 1150
277 : bates 1590 ## evaluate a model frame for fixed and random effects
278 :     mf$formula <- frame.form
279 :     mf$drop.unused.levels <- TRUE
280 :     mf[[1]] <- as.name("model.frame")
281 :     fe <- mf
282 : bates 1605 mf <- eval(mf, parent.frame(2))
283 : bates 1150
284 : bates 1590 ## get the terms for the fixed-effects only
285 :     fe$formula <- fixed.form
286 : bates 1605 fe <- eval(fe, parent.frame(2))
287 : bates 1590 mt <- attr(fe, "terms") # allow model.frame to update them
288 :     ## response vector
289 : bates 1683 Y <- model.response(mf, "any")
290 : bates 1590 ## avoid problems with 1D arrays, but keep names
291 :     if(length(dim(Y)) == 1) {
292 :     nm <- rownames(Y)
293 :     dim(Y) <- NULL
294 :     if(!is.null(nm)) names(Y) <- nm
295 :     }
296 :     ## null model support
297 :     X <- if (!is.empty.model(mt))
298 :     model.matrix(mt, mf, contrasts) else matrix(,NROW(Y),0)
299 : maechler 1188
300 : bates 1590 weights <- model.weights(mf)
301 :     offset <- model.offset(mf)
302 :     ## check weights and offset
303 :     if( !is.null(weights) && any(weights < 0) )
304 :     stop("negative weights not allowed")
305 :     if(!is.null(offset) && length(offset) != NROW(Y))
306 :     stop(gettextf("number of offsets is %d should equal %d (number of observations)",
307 :     length(offset), NROW(Y)), domain = NA)
308 :     list(Y = Y, X = X, weights = weights, offset = offset, mt = mt, mf = mf)
309 :     }
310 : maechler 1188
311 : bates 1605 ## Create the list of grouping factors and corresponding model
312 :     ## matrices. The model matrices are transposed in the Ztl list
313 : bates 1856
314 :     ## FIXME: Should the check on fltype be in this function or an auxillary check?
315 : bates 1683 lmerFactorList <- function(formula, mf, fltype)
316 : bates 1590 {
317 :     ## create factor list for the random effects
318 :     bars <- expandSlash(findbars(formula[[3]]))
319 : bates 1645 if (!length(bars)) stop("No random effects terms specified in formula")
320 : bates 1590 names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
321 :     fl <- lapply(bars,
322 :     function(x)
323 :     eval(substitute(as.factor(fac)[,drop = TRUE],
324 :     list(fac = x[[3]])), mf))
325 :     ## order factor list by decreasing number of levels
326 :     nlev <- sapply(fl, function(x) length(levels(x)))
327 : bates 1605 nobs <- length(fl[[1]])
328 : bates 1590 if(any(nlev == 0))
329 :     stop("resulting factor(s) with 0 levels in random effects part:\n ",
330 :     paste(sQuote(names(nlev[nlev == 0])), collapse=", "))
331 : bates 1683 ## Max # of levels allowed for a grouping factor.
332 :     ## Binomial glmms can have nlev == nobs.
333 : bates 1912 maxlev <- nobs - !(fltype %in% 1:3)
334 : bates 1683 if (any(nlev > maxlev))
335 :     stop("number of levels in grouping factor(s) ",
336 :     paste(sQuote(names(nlev[nlev > maxlev])), collapse=", "),
337 :     " is too large")
338 : bates 1590 if (any(diff(nlev) > 0)) {
339 :     ord <- rev(order(nlev))
340 :     bars <- bars[ord]
341 :     fl <- fl[ord]
342 :     }
343 :     ## create list of transposed model matrices for random effects
344 :     Ztl <- lapply(bars, function(x)
345 :     t(model.matrix(eval(substitute(~ expr,
346 :     list(expr = x[[2]]))),
347 :     mf)))
348 : bates 2084 if (!all(sapply(Ztl, nrow)))
349 :     stop("terms that reduce to (0|grp) are not allowed")
350 : bates 1605 list(fl = fl, Ztl = Ztl)
351 : bates 1590 }
352 : bates 755
353 : bates 1590 lmer <- function(formula, data, family = gaussian,
354 : bates 1605 method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
355 : bates 1590 control = list(), start = NULL,
356 :     subset, weights, na.action, offset, contrasts = NULL,
357 :     model = TRUE, ...)
358 :     {
359 : bates 1683 method <- match.arg(method)
360 : bates 1590 formula <- as.formula(formula)
361 : bates 1607 if (length(formula) < 3) stop("formula must be a two-sided formula")
362 :     cv <- do.call("lmerControl", control)
363 : bates 1166
364 : bates 1590 ## Establish model frame and fixed-effects model matrix and terms
365 :     mc <- match.call()
366 : bates 1856 fr <- lmerFrames(mc, formula, contrasts)
367 : bates 1590 Y <- fr$Y; X <- fr$X; weights <- fr$weights; offset <- fr$offset
368 :     mf <- fr$mf; mt <- fr$mt
369 : bates 1605
370 : bates 1590 ## check and evaluate the family argument
371 : bates 1800 if(is.character(family))
372 :     family <- get(family, mode = "function", envir = parent.frame(2))
373 :     if(is.function(family)) family <- family()
374 :     if(is.null(family$family)) stop("'family' not recognized")
375 : bates 1630 fltype <- mkFltype(family)
376 : maechler 1677
377 : bates 1683 ## establish factor list and Ztl
378 :     FL <- lmerFactorList(formula, mf, fltype)
379 :     cnames <- with(FL, c(lapply(Ztl, rownames), list(.fixed = colnames(X))))
380 :     nc <- with(FL, sapply(Ztl, nrow))
381 :     Ztl <- with(FL, .Call(Ztl_sparse, fl, Ztl))
382 :     ## FIXME: change this when rbind has been fixed.
383 : bates 1811 Zt <- if (length(Ztl) == 1) Ztl[[1]] else do.call("rBind", Ztl)
384 : bates 1683 fl <- FL$fl
385 :    
386 : bates 1590 ## quick return for a linear mixed model
387 : bates 1627 if (fltype < 0) {
388 : bates 1683 mer <- .Call(mer_create, fl, Zt, X, as.double(Y),
389 :     method == "REML", nc, cnames)
390 : bates 1590 if (!is.null(start)) mer <- setOmega(mer, start)
391 :     .Call(mer_ECMEsteps, mer, cv$niterEM, cv$EMverbose)
392 :     LMEoptimize(mer) <- cv
393 :     return(new("lmer", mer,
394 :     frame = if (model) fr$mf else data.frame(),
395 :     terms = mt,
396 :     call = mc))
397 :     }
398 : bates 1150
399 : bates 1590 ## The rest of the function applies to generalized linear mixed models
400 : bates 1605 if (method %in% c("ML", "REML")) method <- "Laplace"
401 : bates 1590 if (method == "AGQ")
402 :     stop('method = "AGQ" not yet implemented for supernodal representation')
403 : bates 1605 if (method == "PQL") cv$usePQL <- TRUE # have to use PQL for method == "PQL"
404 : maechler 1677
405 : bates 1590 ## initial fit of a glm to the fixed-effects only.
406 :     glmFit <- glm.fit(X, Y, weights = weights, offset = offset, family = family,
407 :     intercept = attr(mt, "intercept") > 0)
408 :     Y <- as.double(glmFit$y)
409 :     ## must do this after Y has possibly been reformulated
410 : bates 1605 mer <- .Call(mer_create, fl, Zt, X, Y, 0, nc, cnames)
411 : bates 1590 if (!is.null(start)) mer <- setOmega(mer, start)
412 :    
413 : bates 1605 gVerb <- getOption("verbose")
414 :    
415 : bates 1590 ## extract some of the components of glmFit
416 :     ## weights could have changed
417 :     weights <- glmFit$prior.weights
418 : bates 1800 if (is.null(offset)) offset <- numeric(length(Y))
419 : bates 1590 eta <- glmFit$linear.predictors
420 :     linkinv <- quote(family$linkinv(eta))
421 :     mu.eta <- quote(family$mu.eta(eta))
422 :     mu <- family$linkinv(eta)
423 :     variance <- quote(family$variance(mu))
424 : bates 1636 dev.resids <- quote(family$dev.resids(Y, mu, weights))
425 : bates 1590 LMEopt <- get("LMEoptimize<-")
426 :     doLMEopt <- quote(LMEopt(x = mer, value = cv))
427 :     if (family$family %in% c("binomial", "poisson")) # set the constant scale
428 : bates 1637 mer@devComp[8] <- -1
429 : bates 1605 mer@status["glmm"] <- as.integer(switch(method, PQL = 1, Laplace = 2, AGQ = 3))
430 : bates 1627 GSpt <- .Call(glmer_init, environment(), fltype)
431 : bates 1590 if (cv$usePQL) {
432 :     .Call(glmer_PQL, GSpt) # obtain PQL estimates
433 :     PQLpars <- c(fixef(mer),
434 :     .Call(mer_coef, mer, 2))
435 :     } else {
436 :     PQLpars <- c(coef(glmFit),
437 :     .Call(mer_coef, mer, 2))
438 :     }
439 :     if (method == "PQL") {
440 :     .Call(glmer_devLaplace, PQLpars, GSpt)
441 :     .Call(glmer_finalize, GSpt)
442 :     return(new("glmer",
443 :     new("lmer", mer,
444 :     frame = if (model) mf else data.frame(),
445 : bates 1605 terms = mt, call = match.call()),
446 : bates 1590 weights = weights,
447 :     family=family))
448 :     }
449 :    
450 :     fixInd <- seq(ncol(X))
451 :     ## pars[fixInd] == beta, pars[-fixInd] == theta
452 :     ## indicator of constrained parameters
453 :     const <- c(rep(FALSE, length(fixInd)),
454 :     unlist(lapply(mer@nc[seq(along = fl)],
455 :     function(k) 1:((k*(k+1))/2) <= k)
456 :     ))
457 : bates 1605 devLaplace <- function(pars) .Call(glmer_devLaplace, pars, GSpt)
458 : bates 1668 rel.tol <- abs(0.01/devLaplace(PQLpars))
459 : bates 1683 if (cv$msVerbose) cat(paste("relative tolerance set to", rel.tol, "\n"))
460 : bates 1590
461 : bates 1605 optimRes <- nlminb(PQLpars, devLaplace,
462 :     lower = ifelse(const, 5e-10, -Inf),
463 :     control = list(trace = cv$msVerbose,
464 : bates 1668 iter.max = cv$msMaxIter,
465 :     rel.tol = rel.tol))
466 : bates 1590 .Call(glmer_finalize, GSpt)
467 :     new("glmer",
468 :     new("lmer", mer,
469 :     frame = if (model) mf else data.frame(),
470 : bates 1605 terms = mt, call = match.call()),
471 : bates 1590 weights = weights,
472 :     family=family)
473 : bates 1605
474 : bates 1590 }
475 :    
476 : bates 1150 ## Extract the L matrix
477 :     setAs("mer", "dtCMatrix", function(from)
478 : maechler 1280 .Call(mer_dtCMatrix, from))
479 : bates 1150
480 :     ## Extract the fixed effects
481 :     setMethod("fixef", signature(object = "mer"),
482 : maechler 1222 function(object, ...)
483 : maechler 1280 .Call(mer_fixef, object))
484 : bates 1150
485 :     ## Extract the random effects
486 :     setMethod("ranef", signature(object = "mer"),
487 : bates 1243 function(object, postVar = FALSE, ...) {
488 : bates 1234 ans <- new("ranef.lmer",
489 : maechler 1280 lapply(.Call(mer_ranef, object),
490 : bates 1232 data.frame, check.names = FALSE))
491 :     names(ans) <- names(object@flist)
492 : bates 1243 if (postVar) {
493 : maechler 1280 pV <- .Call(mer_postVar, object)
494 : bates 1243 for (i in seq(along = ans))
495 :     attr(ans[[i]], "postVar") <- pV[[i]]
496 :     }
497 : bates 1232 ans
498 :     })
499 : bates 1709
500 : bates 1150 ## Optimization for mer objects
501 : bates 755 setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
502 : maechler 1222 function(x, value)
503 :     {
504 :     if (value$msMaxIter < 1) return(x)
505 :     nc <- x@nc
506 :     constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
507 :     fn <- function(pars)
508 : maechler 1280 deviance(.Call(mer_coefGets, x, pars, 2))
509 : bates 1636 start <- .Call(mer_coef, x, 2) #starting values
510 :     fval <- fn(start)
511 : maechler 1222 gr <- if (value$gradient)
512 :     function(pars) {
513 :     if (!isTRUE(all.equal(pars,
514 : maechler 1280 .Call(mer_coef, x,
515 :     2))))
516 :     .Call(mer_coefGets, x, pars, 2)
517 :     .Call(mer_gradient, x, 2)
518 : maechler 1222 }
519 :     else NULL
520 : maechler 1280 optimRes <- nlminb(.Call(mer_coef, x, 2),
521 : maechler 1222 fn, gr,
522 :     lower = ifelse(constr, 5e-10, -Inf),
523 :     control = list(iter.max = value$msMaxIter,
524 : bates 1636 trace = as.integer(value$msVerbose),
525 :     rel.tol = abs(0.001/fval)))
526 : bates 1342 estPar <- optimRes$par
527 :     .Call(mer_coefGets, x, estPar, 2)
528 : maechler 1343
529 :     ## check for convergence on boundary
530 :     if (any(bd <- (estPar[constr] < 1e-9))) {
531 :     bpar <- rep.int(FALSE, length(estPar))
532 :     bpar[constr] <- bd
533 :     bgrp <- split(bpar,
534 :     rep(seq(along = nc),
535 :     unlist(lapply(nc,
536 :     function(k) (k*(k+1))/2))))
537 :     bdd <- unlist(lapply(bgrp, any))
538 :     lens <- unlist(lapply(bgrp, length))
539 :     if (all(lens[bdd] == 1)) { # variance components only
540 :     warning("Estimated variance for ",
541 :     factorNames2char(names(x@flist)[bdd]),
542 :     " is effectively zero\n")
543 :     } else {
544 :     warning("Estimated variance-covariance for ",
545 :     factorNames2char(names(x@flist)[bdd]),
546 :     " is singular\n")
547 :     }
548 :     }
549 : maechler 1222 if (optimRes$convergence != 0) {
550 : maechler 1343 warning("nlminb returned message ", optimRes$message,"\n")
551 : maechler 1222 }
552 :     return(x)
553 :     })
554 : bates 316
555 : bates 1241 setMethod("qqmath", signature(x = "ranef.lmer"),
556 :     function(x, data, ...) {
557 : bates 1243 prepanel.ci <- function(x, y, se, subscripts, ...) {
558 :     y <- as.numeric(y)
559 :     se <- as.numeric(se[subscripts])
560 :     hw <- 1.96 * se
561 :     list(ylim = range(y - hw, y + hw, finite = TRUE))
562 :     }
563 :     panel.ci <- function(x, y, se, subscripts, pch = 16, ...) {
564 : bates 1288 panel.grid(h = -1,v = -1)
565 :     panel.abline(h = 0)
566 : bates 1243 x <- as.numeric(x)
567 :     y <- as.numeric(y)
568 :     se <- as.numeric(se[subscripts])
569 :     ly <- y - 1.96 * se
570 : maechler 1280 uy <- y + 1.96 * se
571 : bates 1243 panel.segments(x, y - 1.96*se, x, y + 1.96 * se,
572 :     col = 'black')
573 :     panel.xyplot(x, y, pch = pch, ...)
574 :     }
575 :     f <- function(x) {
576 : bates 1445 if (!is.null(pv <- attr(x, "postVar"))) {
577 : bates 1243 cols <- 1:(dim(pv)[1])
578 :     se <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ])))
579 :     nr <- nrow(x)
580 :     nc <- ncol(x)
581 :     ord <- unlist(lapply(x, order)) +
582 :     rep((0:(nc - 1)) * nr, each = nr)
583 :     rr <- 1:nr
584 :     ind <- gl(ncol(x), nrow(x), labels = names(x))
585 :     xyplot(unlist(x)[ord] ~
586 :     rep(qnorm((rr - 0.5)/nr), ncol(x)) | ind[ord],
587 : bates 1445 se = se[ord], prepanel = prepanel.ci, panel = panel.ci,
588 : bates 1243 scales = list(y = list(relation = "free")),
589 :     xlab = "Standard normal quantiles",
590 :     ylab = NULL, aspect = 1, ...)
591 :     } else {
592 :     qqmath(~values|ind, stack(x),
593 :     scales = list(y = list(relation = "free")),
594 :     xlab = "Standard normal quantiles",
595 :     ylab = NULL, ...)
596 :     }
597 :     }
598 : bates 1241 lapply(x, f)
599 :     })
600 : maechler 1280
601 : bates 1150 setMethod("deviance", signature(object = "mer"),
602 : bates 1303 function(object, REML = NULL, ...) {
603 : bates 1588 if (is.null(REML)) REML <- object@status["REML"]
604 : maechler 1280 .Call(mer_factor, object)
605 : bates 1303 object@deviance[[ifelse(REML, "REML", "ML")]]
606 : maechler 1222 })
607 : bates 316
608 : bates 1306 ## Mangle the names of the columns of the mcmcsamp result ans
609 :     ## This operation is common to the methods for "lmer" and "glmer"
610 : bates 1720 mcmccompnames <- function(ans, object, saveb, trans, glmer, deviance)
611 : bates 1306 {
612 :     gnms <- names(object@flist)
613 :     cnms <- object@cnames
614 :     ff <- fixef(object)
615 :     colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
616 :     unlist(lapply(seq(along = gnms),
617 :     function(i)
618 :     abbrvNms(gnms[i],cnms[[i]]))))
619 :     if (trans) {
620 :     ## parameter type: 0 => fixed effect, 1 => variance,
621 :     ## 2 => covariance
622 :     ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1,
623 :     unlist(lapply(seq(along = gnms),
624 :     function(i)
625 :     {
626 :     k <- length(cnms[[i]])
627 :     rep(1:2, c(k, (k*(k-1))/2))
628 :     })))
629 :     colnms[ptyp == 1] <-
630 :     paste("log(", colnms[ptyp == 1], ")", sep = "")
631 :     colnms[ptyp == 2] <-
632 :     paste("atanh(", colnms[ptyp == 2], ")", sep = "")
633 :     }
634 : bates 1720 if (deviance) colnms <- c(colnms, "deviance")
635 : bates 1784 ## FIXME: this will fail for a mer2 object
636 : bates 1306 if(saveb) {## maybe better colnames, "RE.1","RE.2", ... ?
637 :     rZy <- object@rZy
638 :     colnms <- c(colnms,
639 :     paste("b", sprintf(paste("%0",
640 :     1+floor(log(length(rZy),10)),
641 :     "d", sep = ''),
642 :     seq(along = rZy)),
643 :     sep = '.'))
644 :     }
645 :     colnames(ans) <- colnms
646 :     ans
647 :     }
648 : maechler 1330
649 : bates 1306 setMethod("mcmcsamp", signature(object = "lmer"),
650 : maechler 1222 function(object, n = 1, verbose = FALSE, saveb = FALSE,
651 : bates 1720 trans = TRUE, deviance = FALSE, ...)
652 : bates 1150 {
653 : bates 1720 ans <- t(.Call(mer_MCMCsamp, object, saveb, n, trans, verbose, deviance))
654 : maechler 1222 attr(ans, "mcpar") <- as.integer(c(1, n, 1))
655 :     class(ans) <- "mcmc"
656 : maechler 1730 mcmccompnames(ans, object, saveb, trans,
657 :     glmer=FALSE, deviance=deviance)
658 : bates 1150 })
659 : bates 316
660 : bates 1306 setMethod("mcmcsamp", signature(object = "glmer"),
661 :     function(object, n = 1, verbose = FALSE, saveb = FALSE,
662 : maechler 1730 trans = TRUE, deviance = FALSE, ...)
663 : bates 1306 {
664 :     family <- object@family
665 :     mer <- as(object, "mer")
666 :     weights <- object@weights
667 :     cv <- lmerControl()
668 :     eta <- .Call(mer_fitted, mer)
669 :     offset <- numeric(length(eta)) ## change this, save the offset in mer
670 :     Y <- object@y
671 :     linkinv <- quote(family$linkinv(eta))
672 :     mu.eta <- quote(family$mu.eta(eta))
673 :     mu <- family$linkinv(eta)
674 :     variance <- quote(family$variance(mu))
675 : bates 1637 dev.resids <- quote(family$dev.resids(Y, mu, weights))
676 : bates 1306 LMEopt <- get("LMEoptimize<-")
677 :     doLMEopt <- quote(LMEopt(x = mer, value = cv))
678 : bates 1630 fltype <- mkFltype(family)
679 :     GSpt <- .Call(glmer_init, environment(), fltype)
680 : bates 1743 ans <- t(.Call(glmer_MCMCsamp, GSpt, saveb, n, trans, verbose, deviance))
681 : bates 1306 .Call(glmer_finalize, GSpt)
682 :     attr(ans, "mcpar") <- as.integer(c(1, n, 1))
683 :     class(ans) <- "mcmc"
684 : maechler 1730 mcmccompnames(ans, object, saveb, trans,
685 : bates 1743 glmer=TRUE, deviance=deviance)
686 : bates 1306 })
687 :    
688 : bates 1150 setMethod("simulate", signature(object = "mer"),
689 : bates 2029 function(object, nsim = 1, seed = NULL, ...)
690 : bates 1150 {
691 : bates 2029 if(!exists(".Random.seed", envir = .GlobalEnv))
692 :     runif(1) # initialize the RNG if necessary
693 :     if(is.null(seed))
694 :     RNGstate <- .Random.seed
695 :     else {
696 :     R.seed <- .Random.seed
697 :     set.seed(seed)
698 :     RNGstate <- structure(seed, kind = as.list(RNGkind()))
699 :     on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
700 :     }
701 :    
702 :     stopifnot((nsim <- as.integer(nsim[1])) > 0, is(object, "lmer"))
703 :     ## similate the linear predictors
704 :     lpred <- .Call(lme4:::mer_simulate, object, nsim)
705 :     sc <- abs(object@devComp[8])
706 :    
707 :     ## add fixed-effects contribution
708 :     lpred <- lpred + drop(object@X %*% fixef(object))
709 :     n <- prod(dim(lpred))
710 :    
711 :     if (!is(object, "glmer")) { # and per-observation noise term
712 :     response <- as.data.frame(lpred + rnorm(n, sd = sc))
713 :     attr(response, "seed") <- RNGstate
714 :     return(response)
715 :     }
716 :     fam <- attr(object, "family")
717 :     if ("poisson"==fam$family) {
718 :     response <- as.data.frame(matrix(byrow = FALSE, ncol=nsim,
719 :     rpois(n, fam$linkinv(lpred))))
720 :     attr(response, "seed") <- RNGstate
721 :     return(response)
722 :     }
723 :     stop("simulate() not yet implemented for", fam$family,
724 :     "glmms\n")
725 : bates 1150 })
726 : bates 316
727 : maechler 1677 ### "FIXME": the following has too(?) much cut & paste from above
728 :    
729 : bates 1297 simulestimate <- function(x, FUN, nsim = 1, seed = NULL, control = list())
730 :     {
731 :     FUN <- match.fun(FUN)
732 : maechler 1677 stopifnot((nsim <- as.integer(nsim[1])) > 0,
733 : maechler 2026 is(x, "lmer"))
734 : bates 1336 if (!is.null(seed)) set.seed(seed)
735 : bates 1809 ## simulate the linear predictors
736 : bates 1297 lpred <- .Call(mer_simulate, x, nsim)
737 : bates 1587 sc <- abs(x@devComp[8])
738 : bates 1297 ## add fixed-effects contribution and per-observation noise term
739 :     lpred <- lpred + drop(x@X %*% fixef(x)) + rnorm(prod(dim(lpred)), sd = sc)
740 : bates 316
741 : bates 1297 cv <- do.call(lmerControl, control)
742 :     Omega <- x@Omega
743 :     x@wrkres <- x@y <- lpred[,1]
744 :     .Call(mer_update_ZXy, x)
745 :     LMEoptimize(x) <- cv
746 :     template <- FUN(x)
747 :     if (!is.numeric(template))
748 :     stop("simulestimate currently only handles functions that return numeric vectors")
749 :     ans <- matrix(template, nr = nsim, nc = length(template), byrow = TRUE)
750 :     colnames(ans) <- names(template)
751 :     for (i in 1:nsim) {
752 :     x@wrkres <- x@y <- lpred[,i]
753 :     x@Omega <- Omega
754 :     .Call(mer_update_ZXy, x)
755 :     LMEoptimize(x) <- cv
756 :     foo <- try(FUN(x))
757 :     ans[i,] <- if (inherits(foo, "try-error")) NA else foo
758 :     }
759 :     ans
760 : maechler 1299 }
761 : bates 1297
762 : bates 1336
763 : maechler 1222 formatVC <- function(varc, digits = max(3, getOption("digits") - 2))
764 :     { ## "format()" the 'VarCorr' matrix of the random effects -- for show()ing
765 : bates 1587 sc <- unname(attr(varc, "sc"))
766 : maechler 1222 recorr <- lapply(varc, function(el) el@factors$correlation)
767 :     reStdDev <- c(lapply(recorr, slot, "sd"), list(Residual = sc))
768 :     reLens <- unlist(c(lapply(reStdDev, length)))
769 :     nr <- sum(reLens)
770 :     reMat <- array('', c(nr, 4),
771 :     list(rep.int('', nr),
772 :     c("Groups", "Name", "Variance", "Std.Dev.")))
773 :     reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
774 : bates 2074 reMat[,2] <- c(unlist(lapply(varc, colnames)), "")
775 : maechler 1222 reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
776 :     reMat[,4] <- format(unlist(reStdDev), digits = digits)
777 :     if (any(reLens > 1)) {
778 :     maxlen <- max(reLens)
779 :     corr <-
780 : bates 1811 do.call("rBind",
781 : maechler 1222 lapply(recorr,
782 :     function(x, maxlen) {
783 :     x <- as(x, "matrix")
784 :     cc <- format(round(x, 3), nsmall = 3)
785 :     cc[!lower.tri(cc)] <- ""
786 :     nr <- dim(cc)[1]
787 :     if (nr >= maxlen) return(cc)
788 :     cbind(cc, matrix("", nr, maxlen-nr))
789 :     }, maxlen))
790 :     colnames(corr) <- c("Corr", rep.int("", maxlen - 1))
791 : bates 1811 cbind(reMat, rBind(corr, rep.int("", ncol(corr))))
792 : maechler 1222 } else reMat
793 :     }
794 :    
795 : maechler 1530 ## We need to define an S4 print method, since using an S3 print
796 :     ## method fails as soon as you call print() explicitly, e.g. when
797 :     ## wanting to specify options.
798 : maechler 1165
799 : maechler 1467 ## This is modeled a bit after print.summary.lm :
800 : maechler 1530 printMer <- function(x, digits = max(3, getOption("digits") - 3),
801 : bates 1823 correlation = TRUE, symbolic.cor = FALSE,
802 : maechler 1530 signif.stars = getOption("show.signif.stars"), ...)
803 : maechler 1467 {
804 :     so <- summary(x)
805 : bates 1588 REML <- so@status["REML"]
806 : maechler 1467 llik <- so@logLik
807 :     dev <- so@deviance
808 :     devc <- so@devComp
809 :     glz <- so@isG
810 : maechler 1165
811 : maechler 1467 cat(so@methTitle, "\n")
812 :     if (!is.null(so@call$formula))
813 :     cat("Formula:", deparse(so@call$formula),"\n")
814 :     if (!is.null(so@call$data))
815 :     cat(" Data:", deparse(so@call$data), "\n")
816 :     if (!is.null(so@call$subset))
817 :     cat(" Subset:",
818 :     deparse(asOneSidedFormula(so@call$subset)[[2]]),"\n")
819 :     if (glz)
820 :     cat(" Family: ", so@family$family, "(",
821 :     so@family$link, " link)\n", sep = "")
822 :     print(so@AICtab, digits = digits)
823 : maechler 1222
824 : maechler 1467 cat("Random effects:\n")
825 :     print(so@REmat, quote = FALSE, digits = digits, ...)
826 : maechler 1222
827 : maechler 1467 ngrps <- so@ngrps
828 :     cat(sprintf("number of obs: %d, groups: ", devc[1]))
829 :     cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
830 :     cat("\n")
831 : bates 1587 if (so@devComp[8] < 0)
832 :     cat("\nEstimated scale (compare to ", abs(so@devComp[8]), ") ", so@sigma, "\n")
833 : maechler 1467 if (nrow(so@coefs) > 0) {
834 :     cat("\nFixed effects:\n")
835 :     printCoefmat(so@coefs, zap.ind = 3, #, tst.ind = 4
836 :     digits = digits, signif.stars = signif.stars)
837 :     if(correlation) {
838 :     rn <- rownames(so@coefs)
839 :     corF <- so@vcov@factors$correlation
840 :     if (!is.null(corF)) {
841 :     p <- ncol(corF)
842 :     if (p > 1) {
843 :     cat("\nCorrelation of Fixed Effects:\n")
844 :     if (is.logical(symbolic.cor) && symbolic.cor) {
845 :     print(symnum(as(corF, "matrix"), abbr.col = NULL))
846 :     }
847 :     else {
848 :     corF <- matrix(format(round(corF@x, 3), nsmall = 3),
849 :     nc = p)
850 :     dimnames(corF) <- list(abbreviate(rn, minlen=11),
851 :     abbreviate(rn, minlen=6))
852 :     corF[!lower.tri(corF)] <- ""
853 :     print(corF[-1, -p, drop=FALSE], quote = FALSE)
854 :     }
855 :     }
856 :     }
857 :     }
858 :     }
859 :     invisible(x)
860 :     }
861 : bates 316
862 : maechler 1530 setMethod("print", "mer", printMer)
863 :     setMethod("show", "mer", function(object) printMer(object))
864 : maechler 1467
865 :    
866 : bates 1150 setMethod("vcov", signature(object = "mer"),
867 : bates 1588 function(object, REML = object@status["REML"], ...) {
868 : bates 1587 rr <- as(object@devComp[8]^2 * tcrossprod(solve(object@RXX)), "dpoMatrix")
869 : maechler 1335 rr@factors$correlation <- as(rr, "corMatrix")
870 : maechler 1222 rr
871 :     })
872 : bates 1150
873 : maechler 1222
874 : bates 316 ## calculates degrees of freedom for fixed effects Wald tests
875 :     ## This is a placeholder. The answers are generally wrong. It will
876 :     ## be very tricky to decide what a 'right' answer should be with
877 :     ## crossed random effects.
878 :    
879 : bates 1150 setMethod("getFixDF", signature(object="mer"),
880 : maechler 1222 function(object, ...) {
881 :     devc <- object@devComp
882 :     rep(as.integer(devc[1]- devc[2]), devc[2])
883 :     })
884 : bates 316
885 : bates 755 setMethod("logLik", signature(object="mer"),
886 : bates 1588 function(object, REML = object@status["REML"], ...) {
887 : maechler 1222 val <- -deviance(object, REML = REML)/2
888 :     devc <- as.integer(object@devComp[1:2])
889 :     attr(val, "nall") <- attr(val, "nobs") <- devc[1]
890 :     attr(val, "df") <- abs(devc[2]) +
891 : maechler 1280 length(.Call(mer_coef, object, 0))
892 : maechler 1222 attr(val, "REML") <- REML
893 :     class(val) <- "logLik"
894 :     val
895 :     })
896 : bates 446
897 : bates 1150 setMethod("VarCorr", signature(x = "mer"),
898 : bates 1588 function(x, REML = x@status["REML"], ...)
899 : bates 1150 {
900 : bates 1587 sc2 <- x@devComp[8]^2
901 : maechler 1222 cnames <- x@cnames
902 :     ans <- x@Omega
903 :     for (i in seq(a = ans)) {
904 :     el <- as(sc2 * solve(ans[[i]]), "dpoMatrix")
905 :     el@Dimnames <- list(cnames[[i]], cnames[[i]])
906 : maechler 1335 el@factors$correlation <- as(el, "corMatrix")
907 : maechler 1222 ans[[i]] <- el
908 :     }
909 : bates 1587 attr(ans, "sc") <- sqrt(sc2)
910 : maechler 1222 ans
911 : bates 1150 })
912 : deepayan 721
913 : bates 1150 setMethod("anova", signature(object = "mer"),
914 : maechler 1222 function(object, ...)
915 : bates 446 {
916 : maechler 1222 mCall <- match.call(expand.dots = TRUE)
917 :     dots <- list(...)
918 :     modp <- if (length(dots))
919 :     sapply(dots, is, "mer") | sapply(dots, is, "lm") else logical(0)
920 :     if (any(modp)) { # multiple models - form table
921 :     opts <- dots[!modp]
922 :     mods <- c(list(object), dots[modp])
923 :     names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)],
924 :     as.character)
925 :     mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE),
926 :     attr, "df"))]
927 :     calls <- lapply(mods, slot, "call")
928 :     data <- lapply(calls, "[[", "data")
929 :     if (any(data != data[[1]]))
930 :     stop("all models must be fit to the same data object")
931 :     header <- paste("Data:", data[[1]])
932 :     subset <- lapply(calls, "[[", "subset")
933 :     if (any(subset != subset[[1]]))
934 :     stop("all models must use the same subset")
935 :     if (!is.null(subset[[1]]))
936 :     header <-
937 :     c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
938 :     llks <- lapply(mods, logLik, REML = FALSE)
939 :     Df <- sapply(llks, attr, "df")
940 :     llk <- unlist(llks)
941 :     chisq <- 2 * pmax(0, c(NA, diff(llk)))
942 :     dfChisq <- c(NA, diff(Df))
943 :     val <- data.frame(Df = Df,
944 :     AIC = sapply(llks, AIC),
945 :     BIC = sapply(llks, BIC),
946 :     logLik = llk,
947 :     "Chisq" = chisq,
948 :     "Chi Df" = dfChisq,
949 :     "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
950 :     check.names = FALSE)
951 :     class(val) <- c("anova", class(val))
952 :     attr(val, "heading") <-
953 :     c(header, "Models:",
954 :     paste(names(mods),
955 :     unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
956 :     sep = ": "))
957 :     return(val)
958 :     }
959 :     else { ## ------ single model ---------------------
960 :     foo <- object
961 :     ss <- foo@rXy^2
962 :     ssr <- exp(foo@devComp["logryy2"])
963 :     names(ss) <- object@cnames[[".fixed"]]
964 :     asgn <- attr(foo@X, "assign")
965 :     terms <- foo@terms
966 :     nmeffects <- attr(terms, "term.labels")
967 :     if ("(Intercept)" %in% names(ss))
968 :     nmeffects <- c("(Intercept)", nmeffects)
969 :     ss <- unlist(lapply(split(ss, asgn), sum))
970 :     df <- unlist(lapply(split(asgn, asgn), length))
971 : bates 1856 #dfr <- getFixDF(object)
972 : maechler 1222 ms <- ss/df
973 :     #f <- ms/(ssr/dfr)
974 :     #P <- pf(f, df, dfr, lower.tail = FALSE)
975 :     #table <- data.frame(df, ss, ms, dfr, f, P)
976 :     table <- data.frame(df, ss, ms)
977 :     dimnames(table) <-
978 :     list(nmeffects,
979 :     # c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
980 :     c("Df", "Sum Sq", "Mean Sq"))
981 :     if ("(Intercept)" %in% nmeffects)
982 :     table <- table[-match("(Intercept)", nmeffects), ]
983 :     attr(table, "heading") <- "Analysis of Variance Table"
984 :     class(table) <- c("anova", "data.frame")
985 :     table
986 :     }
987 : bates 316 })
988 : bates 446
989 : bates 1150 setMethod("confint", signature(object = "mer"),
990 : maechler 1222 function(object, parm, level = 0.95, ...)
991 :     .NotYetImplemented()
992 :     )
993 : bates 446
994 : bates 1150 setMethod("fitted", signature(object = "mer"),
995 : maechler 1222 function(object, ...)
996 : maechler 1280 .Call(mer_fitted, object)
997 : maechler 1222 )
998 : bates 446
999 : bates 1150 setMethod("formula", signature(x = "mer"),
1000 : maechler 1222 function(x, ...)
1001 :     x@call$formula
1002 :     )
1003 : bates 449
1004 : bates 1636 setMethod("residuals", signature(object = "glmer"),
1005 :     function(object, ...) .NotYetImplemented())
1006 : bates 689
1007 : bates 1636 setMethod("residuals", signature(object = "lmer"),
1008 :     function(object, ...) object@y - fitted(object))
1009 :    
1010 : bates 1232 ## FIXME: There should not be two identical methods like this but I'm not
1011 :     ## sure how to pass the ... argument to a method for another generic
1012 :     ## cleanly.
1013 : bates 1636 setMethod("resid", signature(object = "glmer"),
1014 :     function(object, ...) .NotYetImplemented())
1015 : bates 689
1016 : bates 1636 setMethod("resid", signature(object = "lmer"),
1017 :     function(object, ...) object@y - fitted(object))
1018 :    
1019 : bates 1150 setMethod("summary", signature(object = "mer"),
1020 : maechler 1222 function(object, ...) {
1021 : bates 689
1022 : maechler 1280 fcoef <- .Call(mer_fixef, object)
1023 : maechler 1222 vcov <- vcov(object)
1024 :     corF <- vcov@factors$correlation
1025 :     ## DF <- getFixDF(object)
1026 :     coefs <- cbind("Estimate" = fcoef, "Std. Error" = corF@sd) #, DF = DF)
1027 : bates 1588 REML <- object@status["REML"]
1028 : maechler 1222 llik <- logLik(object, REML)
1029 :     dev <- object@deviance
1030 :     devc <- object@devComp
1031 :    
1032 : bates 1588 glz <- is(object, "glmer")
1033 : maechler 1222 methTitle <-
1034 :     if (glz)
1035 :     paste("Generalized linear mixed model fit using",
1036 : bates 1588 switch(object@status["glmm"],
1037 :     "PQL", "Laplace", "AGQ"))
1038 : maechler 1222 else paste("Linear mixed-effects model fit by",
1039 :     if(REML) "REML" else "maximum likelihood")
1040 :    
1041 :     AICframe <- {
1042 :     if (glz)
1043 :     data.frame(AIC = AIC(llik), BIC = BIC(llik),
1044 :     logLik = c(llik),
1045 :     deviance = -2*llik,
1046 :     row.names = "")
1047 :     else
1048 :     data.frame(AIC = AIC(llik), BIC = BIC(llik),
1049 :     logLik = c(llik),
1050 :     MLdeviance = dev["ML"],
1051 :     REMLdeviance = dev["REML"],
1052 :     row.names = "")
1053 :     }
1054 :     REmat <- formatVC(VarCorr(object))
1055 : bates 1587 if (object@devComp[8] < 0) REmat <- REmat[-nrow(REmat), , drop = FALSE]
1056 : maechler 1222
1057 :     if (nrow(coefs) > 0) {
1058 : bates 1587 if (object@devComp[8] >= 0) {
1059 : maechler 1222 stat <- coefs[,1]/coefs[,2]
1060 :     ##pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
1061 :     coefs <- cbind(coefs, "t value" = stat) #, "Pr(>|t|)" = pval)
1062 :     } else {
1063 :     coefs <- coefs[, 1:2, drop = FALSE]
1064 :     stat <- coefs[,1]/coefs[,2]
1065 :     pval <- 2*pnorm(abs(stat), lower = FALSE)
1066 :     coefs <- cbind(coefs, "z value" = stat, "Pr(>|z|)" = pval)
1067 :     }
1068 :     } ## else : append columns to 0-row matrix ...
1069 :    
1070 : bates 1586 new(if(is(object, "glmer")) "summary.glmer" else
1071 :     {if(is(object, "lmer")) "summary.lmer" else "summary.mer"},
1072 : maechler 1222 object,
1073 :     isG = glz,
1074 :     methTitle = methTitle,
1075 :     logLik = llik,
1076 :     ngrps = sapply(object@flist, function(x) length(levels(x))),
1077 : maechler 1280 sigma = .Call(mer_sigma, object, REML),
1078 : maechler 1222 coefs = coefs,
1079 :     vcov = vcov,
1080 :     REmat = REmat,
1081 :     AICtab= AICframe
1082 :     )
1083 :     })## summary()
1084 :    
1085 : bates 1704 setMethod("print", "mer", printMer)
1086 :     setMethod("show", "mer", function(object) printMer(object))
1087 :    
1088 :    
1089 : maechler 1222 ## Methods for "summary.*" objects:
1090 :     setMethod("vcov", signature(object = "summary.mer"),
1091 :     function(object) object@vcov)
1092 :     setMethod("logLik", signature(object = "summary.mer"),
1093 :     function(object) object@logLik)
1094 :     setMethod("deviance", signature(object = "summary.mer"),
1095 :     function(object) object@deviance)
1096 :     setMethod("summary", signature(object = "summary.mer"), function(object) object)
1097 :    
1098 :    
1099 : bates 1150 setMethod("update", signature(object = "mer"),
1100 : maechler 1222 function(object, formula., ..., evaluate = TRUE)
1101 : bates 1180 {
1102 : maechler 1222 call <- object@call
1103 :     if (is.null(call))
1104 :     stop("need an object with call slot")
1105 :     extras <- match.call(expand.dots = FALSE)$...
1106 :     if (!missing(formula.))
1107 :     call$formula <- update.formula(formula(object), formula.)
1108 :     if (length(extras) > 0) {
1109 :     existing <- !is.na(match(names(extras), names(call)))
1110 :     for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
1111 :     if (any(!existing)) {
1112 :     call <- c(as.list(call), extras[!existing])
1113 :     call <- as.call(call)
1114 :     }
1115 :     }
1116 :     if (evaluate)
1117 :     eval(call, parent.frame())
1118 :     else call
1119 : bates 1180 })
1120 : bates 888
1121 : bates 1150 simss <- function(fm0, fma, nsim)
1122 : bates 879 {
1123 : bates 1150 ysim <- simulate(fm0, nsim)
1124 : bates 1166 cv <- list(gradient = FALSE, msMaxIter = 200:200,
1125 : maechler 1222 msVerbose = 0:0)
1126 : bates 1150 sapply(ysim, function(yy) {
1127 : maechler 1280 .Call(mer_update_y, fm0, yy)
1128 : maechler 1222 LMEoptimize(fm0) <- cv
1129 : maechler 1280 .Call(mer_update_y, fma, yy)
1130 : maechler 1222 LMEoptimize(fma) <- cv
1131 :     exp(c(H0 = fm0@devComp[["logryy2"]],
1132 :     Ha = fma@devComp[["logryy2"]]))
1133 : bates 1150 })
1134 : bates 879 }
1135 : bates 1166
1136 :     ## Some leftover code from the old AGQ method in lmer.
1137 :     if (FALSE) {
1138 :     ### FIXME: For nf == 1 change this to an AGQ evaluation. Needs
1139 :     ### AGQ for nc > 1 first.
1140 :     fxd <- PQLpars[fixInd]
1141 :     loglik <- logLik(mer)
1142 : maechler 1188
1143 : bates 1166 if (method %in% c("Laplace", "AGQ")) {
1144 : maechler 1222 nAGQ <- 1
1145 :     if (method == "AGQ") { # determine nAGQ at PQL estimates
1146 :     dev11 <- devAGQ(PQLpars, 11)
1147 :     ## FIXME: Should this be an absolute or a relative tolerance?
1148 :     devTol <- sqrt(.Machine$double.eps) * abs(dev11)
1149 :     for (nAGQ in c(9, 7, 5, 3, 1))
1150 :     if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break
1151 :     nAGQ <- nAGQ + 2
1152 :     if (gVerb)
1153 :     cat(paste("Using", nAGQ, "quadrature points per column\n"))
1154 :     }
1155 :     obj <- function(pars)
1156 : maechler 1280 .Call(glmer_devAGQ, pars, GSpt, nAGQ)
1157 : maechler 1222 optimRes <-
1158 :     nlminb(PQLpars, obj,
1159 :     lower = ifelse(const, 5e-10, -Inf),
1160 :     control = list(trace = getOption("verbose"),
1161 :     iter.max = cv$msMaxIter))
1162 :     optpars <- optimRes$par
1163 :     if (optimRes$convergence != 0)
1164 :     warning("nlminb failed to converge")
1165 :     deviance <- optimRes$objective
1166 :     if (gVerb)
1167 :     cat(paste("convergence message", optimRes$message, "\n"))
1168 :     fxd[] <- optpars[fixInd] ## preserve the names
1169 : maechler 1280 .Call(lmer_coefGets, mer, optpars[-fixInd], 2)
1170 : bates 1166 }
1171 : maechler 1188
1172 : maechler 1280 .Call(glmer_finalize, GSpt)
1173 : bates 1166 loglik[] <- -deviance/2
1174 : maechler 1343 }## end{leftover}
1175 : bates 1283
1176 :     setMethod("isNested", "mer",
1177 :     function(x, ...) !(x@L@type[1]),
1178 :     valueClass = "logical")
1179 :    
1180 :     setMethod("denomDF", "mer",
1181 :     function(x, ...)
1182 :     {
1183 :     mm <- x@X
1184 :     aa <- attr(mm, "assign")
1185 :     tt <- x@terms
1186 :     if (!isNested(x))
1187 :     return(list(coef = as.numeric(rep(NA, length(x@fixef))),
1188 :     terms = as.numeric(rep(NA,
1189 :     length(attr(tt, "order"))))))
1190 :     hasintercept <- attr(tt, "intercept") > 0
1191 :     ## check which variables vary within levels of grouping factors
1192 :     vars <- eval(attr(tt, "variables"), x@frame)
1193 :     fl <- x@flist
1194 :     vv <- matrix(0:0, nrow = length(vars), ncol = length(fl),
1195 :     dimnames = list(NULL, names(fl)))
1196 :     ## replace this loop by C code.
1197 :     for (i in 1:nrow(ans)) # check if variables vary within factors
1198 :     for (j in 1:ncol(ans))
1199 :     ans[i,j] <- all(tapply(vars[[i]], fl[[j]],
1200 :     function(x) length(unique(x)) == 1))
1201 :     ## which terms vary within levels of which grouping factors?
1202 :     tv <- crossprod(attr(tt, "factors"), !ans)
1203 :     ## maximum level at which the term is constant
1204 :     ml <- apply(tv, 1, function(rr) max(0, which(as.logical(rr))))
1205 :     ## unravel assignment applied to terms
1206 :     ll <- attr(tt, "term.labels")
1207 :     if (hasintercept)
1208 :     ll <- c("(Intercept)", ll)
1209 :     aaa <- factor(aa, labels = ll)
1210 :     asgn <- split(order(aa), aaa)
1211 :     nco <- lapply(asgn, length) # number of coefficients per term
1212 :     nlev <- lapply(fl, function(x) length(levels(x)))
1213 :     if (hasintercept) asgn$"(Intercept)" <- NULL
1214 :     list(ml = ml, nco = nco, nlev = nlev)
1215 :     })
1216 : bates 1336
1217 :     hatTrace <- function(x)
1218 :     {
1219 :     stopifnot(is(x, "mer"))
1220 :     .Call(mer_hat_trace2, x)
1221 :     }
1222 : bates 1607
1223 : bates 1745 ## Check that the 'start' argument matches the form of the ST
1224 :     ## slot in mer2. If so, install start as the ST slot.
1225 :     setST <- function(mer, start)
1226 :     {
1227 :     ST <- mer@ST
1228 :     if (!is.list(start) || length(start) != length(ST) ||
1229 :     !all.equal(names(ST), names(start)))
1230 :     stop(paste("start must be a list of length", length(ST),
1231 :     "with names\n", paste(names(ST), collapse = ',')))
1232 :     for (i in seq(along = start)) {
1233 :     if (class(start[[i]]) != class(ST[[i]]) ||
1234 :     !all.equal(dim(start[[i]]), dim(ST[[i]])))
1235 :     stop(paste("start[[", i, "]] must be of class '", class(ST[[i]]),
1236 :     "' and dimension ", paste(dim(ST[[i]]), collapse = ','),
1237 :     sep = ''))
1238 :     }
1239 :     mer@ST <- start
1240 :     mer
1241 :     }
1242 :    
1243 : bates 1689 lmer2 <- function(formula, data, family = gaussian,
1244 : bates 1800 method = c("REML", "ML", "Laplace", "AGQ"),
1245 :     control = list(), start = NULL,
1246 :     subset, weights, na.action, offset, contrasts = NULL,
1247 :     model = TRUE, ...)
1248 : bates 1689 {
1249 :     method <- match.arg(method)
1250 :     formula <- as.formula(formula)
1251 :     if (length(formula) < 3) stop("formula must be a two-sided formula")
1252 :     cv <- do.call("lmerControl", control)
1253 : bates 1823
1254 : bates 1689 ## Establish model frame and fixed-effects model matrix and terms
1255 :     mc <- match.call()
1256 : bates 1856 fr <- lmerFrames(mc, formula, contrasts)
1257 : bates 1809 storage.mode(fr$X) <- "double" # when ncol(X) == 0, X is logical
1258 : bates 1823
1259 : bates 1800 ## Fit a glm to the fixed-effects only.
1260 : bates 1797 if(is.character(family))
1261 :     family <- get(family, mode = "function", envir = parent.frame(2))
1262 :     if(is.function(family)) family <- family()
1263 : bates 1800 if(is.null(family$family)) stop("'family' not recognized")
1264 :     glmFit <- glm.fit(fr$X, fr$Y, weights = fr$weights,
1265 :     offset = fr$offset, family = family,
1266 :     intercept = attr(fr$mt, "intercept") > 0)
1267 :     storage.mode(glmFit$y) <- "double"
1268 : bates 1823
1269 : bates 1689 ## establish factor list and Ztl
1270 : bates 1800 FL <- lmerFactorList(formula, fr$mf, !is.null(family))
1271 : bates 1689 Ztl <- with(FL, .Call(Ztl_sparse, fl, Ztl))
1272 : bates 1823
1273 : bates 1800 mer <- .Call(lmer2_create, fr, FL, Ztl, glmFit,
1274 :     method, mc, model)
1275 : bates 1802 rm(fr, FL, Ztl, glmFit)
1276 : bates 1800 if (!is.null(start)) mer <- setST(mer, start)
1277 : maechler 2026 if (!is(mer, "glmer2")) { # linear mixed model
1278 : bates 1800 .Call(lmer2_optimize, mer, cv$msVerbose)
1279 :     .Call(lmer2_update_effects, mer)
1280 : bates 1689 }
1281 : bates 1698 mer
1282 : bates 1689 }
1283 : bates 1704
1284 :     ## FIXME: The fixef and vcov methods should incorporate the permutation.
1285 :     ## Extract the fixed effects
1286 : bates 1800 setMethod("fixef", signature(object = "lmer2"),
1287 : bates 1704 function(object, ...) {
1288 :     ans <- object@fixef
1289 :     names(ans) <- object@cnames$.fixed
1290 :     ans
1291 :     })
1292 :    
1293 :     ## Extract the conditional variance-covariance matrix of the fixed effects
1294 : bates 1800 setMethod("vcov", signature(object = "lmer2"),
1295 : bates 1704 function(object, REML = 0, ...) {
1296 : bates 1800 rr <- as(.Call(lmer2_sigma, object, REML)^2 *
1297 :     crossprod(.Call(lmer2_vcov, object)), "dpoMatrix")
1298 : bates 1704 nms <- object@cnames[[".fixed"]]
1299 :     dimnames(rr) <- list(nms, nms)
1300 :     rr@factors$correlation <- as(rr, "corMatrix")
1301 :     rr
1302 :     })
1303 :    
1304 :     ## This is modeled a bit after print.summary.lm :
1305 :     printMer2 <- function(x, digits = max(3, getOption("digits") - 3),
1306 : bates 1823 correlation = TRUE, symbolic.cor = FALSE,
1307 : bates 1800 signif.stars = getOption("show.signif.stars"), ...)
1308 : bates 1704 {
1309 :     so <- summary(x)
1310 : bates 1778 REML <- so@dims["isREML"]
1311 : bates 1704 llik <- so@logLik
1312 :     dev <- so@deviance
1313 :     dims <- x@dims
1314 :     glz <- so@isG
1315 : bates 1823
1316 : bates 1704 cat(so@methTitle, "\n")
1317 : bates 1800 if (!is.null(x@call$formula))
1318 :     cat("Formula:", deparse(x@call$formula),"\n")
1319 :     if (!is.null(x@call$data))
1320 :     cat(" Data:", deparse(x@call$data), "\n")
1321 :     if (!is.null(x@call$subset))
1322 :     cat(" Subset:",
1323 :     deparse(asOneSidedFormula(x@call$subset)[[2]]),"\n")
1324 : maechler 2026 if (is(x, "glmer2"))
1325 : bates 1704 cat(" Family: ", so@family$family, "(",
1326 :     so@family$link, " link)\n", sep = "")
1327 :     print(so@AICtab, digits = digits)
1328 :    
1329 :     cat("Random effects:\n")
1330 :     print(so@REmat, quote = FALSE, digits = digits, ...)
1331 :    
1332 :     ngrps <- so@ngrps
1333 :     cat(sprintf("Number of obs: %d, groups: ", dims["n"]))
1334 :     cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
1335 :     cat("\n")
1336 :     if (is.na(so@sigma))
1337 :     cat("\nEstimated scale (compare to 1):",
1338 :     sqrt(exp(so@deviance["lr2"])/so@dims["n"]), "\n")
1339 :     if (nrow(so@coefs) > 0) {
1340 :     cat("\nFixed effects:\n")
1341 :     printCoefmat(so@coefs, zap.ind = 3, #, tst.ind = 4
1342 :     digits = digits, signif.stars = signif.stars)
1343 :     if(correlation) {
1344 :     rn <- rownames(so@coefs)
1345 :     corF <- so@vcov@factors$correlation
1346 :     if (!is.null(corF)) {
1347 :     p <- ncol(corF)
1348 :     if (p > 1) {
1349 :     cat("\nCorrelation of Fixed Effects:\n")
1350 :     if (is.logical(symbolic.cor) && symbolic.cor) {
1351 :     print(symnum(as(corF, "matrix"), abbr.col = NULL))
1352 :     }
1353 :     else {
1354 : bates 1800 corf <- matrix(format(round(corF@x, 3), nsmall = 3),
1355 : bates 1704 nc = p)
1356 : bates 1800 dimnames(corf) <- list(abbreviate(rn, minlen=11),
1357 : bates 1704 abbreviate(rn, minlen=6))
1358 : bates 1800 corf[!lower.tri(corf)] <- ""
1359 :     print(corf[-1, -p, drop=FALSE], quote = FALSE)
1360 : bates 1704 }
1361 :     }
1362 :     }
1363 :     }
1364 :     }
1365 :     invisible(x)
1366 :     }
1367 :    
1368 : bates 1800 setMethod("summary", signature(object = "lmer2"),
1369 : bates 1704 function(object, ...) {
1370 :     fcoef <- fixef(object)
1371 :     vcov <- vcov(object)
1372 :     corF <- vcov@factors$correlation
1373 :     dims <- object@dims
1374 :     ## DF <- getFixDF(object)
1375 :     coefs <- cbind("Estimate" = fcoef, "Std. Error" = corF@sd) #, DF = DF)
1376 : bates 1778 REML <- object@dims["isREML"]
1377 : bates 1704 llik <- logLik(object, REML)
1378 :     dev <- object@deviance
1379 : maechler 1730
1380 : bates 1704 glz <- is(object, "glmer")
1381 :     methTitle <-
1382 :     if (glz)
1383 :     paste("Generalized linear mixed model fit using",
1384 :     switch(object@status["glmm"],
1385 :     "PQL", "Laplace", "AGQ"))
1386 :     else paste("Linear mixed-effects model fit by",
1387 :     if(REML) "REML" else "maximum likelihood")
1388 :    
1389 :     AICframe <- {
1390 :     if (glz)
1391 :     data.frame(AIC = AIC(llik), BIC = BIC(llik),
1392 :     logLik = c(llik),
1393 :     deviance = -2*llik,
1394 :     row.names = "")
1395 :     else
1396 :     data.frame(AIC = AIC(llik), BIC = BIC(llik),
1397 :     logLik = c(llik),
1398 :     MLdeviance = dev["ML"],
1399 :     REMLdeviance = dev["REML"],
1400 :     row.names = "")
1401 :     }
1402 :     varcor <- VarCorr(object)
1403 :     REmat <- formatVC(varcor)
1404 :     if (is.na(attr(varcor, "sc")))
1405 :     REmat <- REmat[-nrow(REmat), , drop = FALSE]
1406 :    
1407 :     if (nrow(coefs) > 0) {
1408 : bates 1800 if (dims["famType"] >= 0) {
1409 : bates 1704 coefs <- coefs[, 1:2, drop = FALSE]
1410 :     stat <- coefs[,1]/coefs[,2]
1411 :     pval <- 2*pnorm(abs(stat), lower = FALSE)
1412 :     coefs <- cbind(coefs, "z value" = stat, "Pr(>|z|)" = pval)
1413 :     } else {
1414 :     stat <- coefs[,1]/coefs[,2]
1415 :     ##pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
1416 :     coefs <- cbind(coefs, "t value" = stat) #, "Pr(>|t|)" = pval)
1417 :     }
1418 :     } ## else : append columns to 0-row matrix ...
1419 :    
1420 :     new(if(is(object, "glmer")) "summary.glmer" else
1421 : bates 1800 {if(is(object, "lmer")) "summary.lmer" else "summary.lmer2"},
1422 : bates 1704 object,
1423 :     isG = glz,
1424 :     methTitle = methTitle,
1425 :     logLik = llik,
1426 :     ngrps = sapply(object@flist, function(x) length(levels(x))),
1427 : bates 1800 sigma = .Call(lmer2_sigma, object, REML),
1428 : bates 1704 coefs = coefs,
1429 :     vcov = vcov,
1430 :     REmat = REmat,
1431 :     AICtab= AICframe
1432 :     )
1433 :     })## summary()
1434 :    
1435 :     ## Extract the log-likelihood or restricted log-likelihood
1436 : bates 1800 setMethod("logLik", signature(object="lmer2"),
1437 : bates 1704 function(object, REML = NULL, ...) {
1438 :     dims <- object@dims
1439 :     if (is.null(REML) || is.na(REML[1]))
1440 : bates 1778 REML <- object@dims["isREML"]
1441 : bates 1704 val <- -deviance(object, REML = REML)/2
1442 :     attr(val, "nall") <- attr(val, "nobs") <- dims["n"]
1443 :     attr(val, "df") <-
1444 : bates 1800 dims["p"] + length(.Call(lmer2_getPars, object))
1445 : bates 1704 attr(val, "REML") <- as.logical(REML)
1446 :     class(val) <- "logLik"
1447 :     val
1448 :     })
1449 :    
1450 :     ## Extract the deviance
1451 : bates 1800 setMethod("deviance", signature(object="lmer2"),
1452 : bates 1704 function(object, REML = NULL, ...) {
1453 :     if (missing(REML) || is.null(REML) || is.na(REML[1]))
1454 : bates 1778 REML <- object@dims["isREML"]
1455 : bates 1704 object@deviance[ifelse(REML, "REML", "ML")]
1456 :     })
1457 :    
1458 :     # Create the VarCorr object of variances and covariances
1459 : bates 1800 setMethod("VarCorr", signature(x = "lmer2"),
1460 : bates 1704 function(x, REML = NULL, ...)
1461 :     {
1462 : bates 1800 sc <- .Call(lmer2_sigma, x, REML)
1463 : bates 1704 cnames <- x@cnames
1464 : bates 1715 ans <- x@ST
1465 : bates 1709 for (i in seq(along = ans)) {
1466 : bates 1704 ai <- ans[[i]]
1467 :     dm <- dim(ai)
1468 :     if (dm[1] < 2) {
1469 : bates 1709 el <- (sc * ai)^2
1470 : bates 1704 } else {
1471 : bates 1709 dd <- diag(ai)
1472 :     diag(ai) <- rep(1, dm[1])
1473 : maechler 1730 el <- sc^2 * crossprod(dd * t(ai))
1474 : bates 1704 }
1475 :     el <- as(el, "dpoMatrix")
1476 :     el@Dimnames <- list(cnames[[i]], cnames[[i]])
1477 :     el@factors$correlation <- as(el, "corMatrix")
1478 :     ans[[i]] <- el
1479 :     }
1480 :     attr(ans, "sc") <- sc
1481 :     ans
1482 :     })
1483 :    
1484 : bates 1889 # Create the VarCorr object of variances and covariances
1485 :     setMethod("VarCorr", signature(x = "nlmer"),
1486 : bates 1980 function(x, ...)
1487 : bates 1889 {
1488 :     sc <- sqrt(sum(x@deviance[c("bqd", "Sdr")])/x@dims["n"])
1489 :     cnames <- x@cnames
1490 :     ans <- x@ST
1491 :     for (i in seq(along = ans)) {
1492 :     ai <- ans[[i]]
1493 :     dm <- dim(ai)
1494 :     if (dm[1] < 2) {
1495 :     el <- (sc * ai)^2
1496 :     } else {
1497 :     dd <- diag(ai)
1498 :     diag(ai) <- rep(1, dm[1])
1499 :     el <- sc^2 * crossprod(dd * t(ai))
1500 :     }
1501 :     el <- as(el, "dpoMatrix")
1502 :     el@Dimnames <- list(cnames[[i]], cnames[[i]])
1503 :     el@factors$correlation <- as(el, "corMatrix")
1504 :     ans[[i]] <- el
1505 :     }
1506 :     attr(ans, "sc") <- sc
1507 :     ans
1508 :     })
1509 :    
1510 : bates 1800 setMethod("print", "lmer2", printMer2)
1511 :     setMethod("show", "lmer2", function(object) printMer2(object))
1512 : bates 1704
1513 :     ## Methods for "summary.*" objects:
1514 : bates 1800 setMethod("vcov", signature(object = "summary.lmer2"),
1515 : bates 1704 function(object) object@vcov)
1516 : bates 1800 setMethod("logLik", signature(object = "summary.lmer2"),
1517 : bates 1704 function(object) object@logLik)
1518 : bates 1800 setMethod("deviance", signature(object = "summary.lmer2"),
1519 : bates 1704 function(object) object@deviance)
1520 : bates 1800 setMethod("summary", signature(object = "summary.lmer2"),
1521 : bates 1704 function(object) object)
1522 : bates 1800 setMethod("anova", signature(object = "lmer2"),
1523 : bates 1722 function(object, ...)
1524 :     {
1525 :     mCall <- match.call(expand.dots = TRUE)
1526 :     dots <- list(...)
1527 :     modp <- if (length(dots))
1528 : bates 1800 sapply(dots, is, "lmer2") | sapply(dots, is, "lm") else logical(0)
1529 : bates 1722 if (any(modp)) { # multiple models - form table
1530 :     opts <- dots[!modp]
1531 :     mods <- c(list(object), dots[modp])
1532 :     names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)],
1533 :     as.character)
1534 :     mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE),
1535 :     attr, "df"))]
1536 :     calls <- lapply(mods, slot, "call")
1537 :     data <- lapply(calls, "[[", "data")
1538 :     if (any(data != data[[1]]))
1539 :     stop("all models must be fit to the same data object")
1540 :     header <- paste("Data:", data[[1]])
1541 :     subset <- lapply(calls, "[[", "subset")
1542 :     if (any(subset != subset[[1]]))
1543 :     stop("all models must use the same subset")
1544 :     if (!is.null(subset[[1]]))
1545 :     header <-
1546 :     c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
1547 :     llks <- lapply(mods, logLik, REML = FALSE)
1548 :     Df <- sapply(llks, attr, "df")
1549 :     llk <- unlist(llks)
1550 :     chisq <- 2 * pmax(0, c(NA, diff(llk)))
1551 :     dfChisq <- c(NA, diff(Df))
1552 :     val <- data.frame(Df = Df,
1553 :     AIC = sapply(llks, AIC),
1554 :     BIC = sapply(llks, BIC),
1555 :     logLik = llk,
1556 :     "Chisq" = chisq,
1557 :     "Chi Df" = dfChisq,
1558 :     "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
1559 :     check.names = FALSE)
1560 :     class(val) <- c("anova", class(val))
1561 :     attr(val, "heading") <-
1562 :     c(header, "Models:",
1563 :     paste(names(mods),
1564 :     unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
1565 :     sep = ": "))
1566 :     return(val)
1567 :     }
1568 :     else { ## ------ single model ---------------------
1569 :     foo <- object
1570 :     ss <- foo@rXy^2
1571 :     ssr <- exp(foo@devComp["logryy2"])
1572 :     names(ss) <- object@cnames[[".fixed"]]
1573 :     asgn <- attr(foo@X, "assign")
1574 :     terms <- foo@terms
1575 :     nmeffects <- attr(terms, "term.labels")
1576 :     if ("(Intercept)" %in% names(ss))
1577 :     nmeffects <- c("(Intercept)", nmeffects)
1578 :     ss <- unlist(lapply(split(ss, asgn), sum))
1579 :     df <- unlist(lapply(split(asgn, asgn), length))
1580 :     #dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
1581 :     ms <- ss/df
1582 :     #f <- ms/(ssr/dfr)
1583 :     #P <- pf(f, df, dfr, lower.tail = FALSE)
1584 :     #table <- data.frame(df, ss, ms, dfr, f, P)
1585 :     table <- data.frame(df, ss, ms)
1586 :     dimnames(table) <-
1587 :     list(nmeffects,
1588 :     # c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
1589 :     c("Df", "Sum Sq", "Mean Sq"))
1590 :     if ("(Intercept)" %in% nmeffects)
1591 :     table <- table[-match("(Intercept)", nmeffects), ]
1592 :     attr(table, "heading") <- "Analysis of Variance Table"
1593 :     class(table) <- c("anova", "data.frame")
1594 :     table
1595 :     }
1596 :     })
1597 : bates 1726
1598 :     ## Temporary function to convert the ST representation of the
1599 :     ## relative variance-covariance matrix returned by lmer2 into the
1600 :     ## Omega representation required by lmer
1601 :     ST2Omega <- function(ST)
1602 :     {
1603 :     if (nrow(ST) == 1) return(as(1/ST^2, "dpoMatrix"))
1604 :     dd <- diag(ST)
1605 :     T <- as(ST, "dtrMatrix")
1606 :     T@diag <- "U"
1607 :     crossprod(solve(T)/dd)
1608 :     }
1609 :    
1610 : bates 1762 ## Extract the random effects
1611 : bates 1800 setMethod("ranef", signature(object = "lmer2"),
1612 : bates 1762 function(object, postVar = FALSE, ...) {
1613 :     ans <- new("ranef.lmer",
1614 : bates 1800 lapply(.Call(lmer2_ranef, object),
1615 : bates 1762 data.frame, check.names = FALSE))
1616 :     names(ans) <- names(object@flist)
1617 :     if (postVar) {
1618 : bates 1800 pV <- .Call(lmer2_postVar, object)
1619 : bates 1762 for (i in seq(along = ans))
1620 :     attr(ans[[i]], "postVar") <- pV[[i]]
1621 :     }
1622 :     ans
1623 :     })
1624 : bates 1784
1625 :     ## Generate a Markov chain Monte Carlo sample from the posterior distribution
1626 :     ## of the parameters in a linear mixed model
1627 :     setMethod("mcmcsamp", signature(object = "lmer2"),
1628 :     function(object, n = 1, verbose = FALSE, saveb = FALSE,
1629 :     trans = TRUE, deviance = FALSE, ...)
1630 :     {
1631 : bates 1800 ans <- t(.Call(lmer2_MCMCsamp, object, saveb, n,
1632 : bates 1784 trans, verbose, deviance))
1633 :     attr(ans, "mcpar") <- as.integer(c(1, n, 1))
1634 :     class(ans) <- "mcmc"
1635 :     mcmccompnames(ans, object, saveb, trans,
1636 :     glmer=FALSE, deviance=deviance)
1637 :     })
1638 : bates 1809
1639 :     setMethod("simulate", "lmer2",
1640 :     function(object, nsim = 1, seed = NULL, ...)
1641 :     {
1642 :     if(!is.null(seed)) set.seed(seed)
1643 :     if(!exists(".Random.seed", envir = .GlobalEnv))
1644 :     runif(1) # initialize the RNG if necessary
1645 :     RNGstate <- .Random.seed
1646 :     dims <- object@dims
1647 :     ans <- array(0, c(dims["n"], nsim))
1648 :     fe <- fixef(object)
1649 :     re <- ranef(object)
1650 :     nc <- sapply(re, ncol)
1651 :     nr <- sapply(re, nrow)
1652 :     sz <- nc * nr
1653 :     vc <- VarCorr(object)
1654 : bates 1823
1655 : bates 1809 cholL <- lapply(vc, chol)
1656 :     n <- object@dims["n"]
1657 :     for (i in seq_len(nsim))
1658 : bates 1923 ans[, i] <- crossprod(object@ZXyt,
1659 : bates 1809 c(unlist(lapply(seq_along(re), function(k)
1660 :     (t(cholL[[k]]) %*%
1661 :     matrix(rnorm(sz[k]),
1662 :     nc = nr[k]))@x)),
1663 :     fe, 0))@x
1664 :     ans + rnorm(prod(dim(ans)), sd = attr(vc, "sc"))
1665 :     })
1666 :    
1667 : bates 1856 ## Remove the intercept row from a transposed model matrix
1668 :     ## first checking that it is not the only column
1669 :     rmIntr <- function(mm)
1670 :     {
1671 :     if (is.na(irow <- match("(Intercept)", rownames(mm)))) return(mm)
1672 :     if (ncol(mm) < 2)
1673 :     stop("lhs of a random-effects term cannot be an intercept only")
1674 :     mm[-irow, , drop = FALSE]
1675 :     }
1676 :    
1677 :    
1678 :     ## Fit a nonlinear mixed-effects model
1679 :     nlmer <- function(formula, data,
1680 : bates 1889 control = list(), start = NULL, verbose = FALSE,
1681 : bates 1856 subset, weights, na.action, contrasts = NULL,
1682 : bates 1843 model = TRUE, ...)
1683 :     {
1684 :     mc <- match.call()
1685 :     formula <- as.formula(formula)
1686 : bates 1860 if (length(formula) < 3) stop("formula must be a 3-part formula")
1687 : bates 1856 nlform <- as.formula(formula[[2]])
1688 : bates 1860 if (length(nlform) < 3)
1689 : bates 1843 stop("formula must be a 3-part formula")
1690 : bates 1856 nlmod <- as.call(nlform[[3]])
1691 :    
1692 : bates 1843 cv <- do.call("lmerControl", control)
1693 : bates 1889 if (missing(verbose)) verbose <- cv$msVerbose
1694 : bates 1843 if (is.numeric(start)) start <- list(fixed = start)
1695 : bates 1856 s <- length(pnames <- names(start$fixed))
1696 :     stopifnot(length(start$fixed) > 0, s > 0,
1697 :     inherits(data, "data.frame"), nrow(data) > 1)
1698 :     if (any(pnames %in% names(data)))
1699 :     stop("parameter names must be distinct from names of the variables in data")
1700 :     anms <- all.vars(nlmod)
1701 :     if (!all(pnames %in% anms))
1702 :     stop("not all parameter names are used in the nonlinear model expression")
1703 : bates 1823
1704 : bates 1856 if (!length(vnms <- setdiff(anms, pnames)))
1705 :     stop("there are no variables used in the nonlinear model expression")
1706 :     ## create a frame in which to evaluate the factor list
1707 :     fr <- lmerFrames(mc,
1708 :     eval(substitute(foo ~ bar,
1709 :     list(foo = nlform[[2]],
1710 :     bar = subnms(formula[[3]], lapply(pnames, as.name))))),
1711 :     contrasts, vnms)
1712 :     mf <- fr$mf
1713 :     attr(mf, "terms") <- NULL
1714 :     env <- new.env()
1715 :     lapply(names(mf), function(nm) assign(nm, env = env, mf[[nm]]))
1716 :     n <- nrow(mf)
1717 :     lapply(pnames, function(nm) assign(nm, env = env, rep(start$fixed[[nm]], length.out = n)))
1718 : bates 1843
1719 : bates 1856 mf <- mf[rep(seq_len(nrow(data)), s), ]
1720 :     row.names(mf) <- seq_len(nrow(mf))
1721 :     ss <- rep.int(nrow(data), s)
1722 :     for (nm in pnames) mf[[nm]] <- rep.int(as.numeric(nm == pnames), ss)
1723 :     # factor list and model matrices
1724 :     FL <- lmerFactorList(substitute(foo ~ bar, list(foo = nlform[[2]], bar = formula[[3]])),
1725 :     mf, -1)
1726 :     FL$Ztl <- lapply(FL$Ztl, rmIntr) # Remove any intercept columns
1727 : bates 1843
1728 : bates 1856 Xt <- t(Matrix(as.matrix(mf[,pnames]), sparse = TRUE))
1729 :     xnms <- colnames(fr$X)
1730 :     if (!is.na(icol <- match("(Intercept)",xnms))) xnms <- xnms[-icol]
1731 : bates 1860 Xt@Dimnames[[2]] <- NULL
1732 : bates 1856 ### FIXME: The only times there would be additional columns in the
1733 :     ### fixed effects would be as interactions with parameter names and
1734 :     ### they must be constructed differently
1735 :     if (length(xnms) > 0)
1736 :     Xt <- rBind(Xt, t(Matrix(fr$X[rep.int(seq_len(n), s), xnms, drop = FALSE])))
1737 :    
1738 :     wts <- fr$weights
1739 :     if (is.null(wts)) wts <- numeric(0)
1740 :     Ztl1 <- lapply(with(FL, .Call(Ztl_sparse, fl, Ztl)), drop0)
1741 :     Gp <- unname(c(0L, cumsum(unlist(lapply(Ztl1, nrow)))))
1742 :     Zt <- do.call(rBind, Ztl1)
1743 : bates 1864 attr(fr$mf, "terms") <- NULL
1744 : bates 1889 fixef <- unlist(start$fixed)
1745 :     storage.mode(fixef) <- "double"
1746 : bates 1860 val <- .Call(nlmer_create, env, nlmod, fr$mf, pnames, call = mc,
1747 :     FL$fl, Xt, Zt, unname(fr$Y), wts,
1748 :     cnames = lapply(FL$Ztl, rownames), Gp = Gp,
1749 : bates 1889 fixef = fixef)
1750 :     .Call(nlmer_optimize, val, verbose)
1751 :     .Call(nlmer_update_ranef, val)
1752 : bates 1860 val
1753 : bates 1856 }
1754 : bates 1889
1755 :     setMethod("show", "nlmer", function(object)
1756 :     {
1757 :     dims <- object@dims
1758 :     cat("Nonlinear mixed model fit by Laplace\n")
1759 :     if (!is.null(object@call$formula))
1760 :     cat("Formula:", deparse(object@call$formula),"\n")
1761 :     if (!is.null(object@call$data))
1762 :     cat(" Data:", deparse(object@call$data), "\n")
1763 :     if (!is.null(object@call$subset))
1764 :     cat(" Subset:",
1765 :     deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
1766 :    
1767 :     cat("Random effects:\n")
1768 :     print(formatVC(VarCorr(object)), quote = FALSE,
1769 :     digits = max(3, getOption("digits") - 3))
1770 :    
1771 :     cat(sprintf("Number of obs: %d, groups: ", dims["n"]))
1772 :     ngrps <- sapply(object@flist, function(x) length(levels(x)))
1773 :     cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
1774 :     cat("\n")
1775 :     cat("\nFixed effects:\n")
1776 :     print(object@fixef)
1777 :     invisible(object)
1778 :     })
1779 :    

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge