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