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 1175 - (view) (download)

1 : bates 767 # Methods for lmer and for the objects that it produces
2 : bates 689
3 :     ## Some utilities
4 :    
5 : bates 775 ## Return the pairs of expressions separated by vertical bars
6 : bates 769 findbars <- function(term)
7 :     {
8 :     if (is.name(term) || is.numeric(term)) return(NULL)
9 :     if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
10 :     if (!is.call(term)) stop("term must be of class call")
11 :     if (term[[1]] == as.name('|')) return(term)
12 :     if (length(term) == 2) return(findbars(term[[2]]))
13 :     c(findbars(term[[2]]), findbars(term[[3]]))
14 :     }
15 :    
16 : bates 775 ## Return the formula omitting the pairs of expressions
17 :     ## that are separated by vertical bars
18 : bates 769 nobars <- function(term)
19 :     {
20 : bates 1150 if (!('|' %in% all.names(term))) return(term)
21 : bates 769 if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
22 :     if (length(term) == 2) {
23 :     nb <- nobars(term[[2]])
24 :     if (is.null(nb)) return(NULL)
25 :     term[[2]] <- nb
26 :     return(term)
27 :     }
28 :     nb2 <- nobars(term[[2]])
29 :     nb3 <- nobars(term[[3]])
30 :     if (is.null(nb2)) return(nb3)
31 :     if (is.null(nb3)) return(nb2)
32 :     term[[2]] <- nb2
33 :     term[[3]] <- nb3
34 :     term
35 :     }
36 :    
37 : bates 775 ## Substitute the '+' function for the '|' function
38 : bates 769 subbars <- function(term)
39 :     {
40 :     if (is.name(term) || is.numeric(term)) return(term)
41 :     if (length(term) == 2) {
42 :     term[[2]] <- subbars(term[[2]])
43 :     return(term)
44 :     }
45 :     stopifnot(length(term) == 3)
46 : maechler 832 if (is.call(term) && term[[1]] == as.name('|'))
47 :     term[[1]] <- as.name('+')
48 : bates 769 term[[2]] <- subbars(term[[2]])
49 :     term[[3]] <- subbars(term[[3]])
50 :     term
51 :     }
52 : bates 824
53 : bates 979 ## Expand an expression with colons to the sum of the lhs
54 :     ## and the current expression
55 :     colExpand <- function(term)
56 :     {
57 :     if (is.name(term) || is.numeric(term)) return(term)
58 :     if (length(term) == 2) {
59 :     term[[2]] <- colExpand(term[[2]])
60 :     return(term)
61 :     }
62 :     stopifnot(length(term) == 3)
63 :     if (is.call(term) && term[[1]] == as.name(':')) {
64 :     return(substitute(A+B, list(A = term, B = colExpand(term[[2]]))))
65 :     }
66 :     term[[2]] <- colExpand(term[[2]])
67 :     term[[3]] <- colExpand(term[[3]])
68 :     term
69 :     }
70 :    
71 : bates 824 abbrvNms <- function(gnm, cnms)
72 :     {
73 :     ans <- paste(abbreviate(gnm), abbreviate(cnms), sep = '.')
74 :     if (length(cnms) > 1) {
75 :     anms <- lapply(cnms, abbreviate, minlength = 3)
76 :     nmmat <- outer(anms, anms, paste, sep = '.')
77 :     ans <- c(ans, paste(abbreviate(gnm, minlength = 3),
78 :     nmmat[upper.tri(nmmat)], sep = '.'))
79 :     }
80 :     ans
81 :     }
82 :    
83 : bates 775 ## Control parameters for lmer
84 :     lmerControl <-
85 : maechler 832 function(maxIter = 200, # used in ../src/lmer.c only
86 : bates 888 tolerance = sqrt(.Machine$double.eps),# ditto
87 : bates 769 msMaxIter = 200,
88 : maechler 832 ## msTol = sqrt(.Machine$double.eps),
89 :     ## FIXME: should be able to pass tolerances to nlminb()
90 :     msVerbose = getOption("verbose"),
91 : bates 752 niterEM = 15,
92 : bates 435 EMverbose = getOption("verbose"),
93 : maechler 843 PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead
94 : bates 1166 gradient = TRUE,
95 :     Hessian = FALSE # unused _FIXME_
96 : maechler 832 )
97 : bates 435 {
98 : bates 775 list(maxIter = as.integer(maxIter),
99 : maechler 832 tolerance = as.double(tolerance),
100 : bates 775 msMaxIter = as.integer(msMaxIter),
101 : maechler 832 ## msTol = as.double(msTol),
102 :     msVerbose = as.integer(msVerbose),# "integer" on purpose
103 : bates 775 niterEM = as.integer(niterEM),
104 : maechler 832 EMverbose = as.logical(EMverbose),
105 : bates 775 PQLmaxIt = as.integer(PQLmaxIt),
106 : bates 1166 gradient = as.logical(gradient),
107 :     Hessian = as.logical(Hessian))
108 : bates 435 }
109 :    
110 : bates 1166 rWishart <- function(n, df, invScal)
111 :     .Call("Matrix_rWishart", n, df, invScal, PACKAGE = "Matrix")
112 :    
113 :     setMethod("coef", signature(object = "mer"),
114 : bates 1150 function(object, ...)
115 :     {
116 : bates 1166 fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
117 :     ref <- ranef(object)
118 :     val <- lapply(ref, function(x) fef[rep(1, nrow(x)),,drop = FALSE])
119 :     for (i in seq(a = val)) {
120 :     refi <- ref[[i]]
121 :     row.names(val[[i]]) <- row.names(refi)
122 :     nmsi <- colnames(refi)
123 :     if (!all(nmsi %in% names(fef)))
124 :     stop("unable to align random and fixed effects")
125 :     for (nm in nmsi) val[[i]][[nm]] <- val[[i]][[nm]] + refi[,nm]
126 :     }
127 :     val
128 :     })
129 : bates 1150
130 :     ## setMethod("plot", signature(x = "lmer.coef"),
131 :     ## function(x, y, ...)
132 :     ## {
133 :     ## varying <- unique(do.call("c",
134 :     ## lapply(x, function(el)
135 :     ## names(el)[sapply(el,
136 :     ## function(col)
137 :     ## any(col != col[1]))])))
138 :     ## gf <- do.call("rbind", lapply(x, "[", j = varying))
139 :     ## gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
140 :     ## switch(min(length(varying), 3),
141 :     ## qqmath(eval(substitute(~ x | .grp,
142 :     ## list(x = as.name(varying[1])))), gf, ...),
143 :     ## xyplot(eval(substitute(y ~ x | .grp,
144 :     ## list(y = as.name(varying[1]),
145 :     ## x = as.name(varying[2])))), gf, ...),
146 :     ## splom(~ gf | .grp, ...))
147 :     ## })
148 :    
149 : bates 1166 setMethod("plot", signature(x = "lmer.ranef"),
150 :     function(x, y, ...)
151 :     {
152 :     lapply(x, function(x) {
153 :     cn <- lapply(colnames(x), as.name)
154 :     switch(min(ncol(x), 3),
155 :     qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
156 :     xyplot(eval(substitute(y ~ x,
157 :     list(y = cn[[1]],
158 :     x = cn[[2]]))), x, ...),
159 :     splom(~ x, ...))
160 :     })
161 :     })
162 : bates 1150
163 :     setMethod("with", signature(data = "lmer"),
164 :     function(data, expr, ...) {
165 :     dat <- eval(data@call$data)
166 :     if (!is.null(na.act <- attr(data@frame, "na.action")))
167 :     dat <- dat[-na.act, ]
168 :     lst <- c(list(. = data), data@flist, data@frame, dat)
169 :     eval(substitute(expr), lst[unique(names(lst))])
170 :     })
171 :    
172 :     setMethod("terms", signature(x = "lmer"),
173 :     function(x, ...) x@terms)
174 :    
175 :    
176 : bates 755 setMethod("lmer", signature(formula = "formula"),
177 : bates 1175 function(formula, data, family = gaussian,
178 : bates 689 method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
179 : bates 901 control = list(), start,
180 : bates 1175 subset, weights, na.action, offset, contrasts = NULL,
181 :     model = TRUE, x = TRUE, y = TRUE , ...)
182 : maechler 832 {
183 :     ## match and check parameters
184 : bates 755 if (length(formula) < 3) stop("formula must be a two-sided formula")
185 :     cv <- do.call("lmerControl", control)
186 : bates 1150
187 :     ## Must evaluate the model frame first and then fit the glm using
188 :     ## that frame. Otherwise missing values in the grouping factors
189 :     ## cause inconsistent numbers of observations.
190 : maechler 832 mf <- match.call()
191 : bates 1175 m <- match(c("data", "subset", "weights",
192 : bates 755 "na.action", "offset"), names(mf), 0)
193 : bates 1175 mf <- mf[c(1, m)]
194 : bates 755 frame.form <- subbars(formula) # substitute `+' for `|'
195 :     fixed.form <- nobars(formula) # remove any terms with `|'
196 : bates 767 if (inherits(fixed.form, "name")) # RHS is empty - use a constant
197 : bates 755 fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
198 :     environment(fixed.form) <- environment(frame.form) <- environment(formula)
199 : bates 1150
200 :     ## evaluate a model frame for fixed and random effects
201 :     mf$formula <- frame.form
202 :     mf$drop.unused.levels <- TRUE
203 :     mf[[1]] <- as.name("model.frame")
204 : bates 1175 fe <- mf
205 :     mf <- eval(mf, parent.frame())
206 : bates 1150
207 : bates 1175 ## get the terms for the fixed-effects only
208 :     fe$formula <- fixed.form
209 :     fe <- eval(fe, parent.frame())
210 :     mt <- attr(fe, "terms") # allow model.frame to update them
211 :     ## response vector
212 :     Y <- model.response(mf, "numeric")
213 :     ## avoid problems with 1D arrays, but keep names
214 :     if(length(dim(Y)) == 1) {
215 :     nm <- rownames(Y)
216 :     dim(Y) <- NULL
217 :     if(!is.null(nm)) names(Y) <- nm
218 :     }
219 :     ## null model support
220 :     X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts) else matrix(,NROW(Y),0)
221 :    
222 :     weights <- model.weights(mf)
223 :     offset <- model.offset(mf)
224 :     ## check weights and offset
225 :     if( !is.null(weights) && any(weights < 0) )
226 :     stop("negative weights not allowed")
227 :     if(!is.null(offset) && length(offset) != NROW(Y))
228 :     stop(gettextf("number of offsets is %d should equal %d (number of observations)", length(offset), NROW(Y)), domain = NA)
229 :     if (is.null(weights)) weights <- rep.int(1, length(Y))
230 :     if (is.null(offset)) offset <- numeric(length(Y))
231 :    
232 : bates 1150 ## fit a glm model to the fixed formula
233 : bates 1175 ## fe$formula <- fixed.form
234 :     ## fe$subset <- NULL # subset has already been created in call to data.frame
235 :     ## fe$data <- frm
236 :     ## fe$x <- fe$model <- fe$y <- TRUE
237 :     ## fe[[1]] <- as.name("glm")
238 :     ## glmFit <- eval(fe, parent.frame())
239 :     ## x <- glmFit$x
240 :     ## y <- as.double(glmFit$y)
241 : bates 1150
242 : bates 1175 if(is.character(family))
243 :     family <- get(family, mode = "function", envir = parent.frame())
244 :     if(is.function(family)) family <- family()
245 :     if(is.null(family$family)) {
246 :     print(family)
247 :     stop("'family' not recognized")
248 :     }
249 : bates 939 ## check for a linear mixed model
250 :     lmm <- family$family == "gaussian" && family$link == "identity"
251 : maechler 832 if (lmm) { # linear mixed model
252 :     method <- match.arg(method)
253 :     if (method %in% c("PQL", "Laplace", "AGQ")) {
254 :     warning(paste('Argument method = "', method,
255 :     '" is not meaningful for a linear mixed model.\n',
256 :     'Using method = "REML".\n', sep = ''))
257 :     method <- "REML"
258 :     }
259 :     } else { # generalized linear mixed model
260 :     if (missing(method)) method <- "PQL"
261 :     else {
262 :     method <- match.arg(method)
263 :     if (method == "ML") method <- "PQL"
264 :     if (method == "REML")
265 :     warning('Argument method = "REML" is not meaningful ',
266 :     'for a generalized linear mixed model.',
267 :     '\nUsing method = "PQL".\n')
268 :     }
269 :     }
270 : bates 1166 if (method == "AGQ")
271 :     stop('method = "AGQ" not yet implemented for supernodal representation')
272 : bates 1150 ## create factor list for the random effects
273 : bates 435 bars <- findbars(formula[[3]])
274 : bates 1150 names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
275 :     fl <- lapply(bars,
276 :     function(x)
277 :     eval(substitute(as.factor(fac)[,drop = TRUE],
278 : bates 1175 list(fac = x[[3]])), mf))
279 : bates 435 ## order factor list by decreasing number of levels
280 : bates 1150 nlev <- sapply(fl, function(x) length(levels(x)))
281 : bates 452 if (any(diff(nlev) > 0)) {
282 : bates 1150 ord <- rev(order(nlev))
283 :     bars <- bars[ord]
284 :     fl <- fl[ord]
285 : bates 435 }
286 : bates 1150 ## create list of transposed model matrices for random effects
287 :     Ztl <- lapply(bars, function(x)
288 :     t(model.matrix(eval(substitute(~ expr,
289 :     list(expr = x[[2]]))),
290 : bates 1175 mf)))
291 : bates 1150 ## Create the mixed-effects representation (mer) object
292 :     mer <- .Call("mer_create", fl,
293 :     .Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"),
294 : bates 1175 X, Y, method, sapply(Ztl, nrow),
295 :     c(lapply(Ztl, rownames), list(.fixed = colnames(X))),
296 : bates 1150 !(family$family %in% c("binomial", "poisson")),
297 :     match.call(), family,
298 :     PACKAGE = "Matrix")
299 :     if (lmm) {
300 :     .Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
301 : bates 755 LMEoptimize(mer) <- cv
302 : bates 1166 return(new("lmer", mer,
303 : bates 1175 frame = if (model) mf else data.frame(),
304 :     terms = mt))
305 : bates 755 }
306 :    
307 :     ## The rest of the function applies to generalized linear mixed models
308 :     gVerb <- getOption("verbose")
309 : bates 1175 eta <- glm.fit(X, Y, weights = weights, offset = offset, family = family,
310 :     intercept = attr(mt, "intercept") > 0))$fitted.values
311 :     wtssqr <- weights * weights
312 : bates 767 if (is.null(offset)) offset <- numeric(length(eta))
313 :     linkinv <- quote(family$linkinv(eta))
314 :     mu.eta <- quote(family$mu.eta(eta))
315 : bates 1150 mu <- family$linkinv(eta)
316 : bates 767 variance <- quote(family$variance(mu))
317 : bates 1150 dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
318 : bates 775 LMEopt <- get("LMEoptimize<-")
319 :     doLMEopt <- quote(LMEopt(x = mer, value = cv))
320 : bates 767
321 : bates 809 GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
322 :     .Call("glmer_PQL", GSpt, PACKAGE = "Matrix") # obtain PQL estimates
323 : bates 1166 PQLpars <- c(fixef(mer),
324 :     .Call("mer_coef", mer, 2, PACKAGE = "Matrix"))
325 :     if (method == "PQL") {
326 :     .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")
327 :     .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
328 :     return(new("lmer", mer,
329 : bates 1175 frame = if (model) mf else data.frame(),
330 :     terms = mt))
331 : bates 1166 }
332 :    
333 : bates 774 fixInd <- seq(ncol(x))
334 :     ## pars[fixInd] == beta, pars[-fixInd] == theta
335 : bates 777 ## indicator of constrained parameters
336 : bates 1150 const <- c(rep(FALSE, length(fixInd)),
337 :     unlist(lapply(mer@nc[seq(along = fl)],
338 : bates 777 function(k) 1:((k*(k+1))/2) <= k)
339 :     ))
340 : bates 1150 devLaplace <- function(pars)
341 :     .Call("glmer_devLaplace", pars, GSpt, PACKAGE = "Matrix")
342 : maechler 832
343 : bates 1150 optimRes <-
344 :     nlminb(PQLpars, devLaplace,
345 :     lower = ifelse(const, 5e-10, -Inf),
346 : bates 1166 control = list(trace = cv$msVerbose,
347 : bates 1150 iter.max = cv$msMaxIter))
348 :     .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
349 : bates 1166 return(new("lmer", mer,
350 : bates 1175 frame = if (model) mf else data.frame(),
351 :     terms = mt))
352 : bates 1150
353 : bates 435 })
354 :    
355 : bates 1150 ## Extract the L matrix
356 :     setAs("mer", "dtCMatrix", function(from)
357 :     .Call("mer_dtCMatrix", from, PACKAGE = "Matrix"))
358 :    
359 :     ## Extract the fixed effects
360 :     setMethod("fixef", signature(object = "mer"),
361 :     function(object, ...)
362 :     .Call("mer_fixef", object, PACKAGE = "Matrix"))
363 :    
364 :     ## Extract the random effects
365 :     setMethod("ranef", signature(object = "mer"),
366 :     function(object, ...)
367 : bates 1166 new("lmer.ranef", .Call("mer_ranef", object, PACKAGE = "Matrix"))
368 : bates 1150 )
369 :    
370 :     ## Optimization for mer objects
371 : bates 755 setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
372 : bates 316 function(x, value)
373 :     {
374 :     if (value$msMaxIter < 1) return(x)
375 :     nc <- x@nc
376 : bates 1150 constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
377 : bates 752 fn <- function(pars)
378 : bates 1150 deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
379 : bates 1166 gr <- if (value$gradient)
380 : bates 1150 function(pars) {
381 :     if (!isTRUE(all.equal(pars,
382 :     .Call("mer_coef", x,
383 :     2, PACKAGE = "Matrix"))))
384 :     .Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")
385 :     .Call("mer_gradient", x, 2, PACKAGE = "Matrix")
386 :     }
387 :     else NULL
388 :     optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"),
389 :     fn, gr,
390 :     lower = ifelse(constr, 5e-10, -Inf),
391 :     control = list(iter.max = value$msMaxIter,
392 :     trace = as.integer(value$msVerbose)))
393 :     .Call("mer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
394 : bates 316 if (optimRes$convergence != 0) {
395 : bates 1150 warning(paste("nlminb returned message",
396 : bates 777 optimRes$message,"\n"))
397 : bates 316 }
398 :     return(x)
399 :     })
400 :    
401 : bates 1150 setMethod("deviance", signature(object = "mer"),
402 :     function(object, ...) {
403 :     .Call("mer_factor", object, PACKAGE = "Matrix")
404 :     object@deviance[[ifelse(object@method == "REML", "REML", "ML")]]
405 : bates 316 })
406 :    
407 : bates 1150 setMethod("mcmcsamp", signature(object = "mer"),
408 :     function(object, n = 1, verbose = FALSE, saveb = FALSE,
409 :     trans = TRUE, ...)
410 :     {
411 : bates 1171 family <- object@family
412 : bates 1170 lmm <- family$family == "gaussian" && family$link == "identity"
413 :     if (!lmm)
414 :     stop("mcmcsamp for GLMMs not yet implemented in supernodal representation")
415 : bates 1150 ans <- t(.Call("mer_MCMCsamp", object, saveb, n,
416 :     trans, PACKAGE = "Matrix"))
417 :     attr(ans, "mcpar") <- as.integer(c(1, n, 1))
418 :     class(ans) <- "mcmc"
419 :     glmer <- FALSE
420 :     gnms <- names(object@flist)
421 :     cnms <- object@cnames
422 :     ff <- fixef(object)
423 :     colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
424 :     unlist(lapply(seq(along = gnms),
425 :     function(i)
426 :     abbrvNms(gnms[i],cnms[[i]]))))
427 :     if (trans) {
428 :     ## parameter type: 0 => fixed effect, 1 => variance,
429 :     ## 2 => covariance
430 :     ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1,
431 :     unlist(lapply(seq(along = gnms),
432 :     function(i)
433 :     {
434 :     k <- length(cnms[[i]])
435 :     rep(1:2, c(k, (k*(k-1))/2))
436 :     })))
437 :     colnms[ptyp == 1] <-
438 :     paste("log(", colnms[ptyp == 1], ")", sep = "")
439 :     colnms[ptyp == 2] <-
440 :     paste("atanh(", colnms[ptyp == 2], ")", sep = "")
441 :     }
442 :     colnames(ans) <- colnms
443 :     ans
444 :     })
445 : bates 316
446 : bates 1150 setMethod("simulate", signature(object = "mer"),
447 :     function(object, nsim = 1, seed = NULL, ...)
448 :     {
449 :     if(!exists(".Random.seed", envir = .GlobalEnv))
450 :     runif(1) # initialize the RNG if necessary
451 :     if(is.null(seed))
452 :     RNGstate <- .Random.seed
453 :     else {
454 :     R.seed <- .Random.seed
455 :     set.seed(seed)
456 :     RNGstate <- structure(seed, kind = as.list(RNGkind()))
457 :     on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
458 :     }
459 : deepayan 721
460 : bates 1150 family <- object@family
461 :     if (family$family != "gaussian" ||
462 :     family$link != "identity")
463 :     stop("simulation of generalized linear mixed models not yet implemented")
464 :     ## similate the linear predictors
465 :     lpred <- .Call("mer_simulate", object, nsim, PACKAGE = "Matrix")
466 :     sc <- 1
467 : maechler 1165 if (object@useScale)
468 : bates 1150 sc <- .Call("mer_sigma", object, object@method == "REML",
469 :     PACKAGE = "Matrix")
470 :     ## add fixed-effects contribution and per-observation noise term
471 :     lpred <- as.data.frame(lpred + drop(object@X %*% fixef(object)) +
472 :     rnorm(prod(dim(lpred)), sd = sc))
473 :     ## save the seed
474 :     attr(lpred, "seed") <- RNGstate
475 :     lpred
476 :     })
477 : bates 316
478 :    
479 : bates 1150 setMethod("show", "mer",
480 : bates 316 function(object) {
481 : bates 1150 vcShow <- function(varc, useScale)
482 :     {
483 :     digits <- max(3, getOption("digits") - 2)
484 :     sc <- attr(varc, "sc")
485 :     recorr <- lapply(varc, function(el) el@factors$correlation)
486 :     reStdDev <- c(lapply(recorr, slot, "sd"), list(Residual = sc))
487 :     reLens <- unlist(c(lapply(reStdDev, length)))
488 :     reMat <- array('', c(sum(reLens), 4),
489 :     list(rep('', sum(reLens)),
490 :     c("Groups", "Name", "Variance", "Std.Dev.")))
491 :     reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
492 :     reMat[,2] <- c(unlist(lapply(reStdDev, names)), "")
493 :     reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
494 :     reMat[,4] <- format(unlist(reStdDev), digits = digits)
495 :     if (any(reLens > 1)) {
496 :     maxlen <- max(reLens)
497 :     corr <-
498 :     do.call("rbind",
499 :     lapply(recorr,
500 :     function(x, maxlen) {
501 :     x <- as(x, "matrix")
502 :     cc <- format(round(x, 3), nsmall = 3)
503 :     cc[!lower.tri(cc)] <- ""
504 :     nr <- dim(cc)[1]
505 :     if (nr >= maxlen) return(cc)
506 :     cbind(cc, matrix("", nr, maxlen-nr))
507 :     }, maxlen))
508 :     colnames(corr) <- c("Corr", rep("", maxlen - 1))
509 :     reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr))))
510 :     }
511 :     if (!useScale) reMat <- reMat[-nrow(reMat),]
512 :     print(reMat, quote = FALSE)
513 :     }
514 : maechler 1165
515 : bates 1150 fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")
516 : bates 449 useScale <- object@useScale
517 : bates 1150 corF <- vcov(object)@factors$correlation
518 : bates 1166 #DF <- getFixDF(object)
519 :     coefs <- cbind(fcoef, corF@sd)#, DF)
520 : bates 316 dimnames(coefs) <-
521 : bates 1166 list(names(fcoef), c("Estimate", "Std. Error"))#, "DF"))
522 : bates 1123 digits <- max(3, getOption("digits") - 2)
523 : bates 755 REML <- object@method == "REML"
524 : bates 1150 llik <- logLik(object, REML)
525 : bates 449 dev <- object@deviance
526 : bates 1150 devc <- object@devComp
527 : maechler 1165
528 : bates 449 rdig <- 5
529 : bates 727 if (glz <- !(object@method %in% c("REML", "ML"))) {
530 :     cat(paste("Generalized linear mixed model fit using",
531 :     object@method, "\n"))
532 :     } else {
533 :     cat("Linear mixed-effects model fit by ")
534 : bates 755 cat(if(REML) "REML\n" else "maximum likelihood\n")
535 : bates 727 }
536 : bates 449 if (!is.null(object@call$formula)) {
537 :     cat("Formula:", deparse(object@call$formula),"\n")
538 :     }
539 :     if (!is.null(object@call$data)) {
540 :     cat(" Data:", deparse(object@call$data), "\n")
541 :     }
542 :     if (!is.null(object@call$subset)) {
543 :     cat(" Subset:",
544 :     deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
545 :     }
546 : bates 727 if (glz) {
547 : bates 750 cat(" Family: ", object@family$family, "(",
548 :     object@family$link, " link)\n", sep = "")
549 : bates 727 print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
550 : bates 1150 logLik = c(llik),
551 :     deviance = -2*llik,
552 :     row.names = ""))
553 : bates 727 } else {
554 :     print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
555 : bates 1150 logLik = c(llik),
556 :     MLdeviance = dev["ML"],
557 :     REMLdeviance = dev["REML"],
558 :     row.names = ""))
559 : bates 727 }
560 : bates 449 cat("Random effects:\n")
561 : bates 1150 vcShow(VarCorr(object), useScale)
562 : bates 449 ngrps <- lapply(object@flist, function(x) length(levels(x)))
563 : bates 1150 cat(sprintf("# of obs: %d, groups: ", devc[1]))
564 : bates 449 cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
565 :     cat("\n")
566 :     if (!useScale)
567 :     cat("\nEstimated scale (compare to 1) ",
568 : bates 1150 .Call("mer_sigma", object, FALSE, PACKAGE = "Matrix"),
569 : bates 449 "\n")
570 :     if (nrow(coefs) > 0) {
571 : bates 1150 if (useScale) {
572 :     stat <- coefs[,1]/coefs[,2]
573 : bates 1166 #pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
574 : bates 1150 nms <- colnames(coefs)
575 : bates 1166 coefs <- cbind(coefs, stat) #, pval)
576 :     colnames(coefs) <- c(nms, "t value")#, "Pr(>|t|)")
577 : bates 1150 } else {
578 :     coefs <- coefs[, 1:2, drop = FALSE]
579 :     stat <- coefs[,1]/coefs[,2]
580 :     pval <- 2*pnorm(abs(stat), lower = FALSE)
581 :     nms <- colnames(coefs)
582 :     coefs <- cbind(coefs, stat, pval)
583 :     colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
584 : bates 449 }
585 :     cat("\nFixed effects:\n")
586 : bates 1166 printCoefmat(coefs, zap.ind = 3)#, tst.ind = 4)
587 : bates 1150 rn <- rownames(coefs)
588 :     if (!is.null(corF)) {
589 :     p <- ncol(corF)
590 :     if (p > 1) {
591 :     cat("\nCorrelation of Fixed Effects:\n")
592 :     corF <- matrix(format(round(corF@x, 3), nsmall = 3),
593 :     nc = p)
594 :     dimnames(corF) <- list(
595 :     abbreviate(rn, minlen=11),
596 :     abbreviate(rn, minlen=6))
597 :     corF[!lower.tri(corF)] <- ""
598 :     print(corF[-1, -p, drop=FALSE], quote = FALSE)
599 : bates 449 }
600 :     }
601 :     }
602 :     invisible(object)
603 : bates 316 })
604 :    
605 : bates 1150 setMethod("vcov", signature(object = "mer"),
606 :     function(object, REML = object@method == "REML",
607 :     useScale = object@useScale,...) {
608 :     sc <- if (object@useScale) {
609 :     .Call("mer_sigma", object, REML, PACKAGE = "Matrix")
610 :     } else { 1 }
611 :     rr <- as(sc^2 * tcrossprod(solve(object@RXX)), "dpoMatrix")
612 :     rr@factors$correlation <- as(rr, "correlation")
613 :     rr
614 :     })
615 :    
616 : bates 316 ## calculates degrees of freedom for fixed effects Wald tests
617 :     ## This is a placeholder. The answers are generally wrong. It will
618 :     ## be very tricky to decide what a 'right' answer should be with
619 :     ## crossed random effects.
620 :    
621 : bates 1150 setMethod("getFixDF", signature(object="mer"),
622 :     function(object, ...) {
623 :     devc <- object@devComp
624 :     rep(as.integer(devc[1]- devc[2]), devc[2])
625 :     })
626 : bates 316
627 : bates 755 setMethod("logLik", signature(object="mer"),
628 :     function(object, REML = object@method == "REML", ...) {
629 : bates 446 val <- -deviance(object, REML = REML)/2
630 : bates 1150 devc <- as.integer(object@devComp[1:2])
631 :     attr(val, "nall") <- attr(val, "nobs") <- devc[1]
632 :     attr(val, "df") <- abs(devc[2]) +
633 :     length(.Call("mer_coef", object, 0, PACKAGE = "Matrix"))
634 : maechler 832 attr(val, "REML") <- REML
635 : bates 446 class(val) <- "logLik"
636 :     val
637 :     })
638 :    
639 : bates 1150 setMethod("VarCorr", signature(x = "mer"),
640 :     function(x, REML = x@method == "REML", useScale = x@useScale, ...)
641 :     {
642 :     sc <- 1
643 :     if (useScale)
644 :     sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix")
645 :     sc2 <- sc * sc
646 : bates 1168 cnames <- x@cnames
647 :     ans <- x@Omega
648 :     for (i in seq(a=ans)) {
649 :     el <- as(sc2 * solve(ans[[i]]), "dpoMatrix")
650 :     el@Dimnames <- list(cnames[[i]], cnames[[i]])
651 : bates 1150 el@factors$correlation <- as(el, "correlation")
652 : bates 1168 ans[[i]] <- el
653 :     }
654 : bates 1150 attr(ans, "sc") <- sc
655 :     ans
656 :     })
657 : deepayan 721
658 : bates 1150 setMethod("anova", signature(object = "mer"),
659 : bates 446 function(object, ...)
660 :     {
661 :     mCall <- match.call(expand.dots = TRUE)
662 :     dots <- list(...)
663 :     modp <- logical(0)
664 :     if (length(dots))
665 : bates 1166 modp <- sapply(dots, is, "mer") | sapply(dots, is, "lm")
666 : bates 446 if (any(modp)) { # multiple models - form table
667 :     opts <- dots[!modp]
668 :     mods <- c(list(object), dots[modp])
669 :     names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)
670 :     mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE), attr, "df"))]
671 :     calls <- lapply(mods, slot, "call")
672 :     data <- lapply(calls, "[[", "data")
673 :     if (any(data != data[[1]])) stop("all models must be fit to the same data object")
674 :     header <- paste("Data:", data[[1]])
675 :     subset <- lapply(calls, "[[", "subset")
676 :     if (any(subset != subset[[1]])) stop("all models must use the same subset")
677 :     if (!is.null(subset[[1]]))
678 :     header <-
679 :     c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
680 :     llks <- lapply(mods, logLik, REML = FALSE)
681 :     Df <- sapply(llks, attr, "df")
682 :     llk <- unlist(llks)
683 :     chisq <- 2 * pmax(0, c(NA, diff(llk)))
684 :     dfChisq <- c(NA, diff(Df))
685 :     val <- data.frame(Df = Df,
686 :     AIC = sapply(llks, AIC),
687 :     BIC = sapply(llks, BIC),
688 :     logLik = llk,
689 :     "Chisq" = chisq,
690 :     "Chi Df" = dfChisq,
691 :     "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
692 :     check.names = FALSE)
693 :     class(val) <- c("anova", class(val))
694 :     attr(val, "heading") <-
695 : bates 690 c(header, "Models:",
696 : bates 446 paste(names(mods),
697 :     unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
698 : bates 690 sep = ": "))
699 : bates 446 return(val)
700 :     } else {
701 : bates 571 foo <- object
702 : bates 1166 #foo@status["factored"] <- FALSE
703 :     #.Call("mer_factor", foo, PACKAGE="Matrix")
704 :     #dfr <- getFixDF(foo)
705 : bates 1150 ss <- foo@rXy^2
706 :     ssr <- exp(foo@devComp["logryy2"])
707 : bates 1166 names(ss) <- object@cnames[[".fixed"]]
708 :     asgn <- attr(foo@X, "assign")
709 : bates 571 terms <- foo@terms
710 :     nmeffects <- attr(terms, "term.labels")
711 :     if ("(Intercept)" %in% names(ss))
712 :     nmeffects <- c("(Intercept)", nmeffects)
713 :     ss <- unlist(lapply(split(ss, asgn), sum))
714 :     df <- unlist(lapply(split(asgn, asgn), length))
715 : bates 1123 #dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
716 : bates 571 ms <- ss/df
717 : bates 1123 #f <- ms/(ssr/dfr)
718 :     #P <- pf(f, df, dfr, lower.tail = FALSE)
719 :     #table <- data.frame(df, ss, ms, dfr, f, P)
720 :     table <- data.frame(df, ss, ms)
721 : bates 571 dimnames(table) <-
722 :     list(nmeffects,
723 : bates 1123 # c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
724 :     c("Df", "Sum Sq", "Mean Sq"))
725 : bates 1166 if ("(Intercept)" %in% nmeffects) table <- table[-match("(Intercept)", nmeffects), ]
726 : bates 571 attr(table, "heading") <- "Analysis of Variance Table"
727 :     class(table) <- c("anova", "data.frame")
728 :     table
729 : bates 446 }
730 : bates 316 })
731 : bates 446
732 : bates 1150 setMethod("confint", signature(object = "mer"),
733 :     function(object, parm, level = 0.95, ...)
734 : maechler 1165 .NotYetImplemented()
735 : bates 1150 )
736 : bates 446
737 : bates 1150 setMethod("fitted", signature(object = "mer"),
738 :     function(object, ...)
739 : bates 1162 .Call("mer_fitted", object, PACKAGE = "Matrix")
740 : bates 1123 )
741 : bates 446
742 : bates 1150 setMethod("formula", signature(x = "mer"),
743 :     function(x, ...)
744 :     x@call$formula
745 : bates 449 )
746 :    
747 : bates 1150 setMethod("residuals", signature(object = "mer"),
748 : bates 691 function(object, ...)
749 : maechler 1165 .NotYetImplemented()
750 : bates 1150 )
751 : bates 689
752 : bates 1150 setMethod("resid", signature(object = "mer"),
753 : bates 691 function(object, ...)
754 : maechler 1165 .NotYetImplemented()
755 : bates 1150 )
756 : bates 689
757 : bates 1150 setMethod("summary", signature(object = "mer"),
758 : bates 1166 function(object, ...) object
759 : bates 1150 )
760 : bates 689
761 : bates 1150 setMethod("update", signature(object = "mer"),
762 : bates 888 function(object, ...)
763 : maechler 1165 .NotYetImplemented()
764 : bates 1150 )
765 : bates 888
766 :    
767 : bates 1150 simss <- function(fm0, fma, nsim)
768 : bates 879 {
769 : bates 1150 ysim <- simulate(fm0, nsim)
770 : bates 1166 cv <- list(gradient = FALSE, msMaxIter = 200:200,
771 : bates 1150 msVerbose = 0:0)
772 :     sapply(ysim, function(yy) {
773 :     .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")
774 :     LMEoptimize(fm0) <- cv
775 :     .Call("mer_update_y", fma, yy, PACKAGE = "Matrix")
776 :     LMEoptimize(fma) <- cv
777 :     exp(c(H0 = fm0@devComp[["logryy2"]],
778 :     Ha = fma@devComp[["logryy2"]]))
779 :     })
780 : bates 879 }
781 : bates 1166
782 :     ## Some leftover code from the old AGQ method in lmer.
783 :     if (FALSE) {
784 :     ### FIXME: For nf == 1 change this to an AGQ evaluation. Needs
785 :     ### AGQ for nc > 1 first.
786 :     fxd <- PQLpars[fixInd]
787 :     loglik <- logLik(mer)
788 :    
789 :     if (method %in% c("Laplace", "AGQ")) {
790 :     nAGQ <- 1
791 :     if (method == "AGQ") { # determine nAGQ at PQL estimates
792 :     dev11 <- devAGQ(PQLpars, 11)
793 :     ## FIXME: Should this be an absolute or a relative tolerance?
794 :     devTol <- sqrt(.Machine$double.eps) * abs(dev11)
795 :     for (nAGQ in c(9, 7, 5, 3, 1))
796 :     if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break
797 :     nAGQ <- nAGQ + 2
798 :     if (gVerb)
799 :     cat(paste("Using", nAGQ, "quadrature points per column\n"))
800 :     }
801 :     obj <- function(pars)
802 :     .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
803 :     optimRes <-
804 :     nlminb(PQLpars, obj,
805 :     lower = ifelse(const, 5e-10, -Inf),
806 :     control = list(trace = getOption("verbose"),
807 :     iter.max = cv$msMaxIter))
808 :     optpars <- optimRes$par
809 :     if (optimRes$convergence != 0)
810 :     warning("nlminb failed to converge")
811 :     deviance <- optimRes$objective
812 :     if (gVerb)
813 :     cat(paste("convergence message", optimRes$message, "\n"))
814 :     fxd[] <- optpars[fixInd] ## preserve the names
815 :     .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
816 :     }
817 :    
818 :     .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
819 :     loglik[] <- -deviance/2
820 :     }

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