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 691 - (view) (download)
Original Path: branches/trunk-lme4/R/lmer.R

1 : bates 689 ## Methods for lmer and for the objects that it produces
2 :    
3 :     ## Some utilities
4 :    
5 : bates 446 contr.SAS <- function(n, contrasts = TRUE)
6 : bates 689 ## Eliminate this function after R-2.1.0 is released
7 :     contr.treatment(n,
8 :     if (is.numeric(n) && length(n) == 1) n else length(n),
9 :     contrasts)
10 :    
11 :     Lind <- function(i,j) {
12 :     if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))
13 :     ((i - 1) * i)/2 + j
14 : bates 446 }
15 :    
16 : bates 689 Dhalf <- function(from) {
17 :     D <- from@D
18 :     nf <- length(D)
19 :     Gp <- from@Gp
20 :     res <- array(0, rep(Gp[nf+1],2))
21 :     for (i in 1:nf) {
22 :     DD <- D[[i]]
23 :     dd <- dim(DD)
24 :     for (k in 1:dd[3]) {
25 :     mm <- array(DD[ , , k], dd[1:2])
26 :     base <- Gp[i] + (k - 1)*dd[1]
27 :     res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)
28 :     }
29 :     }
30 :     res
31 :     }
32 :    
33 : bates 435 lmerControl <- # Control parameters for lmer
34 :     function(maxIter = 50,
35 :     msMaxIter = 50,
36 :     tolerance = sqrt((.Machine$double.eps)),
37 :     niterEM = 20,
38 :     msTol = sqrt(.Machine$double.eps),
39 :     msVerbose = getOption("verbose"),
40 :     PQLmaxIt = 20,
41 :     EMverbose = getOption("verbose"),
42 :     analyticGradient = TRUE,
43 :     analyticHessian=FALSE)
44 :     {
45 :     list(maxIter = maxIter,
46 :     msMaxIter = msMaxIter,
47 :     tolerance = tolerance,
48 :     niterEM = niterEM,
49 :     msTol = msTol,
50 :     msVerbose = msVerbose,
51 :     PQLmaxIt = PQLmaxIt,
52 :     EMverbose=EMverbose,
53 :     analyticHessian=analyticHessian,
54 :     analyticGradient=analyticGradient)
55 :     }
56 :    
57 : bates 689 setMethod("lmer", signature(formula = "formula", family = "missing"),
58 :     function(formula, data, family,
59 :     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
60 : bates 435 control = list(),
61 :     subset, weights, na.action, offset,
62 :     model = TRUE, x = FALSE, y = FALSE, ...)
63 :     {
64 :     # match and check parameters
65 : bates 449 REML <- match.arg(method) == "REML"
66 : bates 435 controlvals <- do.call("lmerControl", control)
67 : bates 449 controlvals$REML <- REML
68 : bates 435 if (length(formula) < 3) stop("formula must be a two-sided formula")
69 : bates 449
70 :     mf <- match.call() # create the model frame as frm
71 : bates 435 m <- match(c("data", "subset", "weights", "na.action", "offset"),
72 :     names(mf), 0)
73 :     mf <- mf[c(1, m)]
74 :     mf[[1]] <- as.name("model.frame")
75 :     frame.form <- subbars(formula)
76 :     environment(frame.form) <- environment(formula)
77 :     mf$formula <- frame.form
78 :     mf$drop.unused.levels <- TRUE
79 :     frm <- eval(mf, parent.frame())
80 : bates 449
81 : bates 435 ## grouping factors and model matrices for random effects
82 :     bars <- findbars(formula[[3]])
83 :     random <-
84 :     lapply(bars,
85 :     function(x) list(model.matrix(eval(substitute(~term,
86 :     list(term=x[[2]]))),
87 :     frm),
88 : bates 452 eval(substitute(as.factor(fac)[,drop = TRUE],
89 : bates 435 list(fac = x[[3]])), frm)))
90 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
91 : bates 449
92 : bates 435 ## order factor list by decreasing number of levels
93 : bates 449 nlev <- sapply(random, function(x) length(levels(x[[2]])))
94 : bates 452 if (any(diff(nlev) > 0)) {
95 : bates 449 random <- random[rev(order(nlev))]
96 : bates 435 }
97 : bates 452 fixed.form <- nobars(formula)
98 :     if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula
99 : bates 571 Xmat <- model.matrix(fixed.form, frm)
100 : bates 435 mmats <- c(lapply(random, "[[", 1),
101 : bates 571 .fixed = list(cbind(Xmat, .response = model.response(frm))))
102 : bates 689 obj <- .Call("lmer_create", lapply(random, "[[", 2),
103 :     mmats, PACKAGE = "Matrix")
104 : bates 691 slot(obj, "frame") <- frm
105 : bates 689 slot(obj, "terms") <- attr(model.frame(fixed.form, data), "terms")
106 :     slot(obj, "assign") <- attr(Xmat, "assign")
107 :     slot(obj, "call") <- match.call()
108 :     slot(obj, "REML") <- REML
109 : bates 571 rm(Xmat)
110 : bates 435 .Call("lmer_initial", obj, PACKAGE="Matrix")
111 :     .Call("lmer_ECMEsteps", obj,
112 :     controlvals$niterEM,
113 :     controlvals$REML,
114 :     controlvals$EMverbose,
115 :     PACKAGE = "Matrix")
116 :     LMEoptimize(obj) <- controlvals
117 : bates 689 slot(obj, "residuals") <-
118 :     unname(model.response(frm) -
119 :     (slot(obj, "fitted") <-
120 :     .Call("lmer_fitted", obj, mmats, TRUE, PACKAGE = "Matrix")))
121 : bates 435 obj
122 :     })
123 :    
124 : bates 413 setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),
125 : bates 316 function(x, value)
126 :     {
127 :     if (value$msMaxIter < 1) return(x)
128 :     st <- ccoef(x) # starting values
129 :     nc <- x@nc
130 :     nc <- nc[1:(length(nc) - 2)]
131 :     constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
132 :     fn <- function(pars) {
133 :     ccoef(x) <- pars
134 :     deviance(x, REML = value$REML)
135 :     }
136 : bates 380 gr <- if (value$analyticGradient)
137 :     function(pars) {
138 : bates 468 if (!identical(TRUE,all.equal(pars, ccoef(x)))) ccoef(x) <- pars
139 : bates 457 grad <- gradient(x, REML = value$REML, unconst = TRUE)
140 : bates 380 grad[constr] <- -grad[constr]/pars[constr]
141 :     grad
142 :     } else NULL
143 : bates 316 optimRes <- optim(st, fn, gr,
144 :     method = "L-BFGS-B",
145 :     lower = ifelse(constr, 1e-10, -Inf),
146 : bates 362 control = list(maxit = value$msMaxIter,
147 : bates 411 trace = as.integer(value$msVerbose)))
148 : bates 316 if (optimRes$convergence != 0) {
149 :     warning(paste("optim returned message",optimRes$message,"\n"))
150 :     }
151 : bates 411 ccoef(x) <- optimRes$par
152 : bates 316 return(x)
153 :     })
154 :    
155 : bates 413 setMethod("ranef", signature(object = "lmer"),
156 : bates 689 function(object, accumulate = FALSE, ...) {
157 :     val <- new("lmer.ranef",
158 :     lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),
159 :     data.frame, check.names = FALSE),
160 :     varFac = object@bVar,
161 :     stdErr = .Call("lmer_sigma", object,
162 :     object@REML, PACKAGE = "Matrix"))
163 :     if (!accumulate || length(val@varFac) == 1) return(val)
164 :     ## check for nested factors
165 :     L <- object@L
166 :     if (any(sapply(seq(a = val), function(i) length(L[[Lind(i,i)]]@i))))
167 :     error("Require nested grouping factors to accumulate random effects")
168 :     val
169 : bates 316 })
170 :    
171 : bates 413 setMethod("fixef", signature(object = "lmer"),
172 : bates 316 function(object, ...) {
173 : bates 550 val <- .Call("lmer_fixef", object, PACKAGE = "Matrix")
174 : bates 316 val[-length(val)]
175 :     })
176 :    
177 : bates 413 setMethod("VarCorr", signature(x = "lmer"),
178 : bates 316 function(x, REML = TRUE, useScale = TRUE, ...) {
179 : bates 550 val <- .Call("lmer_variances", x, PACKAGE = "Matrix")
180 : bates 316 for (i in seq(along = val)) {
181 :     dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])
182 :     val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
183 :     }
184 :     new("VarCorr",
185 : bates 449 scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),
186 : bates 316 reSumry = val,
187 :     useScale = useScale)
188 :     })
189 :    
190 : bates 413 setMethod("gradient", signature(x = "lmer"),
191 : bates 316 function(x, REML, unconst, ...)
192 : bates 449 .Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix"))
193 : bates 316
194 : bates 449 setMethod("summary", signature(object = "lmer"),
195 :     function(object, ...)
196 :     new("summary.lmer", object, useScale = TRUE, showCorrelation = TRUE))
197 : bates 316
198 : bates 449 setMethod("show", signature(object = "lmer"),
199 :     function(object)
200 : bates 550 show(new("summary.lmer", object, useScale = TRUE,
201 :     showCorrelation = FALSE))
202 : bates 449 )
203 :    
204 :     setMethod("show", "summary.lmer",
205 : bates 316 function(object) {
206 :     fcoef <- fixef(object)
207 : bates 449 useScale <- object@useScale
208 :     corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"),
209 : bates 316 "corrmatrix")
210 :     DF <- getFixDF(object)
211 :     coefs <- cbind(fcoef, corF@stdDev, DF)
212 :     nc <- object@nc
213 :     dimnames(coefs) <-
214 :     list(names(fcoef), c("Estimate", "Std. Error", "DF"))
215 : bates 449 digits <- max(3, getOption("digits") - 2)
216 :     REML <- length(object@REML) > 0 && object@REML[1]
217 :     llik <- logLik(object)
218 :     dev <- object@deviance
219 :    
220 :     rdig <- 5
221 :     cat("Linear mixed-effects model fit by ")
222 :     cat(ifelse(object@REML, "REML\n", "maximum likelihood\n") )
223 :     if (!is.null(object@call$formula)) {
224 :     cat("Formula:", deparse(object@call$formula),"\n")
225 :     }
226 :     if (!is.null(object@call$data)) {
227 :     cat(" Data:", deparse(object@call$data), "\n")
228 :     }
229 :     if (!is.null(object@call$subset)) {
230 :     cat(" Subset:",
231 :     deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
232 :     }
233 :     print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
234 :     logLik = c(llik),
235 :     MLdeviance = dev["ML"],
236 :     REMLdeviance = dev["REML"],
237 :     row.names = ""))
238 :     cat("Random effects:\n")
239 :     show(VarCorr(object))
240 :     ngrps <- lapply(object@flist, function(x) length(levels(x)))
241 :     cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))
242 :     cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
243 :     cat("\n")
244 :     if (!useScale)
245 :     cat("\nEstimated scale (compare to 1) ",
246 :     .Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"),
247 :     "\n")
248 :     if (nrow(coefs) > 0) {
249 :     if (useScale) {
250 :     stat <- coefs[,1]/coefs[,2]
251 :     pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
252 :     nms <- colnames(coefs)
253 :     coefs <- cbind(coefs, stat, pval)
254 :     colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")
255 :     } else {
256 :     coefs <- coefs[, 1:2, drop = FALSE]
257 :     stat <- coefs[,1]/coefs[,2]
258 :     pval <- 2*pnorm(abs(stat), lower = FALSE)
259 :     nms <- colnames(coefs)
260 :     coefs <- cbind(coefs, stat, pval)
261 :     colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
262 :     }
263 :     cat("\nFixed effects:\n")
264 :     printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
265 :     if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
266 :     rn <- rownames(coefs)
267 :     dimnames(corF) <- list(
268 :     abbreviate(rn, minlen=11),
269 :     abbreviate(rn, minlen=6))
270 :     if (!is.null(corF)) {
271 :     p <- NCOL(corF)
272 :     if (p > 1) {
273 :     cat("\nCorrelation of Fixed Effects:\n")
274 :     corF <- format(round(corF, 3), nsmall = 3)
275 :     corF[!lower.tri(corF)] <- ""
276 :     print(corF[-1, -p, drop=FALSE], quote = FALSE)
277 :     }
278 :     }
279 :     }
280 :     }
281 :     invisible(object)
282 : bates 316 })
283 :    
284 : bates 689 setMethod("lmer", signature(formula = "formula"),
285 :     function(formula, family, data,
286 :     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
287 :     control = list(),
288 :     subset, weights, na.action, offset,
289 :     model = TRUE, x = FALSE, y = FALSE, ...)
290 :     {
291 :     gVerb <- getOption("verbose")
292 :     # match and check parameters
293 :     controlvals <- do.call("lmerControl", control)
294 :     controlvals$REML <- FALSE
295 :     if (length(formula) < 3) stop("formula must be a two-sided formula")
296 :    
297 :     ## initial glm fit
298 :     mf <- match.call()
299 :     m <- match(c("family", "data", "subset", "weights",
300 :     "na.action", "offset"),
301 :     names(mf), 0)
302 :     mf <- mf[c(1, m)]
303 :     mf[[1]] <- as.name("glm")
304 :     fixed.form <- nobars(formula)
305 :     if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula
306 :     environment(fixed.form) <- environment(formula)
307 :     mf$formula <- fixed.form
308 :     mf$x <- mf$model <- mf$y <- TRUE
309 :     glm.fit <- eval(mf, parent.frame())
310 :     family <- glm.fit$family
311 :     ## Note: offset is on the linear scale
312 :     offset <- glm.fit$offset
313 :     if (is.null(offset)) offset <- 0
314 :     weights <- sqrt(abs(glm.fit$prior.weights))
315 :     ## initial 'fitted' values on linear scale
316 :     etaold <- eta <- glm.fit$linear.predictors
317 :    
318 :     ## evaluation of model frame
319 :     mf$x <- mf$model <- mf$y <- mf$family <- NULL
320 :     mf$drop.unused.levels <- TRUE
321 :     this.form <- subbars(formula)
322 :     environment(this.form) <- environment(formula)
323 :     mf$formula <- this.form
324 :     mf[[1]] <- as.name("model.frame")
325 :     frm <- eval(mf, parent.frame())
326 :    
327 :     ## grouping factors and model matrices for random effects
328 :     bars <- findbars(formula[[3]])
329 :     random <-
330 :     lapply(bars,
331 :     function(x) list(model.matrix(eval(substitute(~term,
332 :     list(term=x[[2]]))),
333 :     frm),
334 :     eval(substitute(as.factor(fac)[,drop = TRUE],
335 :     list(fac = x[[3]])), frm)))
336 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
337 :    
338 :     ## order factor list by decreasing number of levels
339 :     nlev <- sapply(random, function(x) length(levels(x[[2]])))
340 :     if (any(diff(nlev) > 0)) {
341 :     random <- random[rev(order(nlev))]
342 :     }
343 :     mmats <- c(lapply(random, "[[", 1),
344 :     .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
345 :     ## FIXME: Use Xfrm and Xmat to get the terms and assign
346 :     ## slots, pass these to lmer_create, then destroy Xfrm, Xmat, etc.
347 : bates 691 obj <- .Call("lmer_create", lapply(random, "[[", 2),
348 :     mmats, PACKAGE = "Matrix")
349 :     slot(obj, "frame") <- frm
350 :     slot(obj, "terms") <- attr(glm.fit$model, "terms")
351 :     slot(obj, "assign") <- attr(glm.fit$x, "assign")
352 :     slot(obj, "call") <- match.call()
353 :     slot(obj, "REML") <- FALSE
354 : bates 689 rm(glm.fit)
355 :     .Call("lmer_initial", obj, PACKAGE="Matrix")
356 :     mmats.unadjusted <- mmats
357 :     mmats[[1]][1,1] <- mmats[[1]][1,1]
358 :     conv <- FALSE
359 :     firstIter <- TRUE
360 :     msMaxIter.orig <- controlvals$msMaxIter
361 :     responseIndex <- ncol(mmats$.fixed)
362 :    
363 :     for (iter in seq(length = controlvals$PQLmaxIt))
364 :     {
365 :     mu <- family$linkinv(eta)
366 :     dmu.deta <- family$mu.eta(eta)
367 :     ## weights (note: weights is already square-rooted)
368 :     w <- weights * dmu.deta / sqrt(family$variance(mu))
369 :     ## adjusted response (should be comparable to X \beta, not including offset
370 :     z <- eta - offset + (mmats.unadjusted$.fixed[, responseIndex] - mu) / dmu.deta
371 :     .Call("nlme_weight_matrix_list",
372 :     mmats.unadjusted, w, z, mmats, PACKAGE="Matrix")
373 :     .Call("lmer_update_mm", obj, mmats, PACKAGE="Matrix")
374 :     if (firstIter) {
375 :     .Call("lmer_initial", obj, PACKAGE="Matrix")
376 :     if (gVerb) cat(" PQL iterations convergence criterion\n")
377 :     }
378 :     .Call("lmer_ECMEsteps", obj,
379 :     controlvals$niterEM,
380 :     FALSE,
381 :     controlvals$EMverbose,
382 :     PACKAGE = "Matrix")
383 :     LMEoptimize(obj) <- controlvals
384 :     eta[] <- offset + ## FIXME: should the offset be here ?
385 :     .Call("lmer_fitted", obj,
386 :     mmats.unadjusted, TRUE, PACKAGE = "Matrix")
387 :     crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))
388 :     if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))
389 :     ## use this to determine convergence
390 :     if (crit < controlvals$tolerance) {
391 :     conv <- TRUE
392 :     break
393 :     }
394 :     etaold[] <- eta
395 :    
396 :     ## Changing number of iterations on second and
397 :     ## subsequent iterations.
398 :     if (firstIter)
399 :     {
400 :     controlvals$niterEM <- 2
401 :     controlvals$msMaxIter <- 10
402 :     firstIter <- FALSE
403 :     }
404 :     }
405 :     if (!conv) warning("IRLS iterations for glmm did not converge")
406 :     obj
407 :     })
408 :    
409 : bates 316 ## calculates degrees of freedom for fixed effects Wald tests
410 :     ## This is a placeholder. The answers are generally wrong. It will
411 :     ## be very tricky to decide what a 'right' answer should be with
412 :     ## crossed random effects.
413 :    
414 : bates 413 setMethod("getFixDF", signature(object="lmer"),
415 : bates 316 function(object, ...)
416 :     {
417 :     nc <- object@nc[-seq(along = object@Omega)]
418 :     p <- nc[1] - 1
419 :     n <- nc[2]
420 :     rep(n-p, p)
421 :     })
422 :    
423 : bates 446 setMethod("logLik", signature(object="lmer"),
424 :     function(object, REML = object@REML, ...) {
425 :     val <- -deviance(object, REML = REML)/2
426 :     nc <- object@nc[-seq(a = object@Omega)]
427 :     attr(val, "nall") <- attr(val, "nobs") <- nc[2]
428 : bates 689 attr(val, "df") <- nc[1] + length(ccoef(object))
429 : bates 446 attr(val, "REML") <- REML
430 :     class(val) <- "logLik"
431 :     val
432 :     })
433 :    
434 :     setMethod("anova", signature(object = "lmer"),
435 :     function(object, ...)
436 :     {
437 :     mCall <- match.call(expand.dots = TRUE)
438 :     dots <- list(...)
439 :     modp <- logical(0)
440 :     if (length(dots))
441 :     modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm")
442 :     if (any(modp)) { # multiple models - form table
443 :     opts <- dots[!modp]
444 :     mods <- c(list(object), dots[modp])
445 :     names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)
446 :     mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE), attr, "df"))]
447 :     calls <- lapply(mods, slot, "call")
448 :     data <- lapply(calls, "[[", "data")
449 :     if (any(data != data[[1]])) stop("all models must be fit to the same data object")
450 :     header <- paste("Data:", data[[1]])
451 :     subset <- lapply(calls, "[[", "subset")
452 :     if (any(subset != subset[[1]])) stop("all models must use the same subset")
453 :     if (!is.null(subset[[1]]))
454 :     header <-
455 :     c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
456 :     llks <- lapply(mods, logLik, REML = FALSE)
457 :     Df <- sapply(llks, attr, "df")
458 :     llk <- unlist(llks)
459 :     chisq <- 2 * pmax(0, c(NA, diff(llk)))
460 :     dfChisq <- c(NA, diff(Df))
461 :     val <- data.frame(Df = Df,
462 :     AIC = sapply(llks, AIC),
463 :     BIC = sapply(llks, BIC),
464 :     logLik = llk,
465 :     "Chisq" = chisq,
466 :     "Chi Df" = dfChisq,
467 :     "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
468 :     check.names = FALSE)
469 :     class(val) <- c("anova", class(val))
470 :     attr(val, "heading") <-
471 : bates 690 c(header, "Models:",
472 : bates 446 paste(names(mods),
473 :     unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
474 : bates 690 sep = ": "))
475 : bates 446 return(val)
476 :     } else {
477 : bates 571 foo <- object
478 :     foo@status["factored"] <- FALSE
479 :     .Call("lmer_factor", foo, PACKAGE="Matrix")
480 :     dfr <- getFixDF(foo)
481 :     rcol <- ncol(foo@RXX)
482 :     ss <- foo@RXX[ , rcol]^2
483 :     ssr <- ss[[rcol]]
484 :     ss <- ss[seq(along = dfr)]
485 :     names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
486 :     asgn <- foo@assign
487 :     terms <- foo@terms
488 :     nmeffects <- attr(terms, "term.labels")
489 :     if ("(Intercept)" %in% names(ss))
490 :     nmeffects <- c("(Intercept)", nmeffects)
491 :     ss <- unlist(lapply(split(ss, asgn), sum))
492 :     df <- unlist(lapply(split(asgn, asgn), length))
493 :     dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
494 :     ms <- ss/df
495 :     f <- ms/(ssr/dfr)
496 :     P <- pf(f, df, dfr, lower.tail = FALSE)
497 :     table <- data.frame(df, ss, ms, dfr, f, P)
498 :     dimnames(table) <-
499 :     list(nmeffects,
500 :     c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
501 :     if ("(Intercept)" %in% nmeffects) table <- table[-1,]
502 :     attr(table, "heading") <- "Analysis of Variance Table"
503 :     class(table) <- c("anova", "data.frame")
504 :     table
505 : bates 446 }
506 : bates 316 })
507 : bates 446
508 :     setMethod("update", signature(object = "lmer"),
509 :     function(object, formula., ..., evaluate = TRUE)
510 :     {
511 :     call <- object@call
512 :     if (is.null(call))
513 :     stop("need an object with call component")
514 :     extras <- match.call(expand.dots = FALSE)$...
515 :     if (!missing(formula.))
516 :     call$formula <- update.formula(formula(object), formula.)
517 :     if (length(extras) > 0) {
518 :     existing <- !is.na(match(names(extras), names(call)))
519 :     for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
520 :     if (any(!existing)) {
521 :     call <- c(as.list(call), extras[!existing])
522 :     call <- as.call(call)
523 :     }
524 :     }
525 :     if (evaluate)
526 :     eval(call, parent.frame())
527 :     else call
528 :     })
529 :    
530 :    
531 :     setMethod("confint", signature(object = "lmer"),
532 :     function (object, parm, level = 0.95, ...)
533 :     {
534 :     cf <- fixef(object)
535 :     pnames <- names(cf)
536 :     if (missing(parm))
537 :     parm <- seq(along = pnames)
538 :     else if (is.character(parm))
539 :     parm <- match(parm, pnames, nomatch = 0)
540 :     a <- (1 - level)/2
541 :     a <- c(a, 1 - a)
542 :     pct <- paste(round(100 * a, 1), "%")
543 :     ci <- array(NA, dim = c(length(parm), 2),
544 :     dimnames = list(pnames[parm], pct))
545 :     ses <- sqrt(diag(vcov(object)))[parm]
546 : bates 449 ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
547 : bates 446 ci
548 :     })
549 :    
550 : bates 689 setMethod("param", signature(object = "lmer"),
551 : bates 449 function(object, unconst = FALSE, ...) {
552 :     .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
553 :     })
554 : bates 446
555 : bates 449 setMethod("deviance", "lmer",
556 :     function(object, REML = NULL, ...) {
557 :     .Call("lmer_factor", object, PACKAGE = "Matrix")
558 :     if (is.null(REML))
559 :     REML <- if (length(oR <- object@REML)) oR else FALSE
560 :     object@deviance[[ifelse(REML, "REML", "ML")]]
561 :     })
562 : bates 446
563 : bates 449 setMethod("chol", signature(x = "lmer"),
564 :     function(x, pivot = FALSE, LINPACK = pivot) {
565 :     x@status["factored"] <- FALSE # force a decomposition
566 :     .Call("lmer_factor", x, PACKAGE = "Matrix")
567 :     })
568 :    
569 :     setMethod("solve", signature(a = "lmer", b = "missing"),
570 :     function(a, b, ...)
571 : bates 562 .Call("lmer_invert", a, PACKAGE = "Matrix")
572 : bates 449 )
573 :    
574 :     setMethod("formula", "lmer", function(x, ...) x@call$formula)
575 :    
576 :     setMethod("vcov", signature(object = "lmer"),
577 :     function(object, REML = object@REML, useScale = TRUE,...) {
578 :     sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
579 :     rr <- object@RXX
580 :     nms <- object@cnames[[".fixed"]]
581 :     dimnames(rr) <- list(nms, nms)
582 :     nr <- nrow(rr)
583 :     rr <- rr[-nr, -nr, drop = FALSE]
584 :     rr <- rr %*% t(rr)
585 :     if (useScale) {
586 :     rr = sc^2 * rr
587 :     }
588 :     rr
589 :     })
590 :    
591 : bates 550 ## Extract the L matrix
592 :     setAs("lmer", "dtTMatrix",
593 :     function(from)
594 :     {
595 :     ## force a refactorization if the factors have been inverted
596 :     if (from@status["inverted"]) from@status["factored"] <- FALSE
597 :     .Call("lmer_factor", from, PACKAGE = "Matrix")
598 :     L <- lapply(from@L, as, "dgTMatrix")
599 :     nf <- length(from@D)
600 :     Gp <- from@Gp
601 :     nL <- Gp[nf + 1]
602 : bates 562 Li <- integer(0)
603 :     Lj <- integer(0)
604 :     Lx <- double(0)
605 : bates 550 for (i in 1:nf) {
606 :     for (j in 1:i) {
607 :     Lij <- L[[Lind(i, j)]]
608 : bates 562 Li <- c(Li, Lij@i + Gp[i])
609 :     Lj <- c(Lj, Lij@j + Gp[j])
610 :     Lx <- c(Lx, Lij@x)
611 : bates 550 }
612 :     }
613 : bates 562 new("dtTMatrix", Dim = as.integer(c(nL, nL)), i = Li, j = Lj, x = Lx,
614 : bates 550 uplo = "L", diag = "U")
615 :     })
616 : bates 562
617 :     ## Extract the ZZX matrix
618 :     setAs("lmer", "dsTMatrix",
619 :     function(from)
620 :     {
621 :     .Call("lmer_inflate", from, PACKAGE = "Matrix")
622 :     ZZpO <- lapply(from@ZZpO, as, "dgTMatrix")
623 :     ZZ <- lapply(from@ZtZ, as, "dgTMatrix")
624 :     nf <- length(ZZpO)
625 :     Gp <- from@Gp
626 :     nZ <- Gp[nf + 1]
627 :     Zi <- integer(0)
628 :     Zj <- integer(0)
629 :     Zx <- double(0)
630 :     for (i in 1:nf) {
631 :     ZZpOi <- ZZpO[[i]]
632 :     Zi <- c(Zi, ZZpOi@i + Gp[i])
633 :     Zj <- c(Zj, ZZpOi@j + Gp[i])
634 :     Zx <- c(Zx, ZZpOi@x)
635 :     if (i > 1) {
636 :     for (j in 1:(i-1)) {
637 :     ZZij <- ZZ[[Lind(i, j)]]
638 :     ## off-diagonal blocks are transposed
639 :     Zi <- c(Zi, ZZij@j + Gp[j])
640 :     Zj <- c(Zj, ZZij@i + Gp[i])
641 :     Zx <- c(Zx, ZZij@x)
642 :     }
643 :     }
644 :     }
645 :     new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx,
646 :     uplo = "U")
647 :     })
648 : bates 689
649 :     setMethod("fitted", signature(object = "lmer"),
650 : bates 691 function(object, ...)
651 :     napredict(attr(object@frame, "na.action"), object@fitted))
652 : bates 689
653 :     setMethod("residuals", signature(object = "lmer"),
654 : bates 691 function(object, ...)
655 :     naresid(attr(object@frame, "na.action"), object@residuals))
656 : bates 689
657 :     setMethod("resid", signature(object = "lmer"),
658 :     function(object, ...) do.call("residuals", c(list(object), list(...))))
659 :    
660 :     setMethod("coef", signature(object = "lmer"),
661 :     function(object, ...)
662 :     {
663 :     fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
664 :     ref <- as(ranef(object), "list")
665 :     names(ref) <- names(object@flist)
666 :     val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
667 :     for (i in seq(a = val)) {
668 :     refi <- ref[[i]]
669 :     row.names(val[[i]]) <- row.names(refi)
670 :     if (!all(names(refi) %in% names(fef)))
671 :     stop("unable to align random and fixed effects")
672 :     val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
673 :     }
674 :     new("lmer.coef", val)
675 :     })
676 :    
677 :     setMethod("plot", signature(x = "lmer.coef"),
678 :     function(x, y, ...)
679 :     {
680 :     varying <- unique(do.call("c",
681 :     lapply(x, function(el)
682 :     names(el)[sapply(el,
683 :     function(col)
684 :     any(col != col[1]))])))
685 :     gf <- do.call("rbind", lapply(x, "[", j = varying))
686 :     gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
687 :     switch(min(length(varying), 3),
688 :     qqmath(eval(substitute(~ x | .grp,
689 :     list(x = as.name(varying[1])))), gf, ...),
690 :     xyplot(eval(substitute(y ~ x | .grp,
691 :     list(y = as.name(varying[1]),
692 :     x = as.name(varying[2])))), gf, ...),
693 :     splom(~ gf | .grp, ...))
694 :     })
695 :    
696 :     setMethod("plot", signature(x = "lmer.ranef"),
697 :     function(x, y, ...)
698 :     {
699 :     lapply(x, function(x) {
700 :     cn <- lapply(colnames(x), as.name)
701 :     switch(min(ncol(x), 3),
702 :     qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
703 :     xyplot(eval(substitute(y ~ x, list(y = cn[[1]], x = cn[[2]]))),
704 :     x, ...),
705 :     splom(~ x, ...))
706 :     })
707 :     })
708 :    
709 :     setMethod("with", signature(data = "lmer"),
710 : bates 690 function(data, expr, ...) {
711 : bates 691 dat <- eval(data@call$data)
712 :     if (!is.null(na.act <- attr(data@frame, "na.action")))
713 :     dat <- dat[-na.act, ]
714 :     lst <- c(list(. = data), data@flist, data@frame, dat)
715 :     eval(substitute(expr), lst[unique(names(lst))])
716 :     })
717 : bates 690
718 : bates 691 setMethod("terms", signature(x = "lmer"),
719 :     function(x, ...) x@terms)

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