SCM

SCM Repository

[matrix] Annotation of /pkg/R/lmer.R
ViewVC logotype

Annotation of /pkg/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 799 - (view) (download)

1 : bates 767 # Methods for lmer and for the objects that it produces
2 : bates 689
3 :     ## Some utilities
4 :    
5 : bates 775 ## Return the index into the packed lower triangle
6 : bates 689 Lind <- function(i,j) {
7 :     if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))
8 :     ((i - 1) * i)/2 + j
9 : bates 446 }
10 :    
11 : bates 775 ## Return the pairs of expressions separated by vertical bars
12 : bates 769 findbars <- function(term)
13 :     {
14 :     if (is.name(term) || is.numeric(term)) return(NULL)
15 :     if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
16 :     if (!is.call(term)) stop("term must be of class call")
17 :     if (term[[1]] == as.name('|')) return(term)
18 :     if (length(term) == 2) return(findbars(term[[2]]))
19 :     c(findbars(term[[2]]), findbars(term[[3]]))
20 :     }
21 :    
22 : bates 775 ## Return the formula omitting the pairs of expressions
23 :     ## that are separated by vertical bars
24 : bates 769 nobars <- function(term)
25 :     {
26 :     # FIXME: is the is.name in the condition redundant?
27 :     # A name won't satisfy the first condition.
28 :     if (!('|' %in% all.names(term)) || is.name(term)) return(term)
29 :     if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
30 :     if (length(term) == 2) {
31 :     nb <- nobars(term[[2]])
32 :     if (is.null(nb)) return(NULL)
33 :     term[[2]] <- nb
34 :     return(term)
35 :     }
36 :     nb2 <- nobars(term[[2]])
37 :     nb3 <- nobars(term[[3]])
38 :     if (is.null(nb2)) return(nb3)
39 :     if (is.null(nb3)) return(nb2)
40 :     term[[2]] <- nb2
41 :     term[[3]] <- nb3
42 :     term
43 :     }
44 :    
45 : bates 775 ## Substitute the '+' function for the '|' function
46 : bates 769 subbars <- function(term)
47 :     {
48 :     if (is.name(term) || is.numeric(term)) return(term)
49 :     if (length(term) == 2) {
50 :     term[[2]] <- subbars(term[[2]])
51 :     return(term)
52 :     }
53 :     stopifnot(length(term) == 3)
54 :     if (is.call(term) && term[[1]] == as.name('|')) term[[1]] <- as.name('+')
55 :     term[[2]] <- subbars(term[[2]])
56 :     term[[3]] <- subbars(term[[3]])
57 :     term
58 :     }
59 :    
60 : bates 775 ## Control parameters for lmer
61 :     lmerControl <-
62 : bates 799 function(maxIter = 200,
63 : bates 769 msMaxIter = 200,
64 : bates 435 tolerance = sqrt((.Machine$double.eps)),
65 : bates 752 niterEM = 15,
66 : bates 435 msTol = sqrt(.Machine$double.eps),
67 :     msVerbose = getOption("verbose"),
68 : bates 769 PQLmaxIt = 30,
69 : bates 435 EMverbose = getOption("verbose"),
70 :     analyticGradient = TRUE,
71 : bates 775 analyticHessian = FALSE)
72 : bates 435 {
73 : bates 775 list(maxIter = as.integer(maxIter),
74 :     msMaxIter = as.integer(msMaxIter),
75 :     tolerance = as.double(tolerance),
76 :     niterEM = as.integer(niterEM),
77 :     msTol = as.double(msTol),
78 :     msVerbose = as.logical(msVerbose),
79 :     PQLmaxIt = as.integer(PQLmaxIt),
80 :     EMverbose = as.logical(EMverbose),
81 :     analyticGradient = as.logical(analyticGradient),
82 :     analyticHessian = as.logical(analyticHessian))
83 : bates 435 }
84 :    
85 : bates 755 setMethod("lmer", signature(formula = "formula"),
86 : bates 689 function(formula, data, family,
87 :     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
88 : bates 435 control = list(),
89 :     subset, weights, na.action, offset,
90 :     model = TRUE, x = FALSE, y = FALSE, ...)
91 : bates 755 { ## match and check parameters
92 :     if (length(formula) < 3) stop("formula must be a two-sided formula")
93 :     cv <- do.call("lmerControl", control)
94 :     if (lmm <- missing(family)) { # linear mixed model
95 :     method <- match.arg(method)
96 :     if (method %in% c("PQL", "Laplace", "AGQ")) {
97 :     warning(paste('Argument method = "', method,
98 :     '" is not meaningful for a linear mixed model.\n',
99 :     'Using method = "REML".\n', sep = ''))
100 :     method <- "REML"
101 :     }
102 :     } else { # generalized linear mixed model
103 :     method <- if (missing(method)) "PQL" else match.arg(method)
104 :     if (method == "ML") method <- "PQL"
105 :     if (method == "REML")
106 :     warning(paste('Argument method = "REML" is not meaningful',
107 :     'for a generalized linear mixed model.',
108 :     '\nUsing method = "PQL".\n'))
109 : bates 704 }
110 :    
111 : bates 755 ## evaluate glm.fit, a generalized linear fit of fixed effects only
112 :     mf <- match.call()
113 :     m <- match(c("family", "data", "subset", "weights",
114 :     "na.action", "offset"), names(mf), 0)
115 :     mf <- mf[c(1, m)]
116 :     frame.form <- subbars(formula) # substitute `+' for `|'
117 :     fixed.form <- nobars(formula) # remove any terms with `|'
118 : bates 767 if (inherits(fixed.form, "name")) # RHS is empty - use a constant
119 : bates 755 fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
120 :     environment(fixed.form) <- environment(frame.form) <- environment(formula)
121 :     mf$formula <- fixed.form
122 :     mf$x <- mf$model <- mf$y <- TRUE
123 :     mf[[1]] <- as.name("glm")
124 :     glm.fit <- eval(mf, parent.frame())
125 : bates 767 x <- glm.fit$x
126 :     y <- as.double(glm.fit$y)
127 : bates 769 family <- glm.fit$family
128 : bates 449
129 : bates 755 ## evaluate a model frame for fixed and random effects
130 : bates 435 mf$formula <- frame.form
131 : bates 755 mf$x <- mf$model <- mf$y <- mf$family <- NULL
132 : bates 435 mf$drop.unused.levels <- TRUE
133 : bates 755 mf[[1]] <- as.name("model.frame")
134 : bates 435 frm <- eval(mf, parent.frame())
135 : bates 755
136 : bates 435 ## grouping factors and model matrices for random effects
137 :     bars <- findbars(formula[[3]])
138 :     random <-
139 :     lapply(bars,
140 :     function(x) list(model.matrix(eval(substitute(~term,
141 :     list(term=x[[2]]))),
142 :     frm),
143 : bates 452 eval(substitute(as.factor(fac)[,drop = TRUE],
144 : bates 435 list(fac = x[[3]])), frm)))
145 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
146 : bates 755
147 : bates 435 ## order factor list by decreasing number of levels
148 : bates 449 nlev <- sapply(random, function(x) length(levels(x[[2]])))
149 : bates 452 if (any(diff(nlev) > 0)) {
150 : bates 449 random <- random[rev(order(nlev))]
151 : bates 435 }
152 : bates 767
153 :     ## Create the model matrices and a mixed-effects representation (mer)
154 : bates 435 mmats <- c(lapply(random, "[[", 1),
155 : bates 755 .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
156 :     mer <- .Call("lmer_create", lapply(random, "[[", 2),
157 :     mmats, method, PACKAGE = "Matrix")
158 : bates 767 if (lmm) { ## linear mixed model
159 : bates 755 .Call("lmer_initial", mer, PACKAGE="Matrix")
160 :     .Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
161 :     LMEoptimize(mer) <- cv
162 :     fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")
163 : bates 767 return(new("lmer",
164 : bates 769 mer,
165 : bates 767 assign = attr(x, "assign"),
166 :     call = match.call(),
167 : bates 769 family = family, fitted = fits,
168 :     fixed = fixef(mer),
169 :     frame = if (model) frm else data.frame(),
170 :     logLik = logLik(mer),
171 : bates 755 residuals = unname(model.response(frm) - fits),
172 : bates 769 terms = glm.fit$terms))
173 : bates 755 }
174 :    
175 :     ## The rest of the function applies to generalized linear mixed models
176 :     gVerb <- getOption("verbose")
177 : bates 776 eta <- glm.fit$linear.predictors
178 : bates 767 wts <- glm.fit$prior.weights
179 : bates 774 wtssqr <- wts * wts
180 : bates 767 offset <- glm.fit$offset
181 :     if (is.null(offset)) offset <- numeric(length(eta))
182 : bates 776 mu <- numeric(length(eta))
183 : bates 767
184 : bates 774 dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
185 : bates 767 linkinv <- quote(family$linkinv(eta))
186 :     mu.eta <- quote(family$mu.eta(eta))
187 :     variance <- quote(family$variance(mu))
188 : bates 775 LMEopt <- get("LMEoptimize<-")
189 :     doLMEopt <- quote(LMEopt(x = mer, value = cv))
190 : bates 767
191 : bates 775 GSpt <- .Call("glmer_init", environment())
192 :     .Call("glmer_PQL", GSpt) # obtain PQL estimates
193 : bates 755
194 : bates 774 fixInd <- seq(ncol(x))
195 :     ## pars[fixInd] == beta, pars[-fixInd] == theta
196 :     PQLpars <- c(fixef(mer),
197 :     .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))
198 : bates 775 ## set flag to skip fixed-effects in subsequent calls
199 :     mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
200 : bates 777 ## indicator of constrained parameters
201 :     const <- c(rep(FALSE, length(fixInd)),
202 :     unlist(lapply(mer@nc[seq(along = random)],
203 :     function(k) 1:((k*(k+1))/2) <= k)
204 :     ))
205 : bates 779 devAGQ <- function(pars, n)
206 :     .Call("glmer_devAGQ", pars, GSpt, n, PACKAGE = "Matrix")
207 :     ## FIXME: Change this to start at 11 and decrease. The
208 :     ## evaluation at the PQL estimates is done once only.
209 :    
210 :     deviance <- devAGQ(PQLpars, 1) # Laplacian approximation
211 : bates 777 fxd <- PQLpars[fixInd]
212 : bates 779 loglik <- logLik(mer)
213 : bates 775
214 : bates 777 if (method %in% c("Laplace", "AGQ")) {
215 : bates 779 nAGQ <- 1
216 :     if (method == "AGQ") { # determine nAGQ at PQL estimates
217 :     dev11 <- devAGQ(PQLpars, 11)
218 : bates 799 ## FIXME: Should this be an absolute or a relative tolerance?
219 : bates 779 devTol <- sqrt(.Machine$double.eps) * abs(dev11)
220 : bates 799 for (nAGQ in c(9, 7, 5, 3, 1))
221 : bates 779 if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break
222 : bates 799 nAGQ <- nAGQ + 2
223 :     if (gVerb)
224 :     cat(paste("Using", nAGQ, "quadrature points per column\n"))
225 : bates 779 }
226 :     obj <- function(pars)
227 :     .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
228 : bates 777 if (exists("nlminb", mode = "function")) {
229 : bates 755 optimRes <-
230 : bates 779 nlminb(PQLpars, obj,
231 : bates 755 lower = ifelse(const, 5e-10, -Inf),
232 :     control = list(trace = getOption("verbose"),
233 :     iter.max = cv$msMaxIter))
234 :     optpars <- optimRes$par
235 :     if (optimRes$convergence != 0)
236 :     warning("nlminb failed to converge")
237 : bates 779 deviance <- optimRes$objective
238 : bates 755 } else {
239 :     optimRes <-
240 : bates 779 optim(PQLpars, obj, method = "L-BFGS-B",
241 : bates 755 lower = ifelse(const, 5e-10, -Inf),
242 :     control = list(trace = getOption("verbose"),
243 : bates 776 maxit = cv$msMaxIter))
244 : bates 755 optpars <- optimRes$par
245 :     if (optimRes$convergence != 0)
246 :     warning("optim failed to converge")
247 : bates 779 deviance <- optimRes$value
248 : bates 755 }
249 : bates 774 if (gVerb) {
250 : bates 772 cat(paste("convergence message", optimRes$message, "\n"))
251 : bates 777 }
252 :     fxd[] <- optpars[fixInd] ## preserve the names
253 : bates 755 }
254 :    
255 : bates 776 .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
256 : bates 779 loglik[] <- -deviance/2
257 : bates 769 new("lmer", mer, frame = frm, terms = glm.fit$terms,
258 : bates 777 assign = attr(glm.fit$x, "assign"),
259 :     call = match.call(), family = family,
260 :     logLik = loglik, fixed = fxd)
261 : bates 435 })
262 :    
263 : bates 755 setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
264 : bates 316 function(x, value)
265 :     {
266 :     if (value$msMaxIter < 1) return(x)
267 :     nc <- x@nc
268 : bates 755 constr <- unlist(lapply(nc[1:(length(nc) - 2)],
269 :     function(k) 1:((k*(k+1))/2) <= k))
270 : bates 752 fn <- function(pars)
271 : bates 755 deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
272 :     gr <- NULL
273 :     if (value$analyticGradient)
274 :     gr <-
275 :     function(pars) {
276 :     if (!isTRUE(all.equal(pars,
277 :     .Call("lmer_coef", x,
278 :     2, PACKAGE = "Matrix"))))
279 :     .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
280 :     .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
281 :     }
282 : bates 777 optimRes <-
283 :     if (exists("nlminb", mode = "function"))
284 :     nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
285 :     fn, gr,
286 :     lower = ifelse(constr, 5e-10, -Inf),
287 :     control = list(iter.max = value$msMaxIter,
288 :     trace = as.integer(value$msVerbose)))
289 :     else
290 :     optim(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
291 :     fn, gr, method = "L-BFGS-B",
292 :     lower = ifelse(constr, 5e-10, -Inf),
293 :     control = list(maxit = value$msMaxIter,
294 :     trace = as.integer(value$msVerbose)))
295 : bates 755 .Call("lmer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
296 : bates 316 if (optimRes$convergence != 0) {
297 : bates 777 warning(paste("optim or nlminb returned message",
298 :     optimRes$message,"\n"))
299 : bates 316 }
300 :     return(x)
301 :     })
302 :    
303 : bates 413 setMethod("ranef", signature(object = "lmer"),
304 : bates 689 function(object, accumulate = FALSE, ...) {
305 :     val <- new("lmer.ranef",
306 :     lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),
307 :     data.frame, check.names = FALSE),
308 :     varFac = object@bVar,
309 :     stdErr = .Call("lmer_sigma", object,
310 : bates 755 object@method == "REML", PACKAGE = "Matrix"))
311 : bates 689 if (!accumulate || length(val@varFac) == 1) return(val)
312 :     ## check for nested factors
313 :     L <- object@L
314 :     if (any(sapply(seq(a = val), function(i) length(L[[Lind(i,i)]]@i))))
315 :     error("Require nested grouping factors to accumulate random effects")
316 :     val
317 : bates 316 })
318 :    
319 : bates 755 setMethod("fixef", signature(object = "mer"),
320 : bates 774 function(object, ...)
321 :     .Call("lmer_fixef", object, PACKAGE = "Matrix"))
322 : bates 316
323 : bates 769 setMethod("fixef", signature(object = "lmer"),
324 :     function(object, ...) object@fixed)
325 : deepayan 721
326 : bates 413 setMethod("VarCorr", signature(x = "lmer"),
327 : bates 316 function(x, REML = TRUE, useScale = TRUE, ...) {
328 : bates 550 val <- .Call("lmer_variances", x, PACKAGE = "Matrix")
329 : bates 316 for (i in seq(along = val)) {
330 :     dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])
331 :     val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
332 :     }
333 :     new("VarCorr",
334 : bates 449 scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),
335 : bates 316 reSumry = val,
336 :     useScale = useScale)
337 :     })
338 :    
339 : bates 413 setMethod("gradient", signature(x = "lmer"),
340 : bates 755 function(x, unconst, ...)
341 :     .Call("lmer_gradient", x, unconst, PACKAGE = "Matrix"))
342 : bates 316
343 : bates 449 setMethod("summary", signature(object = "lmer"),
344 :     function(object, ...)
345 : bates 769 new("summary.lmer", object,
346 : bates 727 showCorrelation = TRUE,
347 : bates 769 useScale = !((object@family)$family %in% c("binomial", "poisson"))))
348 : bates 316
349 : bates 449 setMethod("show", signature(object = "lmer"),
350 :     function(object)
351 : bates 769 show(new("summary.lmer", object,
352 : bates 727 showCorrelation = FALSE,
353 : bates 769 useScale = !((object@family)$family %in% c("binomial", "poisson")))))
354 :    
355 : bates 449 setMethod("show", "summary.lmer",
356 : bates 316 function(object) {
357 : bates 727 fcoef <- object@fixed
358 : bates 449 useScale <- object@useScale
359 :     corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"),
360 : bates 316 "corrmatrix")
361 :     DF <- getFixDF(object)
362 :     coefs <- cbind(fcoef, corF@stdDev, DF)
363 :     nc <- object@nc
364 :     dimnames(coefs) <-
365 :     list(names(fcoef), c("Estimate", "Std. Error", "DF"))
366 : bates 449 digits <- max(3, getOption("digits") - 2)
367 : bates 755 REML <- object@method == "REML"
368 : bates 727 llik <- object@logLik
369 : bates 449 dev <- object@deviance
370 :    
371 :     rdig <- 5
372 : bates 727 if (glz <- !(object@method %in% c("REML", "ML"))) {
373 :     cat(paste("Generalized linear mixed model fit using",
374 :     object@method, "\n"))
375 :     } else {
376 :     cat("Linear mixed-effects model fit by ")
377 : bates 755 cat(if(REML) "REML\n" else "maximum likelihood\n")
378 : bates 727 }
379 : bates 449 if (!is.null(object@call$formula)) {
380 :     cat("Formula:", deparse(object@call$formula),"\n")
381 :     }
382 :     if (!is.null(object@call$data)) {
383 :     cat(" Data:", deparse(object@call$data), "\n")
384 :     }
385 :     if (!is.null(object@call$subset)) {
386 :     cat(" Subset:",
387 :     deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
388 :     }
389 : bates 727 if (glz) {
390 : bates 750 cat(" Family: ", object@family$family, "(",
391 :     object@family$link, " link)\n", sep = "")
392 : bates 727 print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
393 : bates 449 logLik = c(llik),
394 : bates 727 deviance = -2*llik,
395 :     row.names = ""))
396 :     } else {
397 :     print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
398 :     logLik = c(llik),
399 : bates 750 MLdeviance = dev["ML"],
400 : bates 449 REMLdeviance = dev["REML"],
401 :     row.names = ""))
402 : bates 727 }
403 : bates 449 cat("Random effects:\n")
404 : bates 777 show(VarCorr(object, useScale = useScale))
405 : bates 449 ngrps <- lapply(object@flist, function(x) length(levels(x)))
406 :     cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))
407 :     cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
408 :     cat("\n")
409 :     if (!useScale)
410 :     cat("\nEstimated scale (compare to 1) ",
411 : bates 755 .Call("lmer_sigma", object, FALSE, PACKAGE = "Matrix"),
412 : bates 449 "\n")
413 :     if (nrow(coefs) > 0) {
414 :     if (useScale) {
415 :     stat <- coefs[,1]/coefs[,2]
416 :     pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
417 :     nms <- colnames(coefs)
418 :     coefs <- cbind(coefs, stat, pval)
419 :     colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")
420 :     } else {
421 :     coefs <- coefs[, 1:2, drop = FALSE]
422 :     stat <- coefs[,1]/coefs[,2]
423 :     pval <- 2*pnorm(abs(stat), lower = FALSE)
424 :     nms <- colnames(coefs)
425 :     coefs <- cbind(coefs, stat, pval)
426 :     colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
427 :     }
428 :     cat("\nFixed effects:\n")
429 :     printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
430 :     if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
431 :     rn <- rownames(coefs)
432 :     dimnames(corF) <- list(
433 :     abbreviate(rn, minlen=11),
434 :     abbreviate(rn, minlen=6))
435 :     if (!is.null(corF)) {
436 :     p <- NCOL(corF)
437 :     if (p > 1) {
438 :     cat("\nCorrelation of Fixed Effects:\n")
439 :     corF <- format(round(corF, 3), nsmall = 3)
440 :     corF[!lower.tri(corF)] <- ""
441 :     print(corF[-1, -p, drop=FALSE], quote = FALSE)
442 :     }
443 :     }
444 :     }
445 :     }
446 :     invisible(object)
447 : bates 316 })
448 :    
449 : deepayan 707
450 : bates 689
451 :    
452 : bates 316 ## calculates degrees of freedom for fixed effects Wald tests
453 :     ## This is a placeholder. The answers are generally wrong. It will
454 :     ## be very tricky to decide what a 'right' answer should be with
455 :     ## crossed random effects.
456 :    
457 : bates 413 setMethod("getFixDF", signature(object="lmer"),
458 : bates 316 function(object, ...)
459 :     {
460 :     nc <- object@nc[-seq(along = object@Omega)]
461 : bates 777 p <- abs(nc[1]) - 1
462 : bates 316 n <- nc[2]
463 :     rep(n-p, p)
464 :     })
465 :    
466 : bates 755 setMethod("logLik", signature(object="mer"),
467 :     function(object, REML = object@method == "REML", ...) {
468 : bates 446 val <- -deviance(object, REML = REML)/2
469 :     nc <- object@nc[-seq(a = object@Omega)]
470 :     attr(val, "nall") <- attr(val, "nobs") <- nc[2]
471 : bates 782 attr(val, "df") <- abs(nc[1]) +
472 : bates 755 length(.Call("lmer_coef", object, 0, PACKAGE = "Matrix"))
473 : bates 446 attr(val, "REML") <- REML
474 :     class(val) <- "logLik"
475 :     val
476 :     })
477 :    
478 : bates 769 setMethod("logLik", signature(object="lmer"),
479 :     function(object, ...) object@logLik)
480 : deepayan 721
481 : bates 446 setMethod("anova", signature(object = "lmer"),
482 :     function(object, ...)
483 :     {
484 :     mCall <- match.call(expand.dots = TRUE)
485 :     dots <- list(...)
486 :     modp <- logical(0)
487 :     if (length(dots))
488 :     modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm")
489 :     if (any(modp)) { # multiple models - form table
490 :     opts <- dots[!modp]
491 :     mods <- c(list(object), dots[modp])
492 :     names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)
493 :     mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE), attr, "df"))]
494 :     calls <- lapply(mods, slot, "call")
495 :     data <- lapply(calls, "[[", "data")
496 :     if (any(data != data[[1]])) stop("all models must be fit to the same data object")
497 :     header <- paste("Data:", data[[1]])
498 :     subset <- lapply(calls, "[[", "subset")
499 :     if (any(subset != subset[[1]])) stop("all models must use the same subset")
500 :     if (!is.null(subset[[1]]))
501 :     header <-
502 :     c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
503 :     llks <- lapply(mods, logLik, REML = FALSE)
504 :     Df <- sapply(llks, attr, "df")
505 :     llk <- unlist(llks)
506 :     chisq <- 2 * pmax(0, c(NA, diff(llk)))
507 :     dfChisq <- c(NA, diff(Df))
508 :     val <- data.frame(Df = Df,
509 :     AIC = sapply(llks, AIC),
510 :     BIC = sapply(llks, BIC),
511 :     logLik = llk,
512 :     "Chisq" = chisq,
513 :     "Chi Df" = dfChisq,
514 :     "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
515 :     check.names = FALSE)
516 :     class(val) <- c("anova", class(val))
517 :     attr(val, "heading") <-
518 : bates 690 c(header, "Models:",
519 : bates 446 paste(names(mods),
520 :     unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
521 : bates 690 sep = ": "))
522 : bates 446 return(val)
523 :     } else {
524 : bates 571 foo <- object
525 :     foo@status["factored"] <- FALSE
526 :     .Call("lmer_factor", foo, PACKAGE="Matrix")
527 :     dfr <- getFixDF(foo)
528 :     rcol <- ncol(foo@RXX)
529 :     ss <- foo@RXX[ , rcol]^2
530 :     ssr <- ss[[rcol]]
531 :     ss <- ss[seq(along = dfr)]
532 :     names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
533 :     asgn <- foo@assign
534 :     terms <- foo@terms
535 :     nmeffects <- attr(terms, "term.labels")
536 :     if ("(Intercept)" %in% names(ss))
537 :     nmeffects <- c("(Intercept)", nmeffects)
538 :     ss <- unlist(lapply(split(ss, asgn), sum))
539 :     df <- unlist(lapply(split(asgn, asgn), length))
540 :     dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
541 :     ms <- ss/df
542 :     f <- ms/(ssr/dfr)
543 :     P <- pf(f, df, dfr, lower.tail = FALSE)
544 :     table <- data.frame(df, ss, ms, dfr, f, P)
545 :     dimnames(table) <-
546 :     list(nmeffects,
547 :     c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
548 :     if ("(Intercept)" %in% nmeffects) table <- table[-1,]
549 :     attr(table, "heading") <- "Analysis of Variance Table"
550 :     class(table) <- c("anova", "data.frame")
551 :     table
552 : bates 446 }
553 : bates 316 })
554 : bates 446
555 :     setMethod("update", signature(object = "lmer"),
556 :     function(object, formula., ..., evaluate = TRUE)
557 :     {
558 :     call <- object@call
559 :     if (is.null(call))
560 :     stop("need an object with call component")
561 :     extras <- match.call(expand.dots = FALSE)$...
562 :     if (!missing(formula.))
563 :     call$formula <- update.formula(formula(object), formula.)
564 :     if (length(extras) > 0) {
565 :     existing <- !is.na(match(names(extras), names(call)))
566 :     for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
567 :     if (any(!existing)) {
568 :     call <- c(as.list(call), extras[!existing])
569 :     call <- as.call(call)
570 :     }
571 :     }
572 :     if (evaluate)
573 :     eval(call, parent.frame())
574 :     else call
575 :     })
576 :    
577 :    
578 :     setMethod("confint", signature(object = "lmer"),
579 :     function (object, parm, level = 0.95, ...)
580 :     {
581 :     cf <- fixef(object)
582 :     pnames <- names(cf)
583 :     if (missing(parm))
584 :     parm <- seq(along = pnames)
585 :     else if (is.character(parm))
586 :     parm <- match(parm, pnames, nomatch = 0)
587 :     a <- (1 - level)/2
588 :     a <- c(a, 1 - a)
589 :     pct <- paste(round(100 * a, 1), "%")
590 :     ci <- array(NA, dim = c(length(parm), 2),
591 :     dimnames = list(pnames[parm], pct))
592 :     ses <- sqrt(diag(vcov(object)))[parm]
593 : bates 449 ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
594 : bates 446 ci
595 :     })
596 :    
597 : bates 755 setMethod("deviance", "mer",
598 : bates 449 function(object, REML = NULL, ...) {
599 :     .Call("lmer_factor", object, PACKAGE = "Matrix")
600 :     if (is.null(REML))
601 : bates 755 REML <- object@method == "REML"
602 : bates 449 object@deviance[[ifelse(REML, "REML", "ML")]]
603 :     })
604 : bates 446
605 : deepayan 721
606 : bates 769 setMethod("deviance", "lmer",
607 :     function(object, ...) -2 * c(object@logLik))
608 : deepayan 721
609 : bates 769
610 : bates 449 setMethod("chol", signature(x = "lmer"),
611 :     function(x, pivot = FALSE, LINPACK = pivot) {
612 :     x@status["factored"] <- FALSE # force a decomposition
613 :     .Call("lmer_factor", x, PACKAGE = "Matrix")
614 :     })
615 :    
616 :     setMethod("solve", signature(a = "lmer", b = "missing"),
617 :     function(a, b, ...)
618 : bates 562 .Call("lmer_invert", a, PACKAGE = "Matrix")
619 : bates 449 )
620 :    
621 :     setMethod("formula", "lmer", function(x, ...) x@call$formula)
622 :    
623 :     setMethod("vcov", signature(object = "lmer"),
624 : bates 755 function(object, REML = object@method == "REML", useScale = TRUE,...) {
625 : bates 449 sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
626 :     rr <- object@RXX
627 :     nms <- object@cnames[[".fixed"]]
628 :     dimnames(rr) <- list(nms, nms)
629 :     nr <- nrow(rr)
630 :     rr <- rr[-nr, -nr, drop = FALSE]
631 :     rr <- rr %*% t(rr)
632 :     if (useScale) {
633 :     rr = sc^2 * rr
634 :     }
635 :     rr
636 :     })
637 :    
638 : bates 550 ## Extract the L matrix
639 :     setAs("lmer", "dtTMatrix",
640 :     function(from)
641 :     {
642 :     ## force a refactorization if the factors have been inverted
643 :     if (from@status["inverted"]) from@status["factored"] <- FALSE
644 :     .Call("lmer_factor", from, PACKAGE = "Matrix")
645 :     L <- lapply(from@L, as, "dgTMatrix")
646 :     nf <- length(from@D)
647 :     Gp <- from@Gp
648 :     nL <- Gp[nf + 1]
649 : bates 562 Li <- integer(0)
650 :     Lj <- integer(0)
651 :     Lx <- double(0)
652 : bates 550 for (i in 1:nf) {
653 :     for (j in 1:i) {
654 :     Lij <- L[[Lind(i, j)]]
655 : bates 562 Li <- c(Li, Lij@i + Gp[i])
656 :     Lj <- c(Lj, Lij@j + Gp[j])
657 :     Lx <- c(Lx, Lij@x)
658 : bates 550 }
659 :     }
660 : bates 562 new("dtTMatrix", Dim = as.integer(c(nL, nL)), i = Li, j = Lj, x = Lx,
661 : bates 550 uplo = "L", diag = "U")
662 :     })
663 : bates 562
664 :     ## Extract the ZZX matrix
665 :     setAs("lmer", "dsTMatrix",
666 :     function(from)
667 :     {
668 :     .Call("lmer_inflate", from, PACKAGE = "Matrix")
669 :     ZZpO <- lapply(from@ZZpO, as, "dgTMatrix")
670 :     ZZ <- lapply(from@ZtZ, as, "dgTMatrix")
671 :     nf <- length(ZZpO)
672 :     Gp <- from@Gp
673 :     nZ <- Gp[nf + 1]
674 :     Zi <- integer(0)
675 :     Zj <- integer(0)
676 :     Zx <- double(0)
677 :     for (i in 1:nf) {
678 :     ZZpOi <- ZZpO[[i]]
679 :     Zi <- c(Zi, ZZpOi@i + Gp[i])
680 :     Zj <- c(Zj, ZZpOi@j + Gp[i])
681 :     Zx <- c(Zx, ZZpOi@x)
682 :     if (i > 1) {
683 :     for (j in 1:(i-1)) {
684 :     ZZij <- ZZ[[Lind(i, j)]]
685 :     ## off-diagonal blocks are transposed
686 :     Zi <- c(Zi, ZZij@j + Gp[j])
687 :     Zj <- c(Zj, ZZij@i + Gp[i])
688 :     Zx <- c(Zx, ZZij@x)
689 :     }
690 :     }
691 :     }
692 :     new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx,
693 :     uplo = "U")
694 :     })
695 : bates 689
696 :     setMethod("fitted", signature(object = "lmer"),
697 : bates 691 function(object, ...)
698 :     napredict(attr(object@frame, "na.action"), object@fitted))
699 : bates 689
700 :     setMethod("residuals", signature(object = "lmer"),
701 : bates 691 function(object, ...)
702 :     naresid(attr(object@frame, "na.action"), object@residuals))
703 : bates 689
704 :     setMethod("resid", signature(object = "lmer"),
705 :     function(object, ...) do.call("residuals", c(list(object), list(...))))
706 :    
707 :     setMethod("coef", signature(object = "lmer"),
708 :     function(object, ...)
709 :     {
710 : bates 769 fef <- data.frame(rbind(object@fixed), check.names = FALSE)
711 : bates 689 ref <- as(ranef(object), "list")
712 :     names(ref) <- names(object@flist)
713 :     val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
714 :     for (i in seq(a = val)) {
715 :     refi <- ref[[i]]
716 :     row.names(val[[i]]) <- row.names(refi)
717 :     if (!all(names(refi) %in% names(fef)))
718 :     stop("unable to align random and fixed effects")
719 :     val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
720 :     }
721 :     new("lmer.coef", val)
722 :     })
723 :    
724 :     setMethod("plot", signature(x = "lmer.coef"),
725 :     function(x, y, ...)
726 :     {
727 : bates 777 if (require("lattice", quietly = TRUE)) {
728 :     varying <- unique(do.call("c",
729 :     lapply(x, function(el)
730 :     names(el)[sapply(el,
731 :     function(col)
732 :     any(col != col[1]))])))
733 :     gf <- do.call("rbind", lapply(x, "[", j = varying))
734 :     gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
735 :     switch(min(length(varying), 3),
736 :     qqmath(eval(substitute(~ x | .grp,
737 :     list(x = as.name(varying[1])))), gf, ...),
738 :     xyplot(eval(substitute(y ~ x | .grp,
739 :     list(y = as.name(varying[1]),
740 :     x = as.name(varying[2])))), gf, ...),
741 :     splom(~ gf | .grp, ...))
742 :     }
743 : bates 689 })
744 :    
745 :     setMethod("plot", signature(x = "lmer.ranef"),
746 :     function(x, y, ...)
747 :     {
748 : bates 777 if (require("lattice", quietly = TRUE))
749 :     lapply(x, function(x) {
750 :     cn <- lapply(colnames(x), as.name)
751 :     switch(min(ncol(x), 3),
752 :     qqmath(eval(substitute(~ x,
753 :     list(x = cn[[1]]))),
754 :     x, ...),
755 :     xyplot(eval(substitute(y ~ x,
756 :     list(y = cn[[1]],
757 :     x = cn[[2]]))),
758 :     x, ...),
759 :     splom(~ x, ...))
760 :     })
761 : bates 689 })
762 :    
763 :     setMethod("with", signature(data = "lmer"),
764 : bates 690 function(data, expr, ...) {
765 : bates 691 dat <- eval(data@call$data)
766 :     if (!is.null(na.act <- attr(data@frame, "na.action")))
767 :     dat <- dat[-na.act, ]
768 :     lst <- c(list(. = data), data@flist, data@frame, dat)
769 :     eval(substitute(expr), lst[unique(names(lst))])
770 :     })
771 : bates 690
772 : bates 691 setMethod("terms", signature(x = "lmer"),
773 :     function(x, ...) x@terms)
774 : bates 767
775 :     setMethod("show", signature(object="VarCorr"),
776 :     function(object)
777 :     {
778 :     digits <- max(3, getOption("digits") - 2)
779 :     useScale <- length(object@useScale) > 0 && object@useScale[1]
780 :     sc <- ifelse(useScale, object@scale, 1.)
781 :     reStdDev <- c(lapply(object@reSumry,
782 :     function(x, sc)
783 :     sc*x@stdDev,
784 :     sc = sc), list(Residual = sc))
785 :     reLens <- unlist(c(lapply(reStdDev, length)))
786 :     reMat <- array('', c(sum(reLens), 4),
787 :     list(rep('', sum(reLens)),
788 :     c("Groups", "Name", "Variance", "Std.Dev.")))
789 :     reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
790 :     reMat[,2] <- c(unlist(lapply(reStdDev, names)), "")
791 :     reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
792 :     reMat[,4] <- format(unlist(reStdDev), digits = digits)
793 :     if (any(reLens > 1)) {
794 :     maxlen <- max(reLens)
795 :     corr <-
796 :     do.call("rbind",
797 :     lapply(object@reSumry,
798 :     function(x, maxlen) {
799 :     cc <- format(round(x, 3), nsmall = 3)
800 :     cc[!lower.tri(cc)] <- ""
801 :     nr <- dim(cc)[1]
802 :     if (nr >= maxlen) return(cc)
803 :     cbind(cc, matrix("", nr, maxlen-nr))
804 :     }, maxlen))
805 :     colnames(corr) <- c("Corr", rep("", maxlen - 1))
806 :     reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr))))
807 :     }
808 :     if (!useScale) reMat <- reMat[-nrow(reMat),]
809 :     print(reMat, quote = FALSE)
810 :     })
811 : bates 769
812 : bates 782 glmmMCMC <- function(obj, method = c("full"), nsamp = 1)
813 : bates 781 {
814 :     if (!inherits(obj, "lmer")) stop("obj must be of class `lmer'")
815 :     if (obj@family$family == "gaussian" && obj@family$link == "identity")
816 :     warn("glmmMCMC not indended for Gaussian family with identity link")
817 :     cv <- Matrix:::lmerControl()
818 :     family <- obj@family
819 :     frm <- obj@frame
820 :     fixed.form <- Matrix:::nobars(obj@call$formula)
821 :     if (inherits(fixed.form, "name")) # RHS is empty - use a constant
822 :     fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
823 :     glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE,
824 :     y = TRUE)
825 :     x <- glm.fit$x
826 :     y <- as.double(glm.fit$y)
827 :     bars <- Matrix:::findbars(obj@call$formula[[3]])
828 :     random <-
829 :     lapply(bars,
830 :     function(x) list(model.matrix(eval(substitute(~term,
831 :     list(term=x[[2]]))),
832 :     frm),
833 :     eval(substitute(as.factor(fac)[,drop = TRUE],
834 :     list(fac = x[[3]])), frm)))
835 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
836 :     if (any(names(random) != names(obj@flist)))
837 :     random <- random[names(obj@flist)]
838 :     mmats <- c(lapply(random, "[[", 1),
839 :     .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
840 :     mer <- as(obj, "mer")
841 :    
842 :     eta <- glm.fit$linear.predictors # perhaps later change this to obj@fitted?
843 :     wts <- glm.fit$prior.weights
844 :     wtssqr <- wts * wts
845 :     offset <- glm.fit$offset
846 :     if (is.null(offset)) offset <- numeric(length(eta))
847 :     off <- numeric(length(eta))
848 :     mu <- numeric(length(eta))
849 :    
850 :     dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
851 :     linkinv <- quote(family$linkinv(eta))
852 :     mu.eta <- quote(family$mu.eta(eta))
853 :     variance <- quote(family$variance(mu))
854 :     LMEopt <- getAnywhere("LMEoptimize<-")
855 :     doLMEopt <- quote(LMEopt(x = mer, value = cv))
856 :    
857 :     GSpt <- .Call("glmer_init", environment())
858 :     nf <- length(obj@flist)
859 :     fixed <- obj@fixed
860 :     varc <- .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")
861 : bates 784 b <- .Call("lmer_ranef", mer, PACKAGE = "Matrix")
862 :     ans <- list(fixed = matrix(0, nr = length(fixed), nc = nsamp),
863 :     varc = matrix(0, nr = length(varc), nc = nsamp))
864 : bates 781 for (i in 1:nsamp) {
865 : bates 782 ## conditional means and variances of fixed effects
866 : bates 784 fixed <- .Call("glmer_fixed_update", GSpt, b, fixed, PACKAGE = "Matrix")
867 :     ans$fixed[,i] <- fixed
868 : bates 782 ## sample from the conditional distribution of beta given b and y
869 :     ## conditional means and variances of random_effects
870 :     .Call("glmer_bhat", GSpt, fixed, varc, PACKAGE = "Matrix")
871 : bates 784 print(bhat <- .Call("lmer_ranef", mer, PACKAGE = "Matrix"))
872 : bates 782 ## sample from the conditional distribution of b given beta and y
873 :     ## sample from the conditional distribution of varc given b
874 : bates 784 ans$varc[,i] <- varc
875 : bates 781 }
876 : bates 784 ans
877 : bates 781 }
878 :    

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