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 751 - (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 :     Lind <- function(i,j) {
6 :     if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))
7 :     ((i - 1) * i)/2 + j
8 : bates 446 }
9 :    
10 : bates 689 Dhalf <- function(from) {
11 :     D <- from@D
12 :     nf <- length(D)
13 :     Gp <- from@Gp
14 :     res <- array(0, rep(Gp[nf+1],2))
15 :     for (i in 1:nf) {
16 :     DD <- D[[i]]
17 :     dd <- dim(DD)
18 :     for (k in 1:dd[3]) {
19 :     mm <- array(DD[ , , k], dd[1:2])
20 :     base <- Gp[i] + (k - 1)*dd[1]
21 :     res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)
22 :     }
23 :     }
24 :     res
25 :     }
26 :    
27 : bates 435 lmerControl <- # Control parameters for lmer
28 :     function(maxIter = 50,
29 :     msMaxIter = 50,
30 :     tolerance = sqrt((.Machine$double.eps)),
31 :     niterEM = 20,
32 :     msTol = sqrt(.Machine$double.eps),
33 :     msVerbose = getOption("verbose"),
34 :     PQLmaxIt = 20,
35 :     EMverbose = getOption("verbose"),
36 :     analyticGradient = TRUE,
37 :     analyticHessian=FALSE)
38 :     {
39 :     list(maxIter = maxIter,
40 :     msMaxIter = msMaxIter,
41 :     tolerance = tolerance,
42 :     niterEM = niterEM,
43 :     msTol = msTol,
44 :     msVerbose = msVerbose,
45 :     PQLmaxIt = PQLmaxIt,
46 :     EMverbose=EMverbose,
47 :     analyticHessian=analyticHessian,
48 :     analyticGradient=analyticGradient)
49 :     }
50 :    
51 : bates 689 setMethod("lmer", signature(formula = "formula", family = "missing"),
52 :     function(formula, data, family,
53 :     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
54 : bates 435 control = list(),
55 :     subset, weights, na.action, offset,
56 :     model = TRUE, x = FALSE, y = FALSE, ...)
57 :     {
58 :     # match and check parameters
59 : bates 704 method <- match.arg(method)
60 :     if (method %in% c("PQL", "Laplace", "AGQ")) {
61 :     warning(paste('Argument method = "', method,
62 :     '" is not meaningful for a linear mixed model.\n',
63 :     'Using method = "REML".\n', sep = ''))
64 :     method <- "REML"
65 :     }
66 : bates 435 controlvals <- do.call("lmerControl", control)
67 : bates 704 controlvals$REML <- REML <- method == "REML"
68 :    
69 : bates 435 if (length(formula) < 3) stop("formula must be a two-sided formula")
70 : bates 449
71 :     mf <- match.call() # create the model frame as frm
72 : bates 435 m <- match(c("data", "subset", "weights", "na.action", "offset"),
73 :     names(mf), 0)
74 :     mf <- mf[c(1, m)]
75 :     mf[[1]] <- as.name("model.frame")
76 :     frame.form <- subbars(formula)
77 :     environment(frame.form) <- environment(formula)
78 :     mf$formula <- frame.form
79 :     mf$drop.unused.levels <- TRUE
80 :     frm <- eval(mf, parent.frame())
81 : bates 449
82 : bates 435 ## grouping factors and model matrices for random effects
83 :     bars <- findbars(formula[[3]])
84 :     random <-
85 :     lapply(bars,
86 :     function(x) list(model.matrix(eval(substitute(~term,
87 :     list(term=x[[2]]))),
88 :     frm),
89 : bates 452 eval(substitute(as.factor(fac)[,drop = TRUE],
90 : bates 435 list(fac = x[[3]])), frm)))
91 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
92 : bates 449
93 : bates 435 ## order factor list by decreasing number of levels
94 : bates 449 nlev <- sapply(random, function(x) length(levels(x[[2]])))
95 : bates 452 if (any(diff(nlev) > 0)) {
96 : bates 449 random <- random[rev(order(nlev))]
97 : bates 435 }
98 : bates 452 fixed.form <- nobars(formula)
99 :     if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula
100 : bates 571 Xmat <- model.matrix(fixed.form, frm)
101 : bates 435 mmats <- c(lapply(random, "[[", 1),
102 : bates 571 .fixed = list(cbind(Xmat, .response = model.response(frm))))
103 : bates 689 obj <- .Call("lmer_create", lapply(random, "[[", 2),
104 :     mmats, PACKAGE = "Matrix")
105 : bates 691 slot(obj, "frame") <- frm
106 : bates 689 slot(obj, "terms") <- attr(model.frame(fixed.form, data), "terms")
107 :     slot(obj, "assign") <- attr(Xmat, "assign")
108 :     slot(obj, "call") <- match.call()
109 :     slot(obj, "REML") <- REML
110 : bates 571 rm(Xmat)
111 : bates 435 .Call("lmer_initial", obj, PACKAGE="Matrix")
112 :     .Call("lmer_ECMEsteps", obj,
113 :     controlvals$niterEM,
114 :     controlvals$REML,
115 :     controlvals$EMverbose,
116 :     PACKAGE = "Matrix")
117 :     LMEoptimize(obj) <- controlvals
118 : bates 689 slot(obj, "residuals") <-
119 :     unname(model.response(frm) -
120 :     (slot(obj, "fitted") <-
121 :     .Call("lmer_fitted", obj, mmats, TRUE, PACKAGE = "Matrix")))
122 : bates 435 obj
123 :     })
124 :    
125 : bates 413 setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),
126 : bates 316 function(x, value)
127 :     {
128 :     if (value$msMaxIter < 1) return(x)
129 :     st <- ccoef(x) # starting values
130 :     nc <- x@nc
131 :     nc <- nc[1:(length(nc) - 2)]
132 :     constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
133 :     fn <- function(pars) {
134 : bates 751 .Call("lmer_coefGets", x, pars, FALSE)
135 :     #ccoef(x) <- pars
136 : bates 316 deviance(x, REML = value$REML)
137 :     }
138 : bates 380 gr <- if (value$analyticGradient)
139 :     function(pars) {
140 : bates 751 if (!isTRUE(all.equal(pars, ccoef(x))))
141 :     .Call("lmer_coefGets", x, pars, FALSE)
142 :     #ccoef(x) <- pars
143 :     #grad <- gradient(x, REML = value$REML, unconst = TRUE)
144 :     gradient(x, REML = value$REML, unconst = FALSE)
145 :     #grad[constr] <- -grad[constr]/pars[constr]
146 :     #grad
147 : bates 380 } else NULL
148 : bates 751 if (require("port", quietly = TRUE)) {
149 :     optimRes <- portOptim(.Call("lmer_coef", x, FALSE),
150 :     fn, gr,
151 :     lower = ifelse(constr, 1e-10, -Inf),
152 :     control = list(maxit = value$msMaxIter,
153 :     trace = as.integer(value$msVerbose)))
154 :     .Call("lmer_coefGets", x, optimRes$par, FALSE)
155 :     } else {
156 :     optimRes <- optim(ccoef(x), fn, gr,
157 :     method = "L-BFGS-B",
158 :     lower = ifelse(constr, 1e-10, -Inf),
159 :     control = list(maxit = value$msMaxIter,
160 :     trace = as.integer(value$msVerbose)))
161 :     ccoef(x) <- optimRes$par
162 :     }
163 : bates 316 if (optimRes$convergence != 0) {
164 :     warning(paste("optim returned message",optimRes$message,"\n"))
165 :     }
166 :     return(x)
167 :     })
168 :    
169 : bates 413 setMethod("ranef", signature(object = "lmer"),
170 : bates 689 function(object, accumulate = FALSE, ...) {
171 :     val <- new("lmer.ranef",
172 :     lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),
173 :     data.frame, check.names = FALSE),
174 :     varFac = object@bVar,
175 :     stdErr = .Call("lmer_sigma", object,
176 :     object@REML, PACKAGE = "Matrix"))
177 :     if (!accumulate || length(val@varFac) == 1) return(val)
178 :     ## check for nested factors
179 :     L <- object@L
180 :     if (any(sapply(seq(a = val), function(i) length(L[[Lind(i,i)]]@i))))
181 :     error("Require nested grouping factors to accumulate random effects")
182 :     val
183 : bates 316 })
184 :    
185 : bates 413 setMethod("fixef", signature(object = "lmer"),
186 : bates 316 function(object, ...) {
187 : bates 550 val <- .Call("lmer_fixef", object, PACKAGE = "Matrix")
188 : bates 316 val[-length(val)]
189 :     })
190 :    
191 : deepayan 721 setMethod("fixef", signature(object = "glmer"),
192 :     function(object, ...) {
193 :     object@fixed
194 :     })
195 :    
196 : bates 413 setMethod("VarCorr", signature(x = "lmer"),
197 : bates 316 function(x, REML = TRUE, useScale = TRUE, ...) {
198 : bates 550 val <- .Call("lmer_variances", x, PACKAGE = "Matrix")
199 : bates 316 for (i in seq(along = val)) {
200 :     dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])
201 :     val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
202 :     }
203 :     new("VarCorr",
204 : bates 449 scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),
205 : bates 316 reSumry = val,
206 :     useScale = useScale)
207 :     })
208 :    
209 : bates 413 setMethod("gradient", signature(x = "lmer"),
210 : bates 316 function(x, REML, unconst, ...)
211 : bates 449 .Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix"))
212 : bates 316
213 : bates 449 setMethod("summary", signature(object = "lmer"),
214 :     function(object, ...)
215 : bates 727 new("summary.lmer", object, useScale = TRUE,
216 :     showCorrelation = TRUE,
217 :     method = if (object@REML) "REML" else "ML",
218 :     family = gaussian(),
219 :     logLik = logLik(object),
220 :     fixed = fixef(object)))
221 : bates 316
222 : deepayan 721 ## FIXME: glmm-s with scale not handled
223 :     setMethod("summary", signature(object = "glmer"),
224 :     function(object, ...)
225 : bates 727 new("summary.lmer", object, useScale = FALSE,
226 :     showCorrelation = TRUE,
227 :     method = object@method,
228 :     family = object@family,
229 :     logLik = logLik(object),
230 :     fixed = fixef(object)))
231 : deepayan 721
232 : bates 449 setMethod("show", signature(object = "lmer"),
233 :     function(object)
234 : bates 550 show(new("summary.lmer", object, useScale = TRUE,
235 : bates 727 showCorrelation = FALSE,
236 :     method = if (object@REML) "REML" else "ML",
237 :     family = gaussian(),
238 :     logLik = logLik(object),
239 :     fixed = fixef(object)))
240 : bates 449 )
241 : bates 727
242 : deepayan 721 setMethod("show", signature(object = "glmer"),
243 :     function(object)
244 :     show(new("summary.lmer", object, useScale = FALSE,
245 : bates 727 showCorrelation = FALSE,
246 :     method = object@method,
247 :     family = object@family,
248 :     logLik = logLik(object),
249 :     fixed = fixef(object)))
250 : deepayan 721 )
251 : bates 449
252 :     setMethod("show", "summary.lmer",
253 : bates 316 function(object) {
254 : bates 727 fcoef <- object@fixed
255 : bates 449 useScale <- object@useScale
256 :     corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"),
257 : bates 316 "corrmatrix")
258 :     DF <- getFixDF(object)
259 :     coefs <- cbind(fcoef, corF@stdDev, DF)
260 :     nc <- object@nc
261 :     dimnames(coefs) <-
262 :     list(names(fcoef), c("Estimate", "Std. Error", "DF"))
263 : bates 449 digits <- max(3, getOption("digits") - 2)
264 :     REML <- length(object@REML) > 0 && object@REML[1]
265 : bates 727 llik <- object@logLik
266 : bates 449 dev <- object@deviance
267 :    
268 :     rdig <- 5
269 : bates 727 if (glz <- !(object@method %in% c("REML", "ML"))) {
270 :     cat(paste("Generalized linear mixed model fit using",
271 :     object@method, "\n"))
272 :     } else {
273 :     cat("Linear mixed-effects model fit by ")
274 :     cat(if(object@REML) "REML\n" else "maximum likelihood\n")
275 :     }
276 : bates 449 if (!is.null(object@call$formula)) {
277 :     cat("Formula:", deparse(object@call$formula),"\n")
278 :     }
279 :     if (!is.null(object@call$data)) {
280 :     cat(" Data:", deparse(object@call$data), "\n")
281 :     }
282 :     if (!is.null(object@call$subset)) {
283 :     cat(" Subset:",
284 :     deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
285 :     }
286 : bates 727 if (glz) {
287 : bates 750 cat(" Family: ", object@family$family, "(",
288 :     object@family$link, " link)\n", sep = "")
289 : bates 727 print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
290 : bates 449 logLik = c(llik),
291 : bates 727 deviance = -2*llik,
292 :     row.names = ""))
293 :     } else {
294 :     print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
295 :     logLik = c(llik),
296 : bates 750 MLdeviance = dev["ML"],
297 : bates 449 REMLdeviance = dev["REML"],
298 :     row.names = ""))
299 : bates 727 }
300 : bates 449 cat("Random effects:\n")
301 :     show(VarCorr(object))
302 :     ngrps <- lapply(object@flist, function(x) length(levels(x)))
303 :     cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))
304 :     cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
305 :     cat("\n")
306 :     if (!useScale)
307 :     cat("\nEstimated scale (compare to 1) ",
308 :     .Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"),
309 :     "\n")
310 :     if (nrow(coefs) > 0) {
311 :     if (useScale) {
312 :     stat <- coefs[,1]/coefs[,2]
313 :     pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
314 :     nms <- colnames(coefs)
315 :     coefs <- cbind(coefs, stat, pval)
316 :     colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")
317 :     } else {
318 :     coefs <- coefs[, 1:2, drop = FALSE]
319 :     stat <- coefs[,1]/coefs[,2]
320 :     pval <- 2*pnorm(abs(stat), lower = FALSE)
321 :     nms <- colnames(coefs)
322 :     coefs <- cbind(coefs, stat, pval)
323 :     colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
324 :     }
325 :     cat("\nFixed effects:\n")
326 :     printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
327 :     if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
328 :     rn <- rownames(coefs)
329 :     dimnames(corF) <- list(
330 :     abbreviate(rn, minlen=11),
331 :     abbreviate(rn, minlen=6))
332 :     if (!is.null(corF)) {
333 :     p <- NCOL(corF)
334 :     if (p > 1) {
335 :     cat("\nCorrelation of Fixed Effects:\n")
336 :     corF <- format(round(corF, 3), nsmall = 3)
337 :     corF[!lower.tri(corF)] <- ""
338 :     print(corF[-1, -p, drop=FALSE], quote = FALSE)
339 :     }
340 :     }
341 :     }
342 :     }
343 :     invisible(object)
344 : bates 316 })
345 :    
346 : deepayan 707
347 : bates 689 setMethod("lmer", signature(formula = "formula"),
348 :     function(formula, family, data,
349 :     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
350 :     control = list(),
351 :     subset, weights, na.action, offset,
352 :     model = TRUE, x = FALSE, y = FALSE, ...)
353 :     {
354 : bates 704 if (missing(method)) {
355 :     method <- "PQL"
356 :     } else {
357 :     method <- match.arg(method)
358 :     if (method == "ML") method <- "PQL"
359 :     if (method == "REML")
360 :     warning(paste('Argument method = "REML" is not meaningful',
361 :     'for a generalized linear mixed model.',
362 :     '\nUsing method = "PQL".\n'))
363 :     }
364 : deepayan 707 ## if (method %in% c("Laplace", "AGQ"))
365 :     ## stop("Laplace and AGQ methods not yet implemented")
366 :     if (method %in% c("AGQ"))
367 :     stop("AGQ method not yet implemented")
368 : bates 689 gVerb <- getOption("verbose")
369 :     # match and check parameters
370 :     controlvals <- do.call("lmerControl", control)
371 :     controlvals$REML <- FALSE
372 :     if (length(formula) < 3) stop("formula must be a two-sided formula")
373 :     ## initial glm fit
374 :     mf <- match.call()
375 :     m <- match(c("family", "data", "subset", "weights",
376 :     "na.action", "offset"),
377 :     names(mf), 0)
378 :     mf <- mf[c(1, m)]
379 :     mf[[1]] <- as.name("glm")
380 :     fixed.form <- nobars(formula)
381 :     if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula
382 :     environment(fixed.form) <- environment(formula)
383 :     mf$formula <- fixed.form
384 :     mf$x <- mf$model <- mf$y <- TRUE
385 :     glm.fit <- eval(mf, parent.frame())
386 :     family <- glm.fit$family
387 :     ## Note: offset is on the linear scale
388 :     offset <- glm.fit$offset
389 :     if (is.null(offset)) offset <- 0
390 :     weights <- sqrt(abs(glm.fit$prior.weights))
391 :     ## initial 'fitted' values on linear scale
392 :     etaold <- eta <- glm.fit$linear.predictors
393 :    
394 :     ## evaluation of model frame
395 :     mf$x <- mf$model <- mf$y <- mf$family <- NULL
396 :     mf$drop.unused.levels <- TRUE
397 :     this.form <- subbars(formula)
398 :     environment(this.form) <- environment(formula)
399 :     mf$formula <- this.form
400 :     mf[[1]] <- as.name("model.frame")
401 :     frm <- eval(mf, parent.frame())
402 :    
403 :     ## grouping factors and model matrices for random effects
404 :     bars <- findbars(formula[[3]])
405 :     random <-
406 :     lapply(bars,
407 :     function(x) list(model.matrix(eval(substitute(~term,
408 :     list(term=x[[2]]))),
409 :     frm),
410 :     eval(substitute(as.factor(fac)[,drop = TRUE],
411 :     list(fac = x[[3]])), frm)))
412 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
413 :    
414 :     ## order factor list by decreasing number of levels
415 :     nlev <- sapply(random, function(x) length(levels(x[[2]])))
416 :     if (any(diff(nlev) > 0)) {
417 :     random <- random[rev(order(nlev))]
418 :     }
419 :     mmats <- c(lapply(random, "[[", 1),
420 :     .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
421 :     ## FIXME: Use Xfrm and Xmat to get the terms and assign
422 :     ## slots, pass these to lmer_create, then destroy Xfrm, Xmat, etc.
423 : bates 691 obj <- .Call("lmer_create", lapply(random, "[[", 2),
424 :     mmats, PACKAGE = "Matrix")
425 :     slot(obj, "frame") <- frm
426 :     slot(obj, "terms") <- attr(glm.fit$model, "terms")
427 :     slot(obj, "assign") <- attr(glm.fit$x, "assign")
428 :     slot(obj, "call") <- match.call()
429 :     slot(obj, "REML") <- FALSE
430 : bates 689 rm(glm.fit)
431 :     .Call("lmer_initial", obj, PACKAGE="Matrix")
432 :     mmats.unadjusted <- mmats
433 :     mmats[[1]][1,1] <- mmats[[1]][1,1]
434 :     conv <- FALSE
435 :     firstIter <- TRUE
436 :     msMaxIter.orig <- controlvals$msMaxIter
437 :     responseIndex <- ncol(mmats$.fixed)
438 :    
439 :     for (iter in seq(length = controlvals$PQLmaxIt))
440 :     {
441 :     mu <- family$linkinv(eta)
442 :     dmu.deta <- family$mu.eta(eta)
443 :     ## weights (note: weights is already square-rooted)
444 :     w <- weights * dmu.deta / sqrt(family$variance(mu))
445 :     ## adjusted response (should be comparable to X \beta, not including offset
446 :     z <- eta - offset + (mmats.unadjusted$.fixed[, responseIndex] - mu) / dmu.deta
447 :     .Call("nlme_weight_matrix_list",
448 :     mmats.unadjusted, w, z, mmats, PACKAGE="Matrix")
449 :     .Call("lmer_update_mm", obj, mmats, PACKAGE="Matrix")
450 :     if (firstIter) {
451 :     .Call("lmer_initial", obj, PACKAGE="Matrix")
452 :     if (gVerb) cat(" PQL iterations convergence criterion\n")
453 :     }
454 :     .Call("lmer_ECMEsteps", obj,
455 :     controlvals$niterEM,
456 :     FALSE,
457 :     controlvals$EMverbose,
458 :     PACKAGE = "Matrix")
459 :     LMEoptimize(obj) <- controlvals
460 :     eta[] <- offset + ## FIXME: should the offset be here ?
461 :     .Call("lmer_fitted", obj,
462 :     mmats.unadjusted, TRUE, PACKAGE = "Matrix")
463 :     crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))
464 :     if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))
465 :     ## use this to determine convergence
466 :     if (crit < controlvals$tolerance) {
467 :     conv <- TRUE
468 :     break
469 :     }
470 :     etaold[] <- eta
471 :    
472 :     ## Changing number of iterations on second and
473 :     ## subsequent iterations.
474 :     if (firstIter)
475 :     {
476 :     controlvals$niterEM <- 2
477 :     controlvals$msMaxIter <- 10
478 :     firstIter <- FALSE
479 :     }
480 :     }
481 :     if (!conv) warning("IRLS iterations for glmm did not converge")
482 : deepayan 706 controlvals$msMaxIter <- msMaxIter.orig
483 :    
484 :    
485 : deepayan 707 ### if (TRUE) ## Laplace
486 :     ### {
487 :     ## Need to optimize L(theta, beta) using Laplace approximation
488 : deepayan 706
489 : deepayan 707 ## Things needed for that:
490 :     ##
491 :     ## 1. reduced ssclme object, offset, weighted model matrices
492 :     ## 2. facs, reduced model matrices
493 : deepayan 706
494 : deepayan 707 ## Of these, those in 2 will be fixed given theta and beta,
495 :     ## and can be thought of arguments to the L(theta, beta)
496 :     ## function. However, the ones in 1 will have the same
497 :     ## structure. So the plan is to pre-allocate them and pass
498 :     ## them in too so they can be used without creating/copying
499 :     ## them more than once
500 : deepayan 706
501 :    
502 : deepayan 707 ## reduced ssclme
503 : deepayan 706
504 : deepayan 707 reducedObj <- .Call("lmer_collapse", obj, PACKAGE = "Matrix")
505 :     reducedMmats.unadjusted <- mmats.unadjusted
506 :     reducedMmats.unadjusted$.fixed <-
507 :     reducedMmats.unadjusted$.fixed[, responseIndex, drop = FALSE]
508 :     reducedMmats <- mmats
509 :     reducedMmats$.fixed <-
510 :     reducedMmats$.fixed[, responseIndex, drop = FALSE]
511 : deepayan 706
512 : deepayan 707 ## define function that calculates bhats given theta and beta
513 : deepayan 706
514 : deepayan 707 bhat <-
515 :     function(pars = NULL) # 1:(responseIndex-1) - beta, rest - theta
516 :     {
517 :     if (is.null(pars))
518 : deepayan 706 {
519 : deepayan 707 off <- drop(mmats.unadjusted$.fixed %*%
520 :     c(fixef(obj), 0)) + offset
521 :     }
522 :     else
523 :     {
524 :     .Call("lmer_coefGets",
525 :     reducedObj,
526 :     as.double(pars[responseIndex:length(pars)]),
527 :     TRUE,
528 :     PACKAGE = "Matrix")
529 :     off <- drop(mmats.unadjusted$.fixed %*%
530 :     c(pars[1:(responseIndex-1)], 0)) + offset
531 :     }
532 :     niter <- 20
533 :     conv <- FALSE
534 : deepayan 706
535 : deepayan 707 eta <- offset +
536 :     .Call("lmer_fitted",
537 :     obj, mmats.unadjusted, TRUE,
538 :     PACKAGE = "Matrix")
539 :     etaold <- eta
540 :    
541 :     for (iter in seq(length = niter))
542 :     {
543 :     mu <- family$linkinv(eta)
544 :     dmu.deta <- family$mu.eta(eta)
545 :     w <- weights * dmu.deta / sqrt(family$variance(mu))
546 :     z <- eta - off + (reducedMmats.unadjusted$.fixed[, 1]
547 :     - mu) / dmu.deta
548 :     .Call("nlme_weight_matrix_list",
549 :     reducedMmats.unadjusted, w, z, reducedMmats,
550 :     PACKAGE="Matrix")
551 :     .Call("lmer_update_mm",
552 :     reducedObj, reducedMmats,
553 :     PACKAGE="Matrix")
554 :     eta[] <- off +
555 :     .Call("lmer_fitted", reducedObj,
556 :     reducedMmats.unadjusted, TRUE,
557 : deepayan 706 PACKAGE = "Matrix")
558 : deepayan 707 ##cat(paste("bhat Criterion:", max(abs(eta - etaold)) /
559 :     ## (0.1 + max(abs(eta))), "\n"))
560 :     ## use this to determine convergence
561 :     if (max(abs(eta - etaold)) <
562 :     (0.1 + max(abs(eta))) * controlvals$tolerance)
563 : deepayan 706 {
564 : deepayan 707 conv <- TRUE
565 :     break
566 : deepayan 706 }
567 : deepayan 707 etaold[] <- eta
568 :    
569 :     }
570 :     if (!conv) warning("iterations for bhat did not converge")
571 : deepayan 706
572 : deepayan 707 ## bhat doesn't really need to return anything, we
573 :     ## just want the side-effect of modifying reducedObj
574 :     ## In particular, we are interested in
575 :     ## ranef(reducedObj) and reducedObj@bVar (?). But
576 :     ## the mu-scale response will be useful for log-lik
577 :     ## calculations later, so return them anyway
578 : deepayan 706
579 : deepayan 707 invisible(family$linkinv(eta))
580 :     }
581 : deepayan 706
582 : deepayan 707 ## function that calculates log likelihood (the thing that
583 :     ## needs to get evaluated at each Gauss-Hermite location)
584 :    
585 :     ## log scale ? worry about details later, get the pieces in place
586 :    
587 :     ## this is for the Laplace approximation only. GH is more
588 :     ## complicated
589 :    
590 :     devLaplace <- function(pars = NULL)
591 :     {
592 :     ## FIXME: This actually returns half the deviance.
593 : deepayan 706
594 : deepayan 707 ## gets correct values of bhat and bvars. As a side
595 :     ## effect, mu now has fitted values
596 :     mu <- bhat(pars = pars)
597 :    
598 :     ## GLM family log likelihood (upto constant ?)(log scale)
599 :     ## FIXME: need to adjust for sigma^2 for appropriate models (not trivial!)
600 :    
601 :     ## Keep everything on (log) likelihood scale
602 : deepayan 706
603 : deepayan 707 ## log lik from observations given fixed and random effects
604 :     ## get deviance, then multiply by -1/2 (since deviance = -2 log lik)
605 :     ans <- -sum(family$dev.resids(y = mmats.unadjusted$.fixed[, responseIndex],
606 :     mu = mu,
607 :     wt = weights^2))/2
608 :    
609 : deepayan 717 ## if (is.null(getOption("laplaceinR")))
610 :     ## {
611 :     ans <- ans +
612 :     .Call("lmer_laplace_devComp", reducedObj,
613 :     PACKAGE = "Matrix")
614 :     ## }
615 :     ## else
616 :     ## {
617 :     ## ranefs <- .Call("lmer_ranef", reducedObj, PACKAGE = "Matrix")
618 :     ## ## ans <- ans + reducedObj@devComp[2]/2 # log-determinant of Omega
619 : deepayan 706
620 : deepayan 717 ## Omega <- reducedObj@Omega
621 :     ## for (i in seq(along = ranefs))
622 :     ## {
623 :     ## ## contribution for random effects (get it working,
624 :     ## ## optimize later)
625 :     ## ## symmetrize RE variance
626 :     ## Omega[[i]] <- Omega[[i]] + t(Omega[[i]])
627 :     ## diag(Omega[[i]]) <- diag(Omega[[i]]) / 2
628 : deepayan 706
629 : deepayan 717 ## ## want log of `const det(Omega) exp(-1/2 b'
630 :     ## ## Omega b )` i.e., const + log det(Omega) - .5
631 :     ## ## * (b' Omega b)
632 : deepayan 706
633 : deepayan 717 ## ## FIXME: need to adjust for sigma^2 for appropriate
634 :     ## ## models (easy). These are all the b'Omega b,
635 :     ## ## summed as they eventually need to be. Think of
636 :     ## ## this as sum(rowSums((ranefs[[i]] %*% Omega[[i]])
637 :     ## ## * ranefs[[i]]))
638 : deepayan 706
639 : deepayan 717 ## ranef.loglik.det <- nrow(ranefs[[i]]) *
640 :     ## determinant(Omega[[i]], logarithm = TRUE)$modulus/2
641 : deepayan 711
642 : deepayan 717 ## # print(ranef.loglik.det)
643 : deepayan 711
644 : deepayan 717 ## ranef.loglik.re <-
645 :     ## -sum((ranefs[[i]] %*% Omega[[i]]) * ranefs[[i]])/2
646 : deepayan 711
647 : deepayan 717 ## # print(ranef.loglik.re)
648 : deepayan 711
649 : deepayan 717 ## ranef.loglik <- ranef.loglik.det + ranef.loglik.re
650 : deepayan 706
651 : deepayan 717 ## ## Jacobian adjustment
652 :     ## log.jacobian <-
653 :     ## sum(log(abs(apply(reducedObj@bVar[[i]],
654 :     ## 3,
655 : deepayan 706
656 : deepayan 717 ## ## next line depends on
657 :     ## ## whether bVars are variances
658 :     ## ## or Cholesly factors
659 : deepayan 706
660 : deepayan 717 ## ## function(x) sum(diag(x)))
661 :     ## ## Was this a bug? function(x) sum(diag( La.chol( x ) )))
662 :     ## function(x) prod(diag( La.chol( x ) )))
663 :     ## )))
664 : deepayan 706
665 : deepayan 717 ## # print(log.jacobian)
666 : deepayan 711
667 :    
668 : deepayan 717 ## ## the constant terms from the r.e. and the final
669 :     ## ## Laplacian integral cancel out both being:
670 :     ## ## ranef.loglik.constant <- 0.5 * length(ranefs[[i]]) * log(2 * base::pi)
671 : deepayan 710
672 : deepayan 717 ## ans <- ans + ranef.loglik + log.jacobian
673 :     ## }
674 :     ## }
675 : deepayan 707 ## ans is (up to some constant) log of the Laplacian
676 :     ## approximation of the likelihood. Return it's negative
677 :     ## to be minimized
678 : deepayan 706
679 : deepayan 707 ## cat("Parameters: ")
680 :     ## print(pars)
681 :    
682 :     ## cat("Value: ")
683 :     ## print(as.double(-ans))
684 :    
685 :     -ans
686 :     }
687 :    
688 :     if (method == "Laplace")
689 :     {
690 : bates 751 if (require("port")) {
691 :     optimRes <-
692 :     portOptim(fn = devLaplace,
693 :     par =
694 :     c(fixef(obj),
695 :     .Call("lmer_coef",
696 :     obj,
697 :     TRUE,
698 :     PACKAGE = "Matrix")),
699 : deepayan 707 ## WAS: coef(obj, unconst = TRUE)),
700 : bates 751 method = "BFGS", #hessian = TRUE,
701 : deepayan 707 control = list(trace = getOption("verbose"),
702 :     reltol = controlvals$msTol,
703 :     maxit = controlvals$msMaxIter))
704 : bates 751 optpars <- optimRes$par
705 :     } else {
706 :     optimRes <-
707 :     optim(fn =
708 :     c(fixef(obj), .Call("lmer_coef", obj, TRUE,
709 :     PACKAGE = "Matrix")),
710 :     devLaplace,
711 :     method = "BFGS", hessian = TRUE,
712 :     control = list(trace = getOption("verbose"),
713 :     reltol = controlvals$msTol,
714 :     maxit = controlvals$msMaxIter))
715 :     optpars <- optimRes$par
716 :     Hessian <- optimRes$hessian
717 :     }
718 : deepayan 707 if (optimRes$convergence != 0)
719 :     warning("optim failed to converge")
720 : bates 751
721 : deepayan 707 ##fixef(obj) <- optimRes$par[seq(length = responseIndex - 1)]
722 :     if (getOption("verbose")) {
723 :     cat(paste("optim convergence code",
724 :     optimRes$convergence, "\n"))
725 :     cat("Fixed effects:\n")
726 :     print(fixef(obj))
727 :     print(optimRes$par[seq(length = responseIndex - 1)])
728 :     cat("(Unconstrained) variance coefficients:\n")
729 : deepayan 706
730 : deepayan 707 print(
731 : deepayan 706 .Call("lmer_coef",
732 :     obj,
733 :     TRUE,
734 :     PACKAGE = "Matrix"))
735 : deepayan 707
736 :     ##coef(obj, unconst = TRUE) <-
737 :     ## optimRes$par[responseIndex:length(optimRes$par)]
738 :     ##print(coef(obj, unconst = TRUE))
739 :     print( optimRes$par[responseIndex:length(optimRes$par)] )
740 : deepayan 706 }
741 : deepayan 707
742 :     ## need to calculate likelihood. also need to store
743 :     ## new estimates of fixed effects somewhere
744 :     ## (probably cannot update standard errors)
745 : deepayan 717 ###Rprof(NULL)
746 : deepayan 707 }
747 :     else
748 :     {
749 :     optpars <-
750 :     c(fixef(obj),
751 :     .Call("lmer_coef",
752 :     obj,
753 :     TRUE,
754 :     PACKAGE = "Matrix"))
755 :     Hessian <- new("matrix")
756 :     }
757 : deepayan 706
758 :    
759 : deepayan 707 ## Before finishing, we need to call devLaplace with the
760 :     ## optimum pars to get the final log likelihood (still need
761 :     ## to make sure it's the actual likelihood and not a
762 :     ## multiple). This would automatically call bhat() and hence
763 :     ## have the 'correct' random effects in reducedObj.
764 : deepayan 706
765 : deepayan 721 loglik <- -devLaplace(optpars)
766 : deepayan 717 ##print(loglik)
767 : deepayan 707 ff <- optpars[1:(responseIndex-1)]
768 :     names(ff) <- names(fixef(obj))
769 : deepayan 706
770 : deepayan 707 if (!x) mmats <- list()
771 : deepayan 706
772 : deepayan 707 ### }
773 : deepayan 706
774 : deepayan 721 ##obj
775 :    
776 :     new("glmer", obj,
777 :     family = family,
778 :     glmmll = loglik,
779 :     method = method,
780 :     fixed = ff)
781 : bates 689 })
782 :    
783 : deepayan 707
784 : bates 316 ## calculates degrees of freedom for fixed effects Wald tests
785 :     ## This is a placeholder. The answers are generally wrong. It will
786 :     ## be very tricky to decide what a 'right' answer should be with
787 :     ## crossed random effects.
788 :    
789 : bates 413 setMethod("getFixDF", signature(object="lmer"),
790 : bates 316 function(object, ...)
791 :     {
792 :     nc <- object@nc[-seq(along = object@Omega)]
793 :     p <- nc[1] - 1
794 :     n <- nc[2]
795 :     rep(n-p, p)
796 :     })
797 :    
798 : bates 446 setMethod("logLik", signature(object="lmer"),
799 :     function(object, REML = object@REML, ...) {
800 :     val <- -deviance(object, REML = REML)/2
801 :     nc <- object@nc[-seq(a = object@Omega)]
802 :     attr(val, "nall") <- attr(val, "nobs") <- nc[2]
803 : bates 689 attr(val, "df") <- nc[1] + length(ccoef(object))
804 : bates 446 attr(val, "REML") <- REML
805 :     class(val) <- "logLik"
806 :     val
807 :     })
808 :    
809 : deepayan 721 setMethod("logLik", signature(object="glmer"),
810 :     function(object, ...) {
811 :     val <- object@glmmll
812 :     nc <- object@nc[-seq(a = object@Omega)]
813 :     attr(val, "nall") <- attr(val, "nobs") <- nc[2]
814 :     attr(val, "df") <- nc[1] + length(ccoef(object))
815 :     ## attr(val, "REML") <- REML
816 :     class(val) <- "logLik"
817 :     val
818 :     })
819 :    
820 : bates 446 setMethod("anova", signature(object = "lmer"),
821 :     function(object, ...)
822 :     {
823 :     mCall <- match.call(expand.dots = TRUE)
824 :     dots <- list(...)
825 :     modp <- logical(0)
826 :     if (length(dots))
827 :     modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm")
828 :     if (any(modp)) { # multiple models - form table
829 :     opts <- dots[!modp]
830 :     mods <- c(list(object), dots[modp])
831 :     names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)
832 :     mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE), attr, "df"))]
833 :     calls <- lapply(mods, slot, "call")
834 :     data <- lapply(calls, "[[", "data")
835 :     if (any(data != data[[1]])) stop("all models must be fit to the same data object")
836 :     header <- paste("Data:", data[[1]])
837 :     subset <- lapply(calls, "[[", "subset")
838 :     if (any(subset != subset[[1]])) stop("all models must use the same subset")
839 :     if (!is.null(subset[[1]]))
840 :     header <-
841 :     c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
842 :     llks <- lapply(mods, logLik, REML = FALSE)
843 :     Df <- sapply(llks, attr, "df")
844 :     llk <- unlist(llks)
845 :     chisq <- 2 * pmax(0, c(NA, diff(llk)))
846 :     dfChisq <- c(NA, diff(Df))
847 :     val <- data.frame(Df = Df,
848 :     AIC = sapply(llks, AIC),
849 :     BIC = sapply(llks, BIC),
850 :     logLik = llk,
851 :     "Chisq" = chisq,
852 :     "Chi Df" = dfChisq,
853 :     "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
854 :     check.names = FALSE)
855 :     class(val) <- c("anova", class(val))
856 :     attr(val, "heading") <-
857 : bates 690 c(header, "Models:",
858 : bates 446 paste(names(mods),
859 :     unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
860 : bates 690 sep = ": "))
861 : bates 446 return(val)
862 :     } else {
863 : bates 571 foo <- object
864 :     foo@status["factored"] <- FALSE
865 :     .Call("lmer_factor", foo, PACKAGE="Matrix")
866 :     dfr <- getFixDF(foo)
867 :     rcol <- ncol(foo@RXX)
868 :     ss <- foo@RXX[ , rcol]^2
869 :     ssr <- ss[[rcol]]
870 :     ss <- ss[seq(along = dfr)]
871 :     names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
872 :     asgn <- foo@assign
873 :     terms <- foo@terms
874 :     nmeffects <- attr(terms, "term.labels")
875 :     if ("(Intercept)" %in% names(ss))
876 :     nmeffects <- c("(Intercept)", nmeffects)
877 :     ss <- unlist(lapply(split(ss, asgn), sum))
878 :     df <- unlist(lapply(split(asgn, asgn), length))
879 :     dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
880 :     ms <- ss/df
881 :     f <- ms/(ssr/dfr)
882 :     P <- pf(f, df, dfr, lower.tail = FALSE)
883 :     table <- data.frame(df, ss, ms, dfr, f, P)
884 :     dimnames(table) <-
885 :     list(nmeffects,
886 :     c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
887 :     if ("(Intercept)" %in% nmeffects) table <- table[-1,]
888 :     attr(table, "heading") <- "Analysis of Variance Table"
889 :     class(table) <- c("anova", "data.frame")
890 :     table
891 : bates 446 }
892 : bates 316 })
893 : bates 446
894 :     setMethod("update", signature(object = "lmer"),
895 :     function(object, formula., ..., evaluate = TRUE)
896 :     {
897 :     call <- object@call
898 :     if (is.null(call))
899 :     stop("need an object with call component")
900 :     extras <- match.call(expand.dots = FALSE)$...
901 :     if (!missing(formula.))
902 :     call$formula <- update.formula(formula(object), formula.)
903 :     if (length(extras) > 0) {
904 :     existing <- !is.na(match(names(extras), names(call)))
905 :     for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
906 :     if (any(!existing)) {
907 :     call <- c(as.list(call), extras[!existing])
908 :     call <- as.call(call)
909 :     }
910 :     }
911 :     if (evaluate)
912 :     eval(call, parent.frame())
913 :     else call
914 :     })
915 :    
916 :    
917 :     setMethod("confint", signature(object = "lmer"),
918 :     function (object, parm, level = 0.95, ...)
919 :     {
920 :     cf <- fixef(object)
921 :     pnames <- names(cf)
922 :     if (missing(parm))
923 :     parm <- seq(along = pnames)
924 :     else if (is.character(parm))
925 :     parm <- match(parm, pnames, nomatch = 0)
926 :     a <- (1 - level)/2
927 :     a <- c(a, 1 - a)
928 :     pct <- paste(round(100 * a, 1), "%")
929 :     ci <- array(NA, dim = c(length(parm), 2),
930 :     dimnames = list(pnames[parm], pct))
931 :     ses <- sqrt(diag(vcov(object)))[parm]
932 : bates 449 ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
933 : bates 446 ci
934 :     })
935 :    
936 : bates 689 setMethod("param", signature(object = "lmer"),
937 : bates 449 function(object, unconst = FALSE, ...) {
938 :     .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
939 :     })
940 : bates 446
941 : bates 449 setMethod("deviance", "lmer",
942 :     function(object, REML = NULL, ...) {
943 :     .Call("lmer_factor", object, PACKAGE = "Matrix")
944 :     if (is.null(REML))
945 :     REML <- if (length(oR <- object@REML)) oR else FALSE
946 :     object@deviance[[ifelse(REML, "REML", "ML")]]
947 :     })
948 : bates 446
949 : deepayan 721
950 :     setMethod("deviance", "glmer",
951 :     function(object, ...) {
952 :     -2 * object@glmmll
953 :     })
954 :    
955 : bates 449 setMethod("chol", signature(x = "lmer"),
956 :     function(x, pivot = FALSE, LINPACK = pivot) {
957 :     x@status["factored"] <- FALSE # force a decomposition
958 :     .Call("lmer_factor", x, PACKAGE = "Matrix")
959 :     })
960 :    
961 :     setMethod("solve", signature(a = "lmer", b = "missing"),
962 :     function(a, b, ...)
963 : bates 562 .Call("lmer_invert", a, PACKAGE = "Matrix")
964 : bates 449 )
965 :    
966 :     setMethod("formula", "lmer", function(x, ...) x@call$formula)
967 :    
968 :     setMethod("vcov", signature(object = "lmer"),
969 :     function(object, REML = object@REML, useScale = TRUE,...) {
970 :     sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
971 :     rr <- object@RXX
972 :     nms <- object@cnames[[".fixed"]]
973 :     dimnames(rr) <- list(nms, nms)
974 :     nr <- nrow(rr)
975 :     rr <- rr[-nr, -nr, drop = FALSE]
976 :     rr <- rr %*% t(rr)
977 :     if (useScale) {
978 :     rr = sc^2 * rr
979 :     }
980 :     rr
981 :     })
982 :    
983 : bates 550 ## Extract the L matrix
984 :     setAs("lmer", "dtTMatrix",
985 :     function(from)
986 :     {
987 :     ## force a refactorization if the factors have been inverted
988 :     if (from@status["inverted"]) from@status["factored"] <- FALSE
989 :     .Call("lmer_factor", from, PACKAGE = "Matrix")
990 :     L <- lapply(from@L, as, "dgTMatrix")
991 :     nf <- length(from@D)
992 :     Gp <- from@Gp
993 :     nL <- Gp[nf + 1]
994 : bates 562 Li <- integer(0)
995 :     Lj <- integer(0)
996 :     Lx <- double(0)
997 : bates 550 for (i in 1:nf) {
998 :     for (j in 1:i) {
999 :     Lij <- L[[Lind(i, j)]]
1000 : bates 562 Li <- c(Li, Lij@i + Gp[i])
1001 :     Lj <- c(Lj, Lij@j + Gp[j])
1002 :     Lx <- c(Lx, Lij@x)
1003 : bates 550 }
1004 :     }
1005 : bates 562 new("dtTMatrix", Dim = as.integer(c(nL, nL)), i = Li, j = Lj, x = Lx,
1006 : bates 550 uplo = "L", diag = "U")
1007 :     })
1008 : bates 562
1009 :     ## Extract the ZZX matrix
1010 :     setAs("lmer", "dsTMatrix",
1011 :     function(from)
1012 :     {
1013 :     .Call("lmer_inflate", from, PACKAGE = "Matrix")
1014 :     ZZpO <- lapply(from@ZZpO, as, "dgTMatrix")
1015 :     ZZ <- lapply(from@ZtZ, as, "dgTMatrix")
1016 :     nf <- length(ZZpO)
1017 :     Gp <- from@Gp
1018 :     nZ <- Gp[nf + 1]
1019 :     Zi <- integer(0)
1020 :     Zj <- integer(0)
1021 :     Zx <- double(0)
1022 :     for (i in 1:nf) {
1023 :     ZZpOi <- ZZpO[[i]]
1024 :     Zi <- c(Zi, ZZpOi@i + Gp[i])
1025 :     Zj <- c(Zj, ZZpOi@j + Gp[i])
1026 :     Zx <- c(Zx, ZZpOi@x)
1027 :     if (i > 1) {
1028 :     for (j in 1:(i-1)) {
1029 :     ZZij <- ZZ[[Lind(i, j)]]
1030 :     ## off-diagonal blocks are transposed
1031 :     Zi <- c(Zi, ZZij@j + Gp[j])
1032 :     Zj <- c(Zj, ZZij@i + Gp[i])
1033 :     Zx <- c(Zx, ZZij@x)
1034 :     }
1035 :     }
1036 :     }
1037 :     new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx,
1038 :     uplo = "U")
1039 :     })
1040 : bates 689
1041 :     setMethod("fitted", signature(object = "lmer"),
1042 : bates 691 function(object, ...)
1043 :     napredict(attr(object@frame, "na.action"), object@fitted))
1044 : bates 689
1045 :     setMethod("residuals", signature(object = "lmer"),
1046 : bates 691 function(object, ...)
1047 :     naresid(attr(object@frame, "na.action"), object@residuals))
1048 : bates 689
1049 :     setMethod("resid", signature(object = "lmer"),
1050 :     function(object, ...) do.call("residuals", c(list(object), list(...))))
1051 :    
1052 :     setMethod("coef", signature(object = "lmer"),
1053 :     function(object, ...)
1054 :     {
1055 :     fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
1056 :     ref <- as(ranef(object), "list")
1057 :     names(ref) <- names(object@flist)
1058 :     val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
1059 :     for (i in seq(a = val)) {
1060 :     refi <- ref[[i]]
1061 :     row.names(val[[i]]) <- row.names(refi)
1062 :     if (!all(names(refi) %in% names(fef)))
1063 :     stop("unable to align random and fixed effects")
1064 :     val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
1065 :     }
1066 :     new("lmer.coef", val)
1067 :     })
1068 :    
1069 :     setMethod("plot", signature(x = "lmer.coef"),
1070 :     function(x, y, ...)
1071 :     {
1072 :     varying <- unique(do.call("c",
1073 :     lapply(x, function(el)
1074 :     names(el)[sapply(el,
1075 :     function(col)
1076 :     any(col != col[1]))])))
1077 :     gf <- do.call("rbind", lapply(x, "[", j = varying))
1078 :     gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
1079 :     switch(min(length(varying), 3),
1080 :     qqmath(eval(substitute(~ x | .grp,
1081 :     list(x = as.name(varying[1])))), gf, ...),
1082 :     xyplot(eval(substitute(y ~ x | .grp,
1083 :     list(y = as.name(varying[1]),
1084 :     x = as.name(varying[2])))), gf, ...),
1085 :     splom(~ gf | .grp, ...))
1086 :     })
1087 :    
1088 :     setMethod("plot", signature(x = "lmer.ranef"),
1089 :     function(x, y, ...)
1090 :     {
1091 :     lapply(x, function(x) {
1092 :     cn <- lapply(colnames(x), as.name)
1093 :     switch(min(ncol(x), 3),
1094 :     qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
1095 :     xyplot(eval(substitute(y ~ x, list(y = cn[[1]], x = cn[[2]]))),
1096 :     x, ...),
1097 :     splom(~ x, ...))
1098 :     })
1099 :     })
1100 :    
1101 :     setMethod("with", signature(data = "lmer"),
1102 : bates 690 function(data, expr, ...) {
1103 : bates 691 dat <- eval(data@call$data)
1104 :     if (!is.null(na.act <- attr(data@frame, "na.action")))
1105 :     dat <- dat[-na.act, ]
1106 :     lst <- c(list(. = data), data@flist, data@frame, dat)
1107 :     eval(substitute(expr), lst[unique(names(lst))])
1108 :     })
1109 : bates 690
1110 : bates 691 setMethod("terms", signature(x = "lmer"),
1111 :     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