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

1 : bates 767 # Methods for lmer and for the objects that it produces
2 : bates 689
3 :     ## Some utilities
4 :    
5 : bates 775 ## Return the index into the packed lower triangle
6 : bates 689 Lind <- function(i,j) {
7 :     if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))
8 :     ((i - 1) * i)/2 + j
9 : bates 446 }
10 :    
11 : bates 775 ## Return the pairs of expressions separated by vertical bars
12 : bates 769 findbars <- function(term)
13 :     {
14 :     if (is.name(term) || is.numeric(term)) return(NULL)
15 :     if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
16 :     if (!is.call(term)) stop("term must be of class call")
17 :     if (term[[1]] == as.name('|')) return(term)
18 :     if (length(term) == 2) return(findbars(term[[2]]))
19 :     c(findbars(term[[2]]), findbars(term[[3]]))
20 :     }
21 :    
22 : bates 775 ## Return the formula omitting the pairs of expressions
23 :     ## that are separated by vertical bars
24 : bates 769 nobars <- function(term)
25 :     {
26 :     # FIXME: is the is.name in the condition redundant?
27 :     # A name won't satisfy the first condition.
28 :     if (!('|' %in% all.names(term)) || is.name(term)) return(term)
29 :     if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
30 :     if (length(term) == 2) {
31 :     nb <- nobars(term[[2]])
32 :     if (is.null(nb)) return(NULL)
33 :     term[[2]] <- nb
34 :     return(term)
35 :     }
36 :     nb2 <- nobars(term[[2]])
37 :     nb3 <- nobars(term[[3]])
38 :     if (is.null(nb2)) return(nb3)
39 :     if (is.null(nb3)) return(nb2)
40 :     term[[2]] <- nb2
41 :     term[[3]] <- nb3
42 :     term
43 :     }
44 :    
45 : bates 775 ## Substitute the '+' function for the '|' function
46 : bates 769 subbars <- function(term)
47 :     {
48 :     if (is.name(term) || is.numeric(term)) return(term)
49 :     if (length(term) == 2) {
50 :     term[[2]] <- subbars(term[[2]])
51 :     return(term)
52 :     }
53 :     stopifnot(length(term) == 3)
54 : maechler 832 if (is.call(term) && term[[1]] == as.name('|'))
55 :     term[[1]] <- as.name('+')
56 : bates 769 term[[2]] <- subbars(term[[2]])
57 :     term[[3]] <- subbars(term[[3]])
58 :     term
59 :     }
60 : bates 824
61 : bates 979 ## Expand an expression with colons to the sum of the lhs
62 :     ## and the current expression
63 :    
64 :     colExpand <- function(term)
65 :     {
66 :     if (is.name(term) || is.numeric(term)) return(term)
67 :     if (length(term) == 2) {
68 :     term[[2]] <- colExpand(term[[2]])
69 :     return(term)
70 :     }
71 :     stopifnot(length(term) == 3)
72 :     if (is.call(term) && term[[1]] == as.name(':')) {
73 :     return(substitute(A+B, list(A = term, B = colExpand(term[[2]]))))
74 :     }
75 :     term[[2]] <- colExpand(term[[2]])
76 :     term[[3]] <- colExpand(term[[3]])
77 :     term
78 :     }
79 :    
80 : bates 824 abbrvNms <- function(gnm, cnms)
81 :     {
82 :     ans <- paste(abbreviate(gnm), abbreviate(cnms), sep = '.')
83 :     if (length(cnms) > 1) {
84 :     anms <- lapply(cnms, abbreviate, minlength = 3)
85 :     nmmat <- outer(anms, anms, paste, sep = '.')
86 :     ans <- c(ans, paste(abbreviate(gnm, minlength = 3),
87 :     nmmat[upper.tri(nmmat)], sep = '.'))
88 :     }
89 :     ans
90 :     }
91 :    
92 : bates 775 ## Control parameters for lmer
93 :     lmerControl <-
94 : maechler 832 function(maxIter = 200, # used in ../src/lmer.c only
95 : bates 888 tolerance = sqrt(.Machine$double.eps),# ditto
96 : bates 769 msMaxIter = 200,
97 : maechler 832 ## msTol = sqrt(.Machine$double.eps),
98 :     ## FIXME: should be able to pass tolerances to nlminb()
99 :     msVerbose = getOption("verbose"),
100 : bates 752 niterEM = 15,
101 : bates 435 EMverbose = getOption("verbose"),
102 : maechler 843 PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead
103 : bates 435 analyticGradient = TRUE,
104 : maechler 832 analyticHessian = FALSE # unused _FIXME_
105 :     )
106 : bates 435 {
107 : bates 775 list(maxIter = as.integer(maxIter),
108 : maechler 832 tolerance = as.double(tolerance),
109 : bates 775 msMaxIter = as.integer(msMaxIter),
110 : maechler 832 ## msTol = as.double(msTol),
111 :     msVerbose = as.integer(msVerbose),# "integer" on purpose
112 : bates 775 niterEM = as.integer(niterEM),
113 : maechler 832 EMverbose = as.logical(EMverbose),
114 : bates 775 PQLmaxIt = as.integer(PQLmaxIt),
115 :     analyticGradient = as.logical(analyticGradient),
116 :     analyticHessian = as.logical(analyticHessian))
117 : bates 435 }
118 :    
119 : bates 755 setMethod("lmer", signature(formula = "formula"),
120 : bates 689 function(formula, data, family,
121 :     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
122 : bates 901 control = list(), start,
123 : bates 435 subset, weights, na.action, offset,
124 : maechler 832 model = TRUE, x = FALSE, y = FALSE , ...)
125 :     ## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
126 :     {
127 :     ## match and check parameters
128 : bates 755 if (length(formula) < 3) stop("formula must be a two-sided formula")
129 :     cv <- do.call("lmerControl", control)
130 :     ## evaluate glm.fit, a generalized linear fit of fixed effects only
131 : maechler 832 mf <- match.call()
132 : bates 755 m <- match(c("family", "data", "subset", "weights",
133 :     "na.action", "offset"), names(mf), 0)
134 :     mf <- mf[c(1, m)]
135 :     frame.form <- subbars(formula) # substitute `+' for `|'
136 :     fixed.form <- nobars(formula) # remove any terms with `|'
137 : bates 767 if (inherits(fixed.form, "name")) # RHS is empty - use a constant
138 : bates 755 fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
139 :     environment(fixed.form) <- environment(frame.form) <- environment(formula)
140 :     mf$formula <- fixed.form
141 :     mf$x <- mf$model <- mf$y <- TRUE
142 :     mf[[1]] <- as.name("glm")
143 :     glm.fit <- eval(mf, parent.frame())
144 : bates 767 x <- glm.fit$x
145 :     y <- as.double(glm.fit$y)
146 : bates 769 family <- glm.fit$family
147 : bates 939 ## check for a linear mixed model
148 :     lmm <- family$family == "gaussian" && family$link == "identity"
149 : maechler 832 if (lmm) { # linear mixed model
150 :     method <- match.arg(method)
151 :     if (method %in% c("PQL", "Laplace", "AGQ")) {
152 :     warning(paste('Argument method = "', method,
153 :     '" is not meaningful for a linear mixed model.\n',
154 :     'Using method = "REML".\n', sep = ''))
155 :     method <- "REML"
156 :     }
157 :     } else { # generalized linear mixed model
158 :     if (missing(method)) method <- "PQL"
159 :     else {
160 :     method <- match.arg(method)
161 :     if (method == "ML") method <- "PQL"
162 :     if (method == "REML")
163 :     warning('Argument method = "REML" is not meaningful ',
164 :     'for a generalized linear mixed model.',
165 :     '\nUsing method = "PQL".\n')
166 :     }
167 :     }
168 :    
169 : bates 755 ## evaluate a model frame for fixed and random effects
170 : bates 435 mf$formula <- frame.form
171 : bates 755 mf$x <- mf$model <- mf$y <- mf$family <- NULL
172 : bates 435 mf$drop.unused.levels <- TRUE
173 : bates 755 mf[[1]] <- as.name("model.frame")
174 : bates 435 frm <- eval(mf, parent.frame())
175 : bates 755
176 : bates 435 ## grouping factors and model matrices for random effects
177 :     bars <- findbars(formula[[3]])
178 : maechler 832 random <-
179 :     lapply(bars, function(x)
180 :     list(model.matrix(eval(substitute(~ T, list(T = x[[2]]))),
181 :     frm),
182 :     eval(substitute(as.factor(fac)[,drop = TRUE],
183 :     list(fac = x[[3]])), frm)))
184 : bates 435 names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
185 : bates 755
186 : bates 435 ## order factor list by decreasing number of levels
187 : bates 449 nlev <- sapply(random, function(x) length(levels(x[[2]])))
188 : bates 452 if (any(diff(nlev) > 0)) {
189 : bates 449 random <- random[rev(order(nlev))]
190 : bates 435 }
191 : bates 767
192 :     ## Create the model matrices and a mixed-effects representation (mer)
193 : bates 435 mmats <- c(lapply(random, "[[", 1),
194 : bates 755 .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
195 :     mer <- .Call("lmer_create", lapply(random, "[[", 2),
196 :     mmats, method, PACKAGE = "Matrix")
197 : bates 767 if (lmm) { ## linear mixed model
198 : bates 901 if (missing(start)) .Call("lmer_initial", mer, PACKAGE="Matrix")
199 : maechler 973 else .Call("lmer_set_initial", mer, start, PACKAGE = "Matrix")
200 : bates 755 .Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
201 :     LMEoptimize(mer) <- cv
202 :     fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")
203 : bates 767 return(new("lmer",
204 : bates 769 mer,
205 : bates 767 assign = attr(x, "assign"),
206 :     call = match.call(),
207 : bates 769 family = family, fitted = fits,
208 :     fixed = fixef(mer),
209 :     frame = if (model) frm else data.frame(),
210 :     logLik = logLik(mer),
211 : bates 755 residuals = unname(model.response(frm) - fits),
212 : bates 769 terms = glm.fit$terms))
213 : bates 755 }
214 :    
215 :     ## The rest of the function applies to generalized linear mixed models
216 :     gVerb <- getOption("verbose")
217 : bates 776 eta <- glm.fit$linear.predictors
218 : bates 767 wts <- glm.fit$prior.weights
219 : bates 774 wtssqr <- wts * wts
220 : bates 767 offset <- glm.fit$offset
221 :     if (is.null(offset)) offset <- numeric(length(eta))
222 : bates 776 mu <- numeric(length(eta))
223 : bates 767
224 : bates 774 dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
225 : bates 767 linkinv <- quote(family$linkinv(eta))
226 :     mu.eta <- quote(family$mu.eta(eta))
227 :     variance <- quote(family$variance(mu))
228 : bates 775 LMEopt <- get("LMEoptimize<-")
229 :     doLMEopt <- quote(LMEopt(x = mer, value = cv))
230 : bates 767
231 : bates 809 GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
232 :     .Call("glmer_PQL", GSpt, PACKAGE = "Matrix") # obtain PQL estimates
233 : bates 755
234 : bates 774 fixInd <- seq(ncol(x))
235 :     ## pars[fixInd] == beta, pars[-fixInd] == theta
236 :     PQLpars <- c(fixef(mer),
237 :     .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))
238 : bates 775 ## set flag to skip fixed-effects in subsequent calls
239 :     mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
240 : bates 777 ## indicator of constrained parameters
241 :     const <- c(rep(FALSE, length(fixInd)),
242 :     unlist(lapply(mer@nc[seq(along = random)],
243 :     function(k) 1:((k*(k+1))/2) <= k)
244 :     ))
245 : bates 779 devAGQ <- function(pars, n)
246 :     .Call("glmer_devAGQ", pars, GSpt, n, PACKAGE = "Matrix")
247 : maechler 832
248 : bates 801 deviance <- devAGQ(PQLpars, 1)
249 : bates 804 ### FIXME: For nf == 1 change this to an AGQ evaluation. Needs
250 : bates 801 ### AGQ for nc > 1 first.
251 : bates 777 fxd <- PQLpars[fixInd]
252 : bates 779 loglik <- logLik(mer)
253 : bates 775
254 : bates 777 if (method %in% c("Laplace", "AGQ")) {
255 : bates 779 nAGQ <- 1
256 :     if (method == "AGQ") { # determine nAGQ at PQL estimates
257 :     dev11 <- devAGQ(PQLpars, 11)
258 : bates 799 ## FIXME: Should this be an absolute or a relative tolerance?
259 : bates 779 devTol <- sqrt(.Machine$double.eps) * abs(dev11)
260 : bates 799 for (nAGQ in c(9, 7, 5, 3, 1))
261 : bates 779 if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break
262 : bates 799 nAGQ <- nAGQ + 2
263 :     if (gVerb)
264 :     cat(paste("Using", nAGQ, "quadrature points per column\n"))
265 : bates 779 }
266 :     obj <- function(pars)
267 :     .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
268 : bates 777 if (exists("nlminb", mode = "function")) {
269 : bates 755 optimRes <-
270 : bates 779 nlminb(PQLpars, obj,
271 : bates 755 lower = ifelse(const, 5e-10, -Inf),
272 :     control = list(trace = getOption("verbose"),
273 : maechler 832 iter.max = cv$msMaxIter))
274 : bates 755 optpars <- optimRes$par
275 :     if (optimRes$convergence != 0)
276 :     warning("nlminb failed to converge")
277 : bates 779 deviance <- optimRes$objective
278 : bates 755 } else {
279 :     optimRes <-
280 : bates 779 optim(PQLpars, obj, method = "L-BFGS-B",
281 : bates 755 lower = ifelse(const, 5e-10, -Inf),
282 :     control = list(trace = getOption("verbose"),
283 : bates 776 maxit = cv$msMaxIter))
284 : bates 755 optpars <- optimRes$par
285 :     if (optimRes$convergence != 0)
286 :     warning("optim failed to converge")
287 : bates 779 deviance <- optimRes$value
288 : bates 755 }
289 : bates 774 if (gVerb) {
290 : bates 772 cat(paste("convergence message", optimRes$message, "\n"))
291 : bates 777 }
292 :     fxd[] <- optpars[fixInd] ## preserve the names
293 : bates 809 .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
294 : bates 755 }
295 :    
296 : bates 776 .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
297 : bates 779 loglik[] <- -deviance/2
298 : maechler 832 new("lmer", mer,
299 :     frame = if (model) frm else data.frame(),
300 :     terms = glm.fit$terms,
301 : bates 777 assign = attr(glm.fit$x, "assign"),
302 :     call = match.call(), family = family,
303 :     logLik = loglik, fixed = fxd)
304 : bates 435 })
305 : maechler 832 ## end{ "lmer . formula " }
306 : bates 435
307 : bates 755 setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
308 : bates 316 function(x, value)
309 :     {
310 :     if (value$msMaxIter < 1) return(x)
311 :     nc <- x@nc
312 : bates 755 constr <- unlist(lapply(nc[1:(length(nc) - 2)],
313 :     function(k) 1:((k*(k+1))/2) <= k))
314 : bates 752 fn <- function(pars)
315 : bates 755 deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
316 : maechler 832 gr <-
317 :     if (value$analyticGradient)
318 : bates 755 function(pars) {
319 :     if (!isTRUE(all.equal(pars,
320 :     .Call("lmer_coef", x,
321 :     2, PACKAGE = "Matrix"))))
322 :     .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
323 :     .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
324 :     }
325 : maechler 832 ## else NULL
326 :     optimRes <-
327 :     if (exists("nlminb", mode = "function"))
328 :     nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
329 :     fn, gr,
330 :     lower = ifelse(constr, 5e-10, -Inf),
331 :     control = list(iter.max = value$msMaxIter,
332 :     trace = as.integer(value$msVerbose)))
333 :     else
334 :     optim(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
335 :     fn, gr, method = "L-BFGS-B",
336 :     lower = ifelse(constr, 5e-10, -Inf),
337 :     control = list(maxit = value$msMaxIter,
338 :     trace = as.integer(value$msVerbose)))
339 : bates 755 .Call("lmer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
340 : bates 316 if (optimRes$convergence != 0) {
341 : bates 777 warning(paste("optim or nlminb returned message",
342 :     optimRes$message,"\n"))
343 : bates 316 }
344 :     return(x)
345 :     })
346 :    
347 : bates 413 setMethod("ranef", signature(object = "lmer"),
348 : bates 689 function(object, accumulate = FALSE, ...) {
349 :     val <- new("lmer.ranef",
350 :     lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),
351 :     data.frame, check.names = FALSE),
352 :     varFac = object@bVar,
353 :     stdErr = .Call("lmer_sigma", object,
354 : bates 755 object@method == "REML", PACKAGE = "Matrix"))
355 : bates 689 if (!accumulate || length(val@varFac) == 1) return(val)
356 :     ## check for nested factors
357 :     L <- object@L
358 :     if (any(sapply(seq(a = val), function(i) length(L[[Lind(i,i)]]@i))))
359 :     error("Require nested grouping factors to accumulate random effects")
360 :     val
361 : bates 316 })
362 :    
363 : bates 755 setMethod("fixef", signature(object = "mer"),
364 : bates 774 function(object, ...)
365 :     .Call("lmer_fixef", object, PACKAGE = "Matrix"))
366 : bates 316
367 : bates 769 setMethod("fixef", signature(object = "lmer"),
368 :     function(object, ...) object@fixed)
369 : deepayan 721
370 : bates 413 setMethod("VarCorr", signature(x = "lmer"),
371 : bates 902 ##FIXME - change this for reasonable defaults of useScale according to
372 :     ##the family slot.
373 : bates 316 function(x, REML = TRUE, useScale = TRUE, ...) {
374 : bates 550 val <- .Call("lmer_variances", x, PACKAGE = "Matrix")
375 : bates 316 for (i in seq(along = val)) {
376 :     dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])
377 :     val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
378 :     }
379 :     new("VarCorr",
380 : bates 449 scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),
381 : bates 316 reSumry = val,
382 :     useScale = useScale)
383 :     })
384 :    
385 : bates 413 setMethod("gradient", signature(x = "lmer"),
386 : bates 755 function(x, unconst, ...)
387 :     .Call("lmer_gradient", x, unconst, PACKAGE = "Matrix"))
388 : bates 316
389 : bates 449 setMethod("summary", signature(object = "lmer"),
390 :     function(object, ...)
391 : bates 769 new("summary.lmer", object,
392 : bates 727 showCorrelation = TRUE,
393 : bates 769 useScale = !((object@family)$family %in% c("binomial", "poisson"))))
394 : bates 316
395 : bates 449 setMethod("show", signature(object = "lmer"),
396 :     function(object)
397 : bates 769 show(new("summary.lmer", object,
398 : bates 727 showCorrelation = FALSE,
399 : bates 769 useScale = !((object@family)$family %in% c("binomial", "poisson")))))
400 : maechler 832
401 : bates 449 setMethod("show", "summary.lmer",
402 : bates 316 function(object) {
403 : bates 727 fcoef <- object@fixed
404 : bates 449 useScale <- object@useScale
405 :     corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"),
406 : bates 316 "corrmatrix")
407 :     DF <- getFixDF(object)
408 :     coefs <- cbind(fcoef, corF@stdDev, DF)
409 :     nc <- object@nc
410 :     dimnames(coefs) <-
411 :     list(names(fcoef), c("Estimate", "Std. Error", "DF"))
412 : bates 449 digits <- max(3, getOption("digits") - 2)
413 : bates 755 REML <- object@method == "REML"
414 : bates 727 llik <- object@logLik
415 : bates 449 dev <- object@deviance
416 : maechler 832
417 : bates 449 rdig <- 5
418 : bates 727 if (glz <- !(object@method %in% c("REML", "ML"))) {
419 :     cat(paste("Generalized linear mixed model fit using",
420 :     object@method, "\n"))
421 :     } else {
422 :     cat("Linear mixed-effects model fit by ")
423 : bates 755 cat(if(REML) "REML\n" else "maximum likelihood\n")
424 : bates 727 }
425 : bates 449 if (!is.null(object@call$formula)) {
426 :     cat("Formula:", deparse(object@call$formula),"\n")
427 :     }
428 :     if (!is.null(object@call$data)) {
429 :     cat(" Data:", deparse(object@call$data), "\n")
430 :     }
431 :     if (!is.null(object@call$subset)) {
432 :     cat(" Subset:",
433 :     deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
434 :     }
435 : bates 727 if (glz) {
436 : bates 750 cat(" Family: ", object@family$family, "(",
437 :     object@family$link, " link)\n", sep = "")
438 : bates 727 print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
439 : bates 449 logLik = c(llik),
440 : bates 727 deviance = -2*llik,
441 :     row.names = ""))
442 :     } else {
443 :     print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
444 :     logLik = c(llik),
445 : bates 750 MLdeviance = dev["ML"],
446 : bates 449 REMLdeviance = dev["REML"],
447 :     row.names = ""))
448 : bates 727 }
449 : bates 449 cat("Random effects:\n")
450 : bates 777 show(VarCorr(object, useScale = useScale))
451 : bates 449 ngrps <- lapply(object@flist, function(x) length(levels(x)))
452 :     cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))
453 :     cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
454 :     cat("\n")
455 :     if (!useScale)
456 :     cat("\nEstimated scale (compare to 1) ",
457 : bates 755 .Call("lmer_sigma", object, FALSE, PACKAGE = "Matrix"),
458 : bates 449 "\n")
459 :     if (nrow(coefs) > 0) {
460 :     if (useScale) {
461 :     stat <- coefs[,1]/coefs[,2]
462 :     pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
463 :     nms <- colnames(coefs)
464 :     coefs <- cbind(coefs, stat, pval)
465 :     colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")
466 :     } else {
467 :     coefs <- coefs[, 1:2, drop = FALSE]
468 :     stat <- coefs[,1]/coefs[,2]
469 :     pval <- 2*pnorm(abs(stat), lower = FALSE)
470 :     nms <- colnames(coefs)
471 :     coefs <- cbind(coefs, stat, pval)
472 :     colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
473 :     }
474 :     cat("\nFixed effects:\n")
475 :     printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
476 :     if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
477 :     rn <- rownames(coefs)
478 :     dimnames(corF) <- list(
479 :     abbreviate(rn, minlen=11),
480 :     abbreviate(rn, minlen=6))
481 :     if (!is.null(corF)) {
482 :     p <- NCOL(corF)
483 :     if (p > 1) {
484 :     cat("\nCorrelation of Fixed Effects:\n")
485 :     corF <- format(round(corF, 3), nsmall = 3)
486 :     corF[!lower.tri(corF)] <- ""
487 :     print(corF[-1, -p, drop=FALSE], quote = FALSE)
488 :     }
489 :     }
490 :     }
491 :     }
492 :     invisible(object)
493 : bates 316 })
494 :    
495 :     ## calculates degrees of freedom for fixed effects Wald tests
496 :     ## This is a placeholder. The answers are generally wrong. It will
497 :     ## be very tricky to decide what a 'right' answer should be with
498 :     ## crossed random effects.
499 :    
500 : bates 413 setMethod("getFixDF", signature(object="lmer"),
501 : bates 316 function(object, ...)
502 :     {
503 :     nc <- object@nc[-seq(along = object@Omega)]
504 : bates 777 p <- abs(nc[1]) - 1
505 : bates 316 n <- nc[2]
506 :     rep(n-p, p)
507 :     })
508 :    
509 : bates 755 setMethod("logLik", signature(object="mer"),
510 :     function(object, REML = object@method == "REML", ...) {
511 : bates 446 val <- -deviance(object, REML = REML)/2
512 :     nc <- object@nc[-seq(a = object@Omega)]
513 :     attr(val, "nall") <- attr(val, "nobs") <- nc[2]
514 : bates 782 attr(val, "df") <- abs(nc[1]) +
515 : bates 755 length(.Call("lmer_coef", object, 0, PACKAGE = "Matrix"))
516 : maechler 832 attr(val, "REML") <- REML
517 : bates 446 class(val) <- "logLik"
518 :     val
519 :     })
520 :    
521 : bates 769 setMethod("logLik", signature(object="lmer"),
522 :     function(object, ...) object@logLik)
523 : deepayan 721
524 : bates 446 setMethod("anova", signature(object = "lmer"),
525 :     function(object, ...)
526 :     {
527 :     mCall <- match.call(expand.dots = TRUE)
528 :     dots <- list(...)
529 :     modp <- logical(0)
530 :     if (length(dots))
531 :     modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm")
532 :     if (any(modp)) { # multiple models - form table
533 :     opts <- dots[!modp]
534 :     mods <- c(list(object), dots[modp])
535 :     names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)
536 :     mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE), attr, "df"))]
537 :     calls <- lapply(mods, slot, "call")
538 :     data <- lapply(calls, "[[", "data")
539 :     if (any(data != data[[1]])) stop("all models must be fit to the same data object")
540 :     header <- paste("Data:", data[[1]])
541 :     subset <- lapply(calls, "[[", "subset")
542 :     if (any(subset != subset[[1]])) stop("all models must use the same subset")
543 :     if (!is.null(subset[[1]]))
544 :     header <-
545 :     c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
546 :     llks <- lapply(mods, logLik, REML = FALSE)
547 :     Df <- sapply(llks, attr, "df")
548 :     llk <- unlist(llks)
549 :     chisq <- 2 * pmax(0, c(NA, diff(llk)))
550 :     dfChisq <- c(NA, diff(Df))
551 :     val <- data.frame(Df = Df,
552 :     AIC = sapply(llks, AIC),
553 :     BIC = sapply(llks, BIC),
554 :     logLik = llk,
555 :     "Chisq" = chisq,
556 :     "Chi Df" = dfChisq,
557 :     "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
558 :     check.names = FALSE)
559 :     class(val) <- c("anova", class(val))
560 :     attr(val, "heading") <-
561 : bates 690 c(header, "Models:",
562 : bates 446 paste(names(mods),
563 :     unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
564 : bates 690 sep = ": "))
565 : bates 446 return(val)
566 :     } else {
567 : bates 571 foo <- object
568 :     foo@status["factored"] <- FALSE
569 :     .Call("lmer_factor", foo, PACKAGE="Matrix")
570 :     dfr <- getFixDF(foo)
571 :     rcol <- ncol(foo@RXX)
572 :     ss <- foo@RXX[ , rcol]^2
573 :     ssr <- ss[[rcol]]
574 :     ss <- ss[seq(along = dfr)]
575 :     names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
576 :     asgn <- foo@assign
577 :     terms <- foo@terms
578 :     nmeffects <- attr(terms, "term.labels")
579 :     if ("(Intercept)" %in% names(ss))
580 :     nmeffects <- c("(Intercept)", nmeffects)
581 :     ss <- unlist(lapply(split(ss, asgn), sum))
582 :     df <- unlist(lapply(split(asgn, asgn), length))
583 :     dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
584 :     ms <- ss/df
585 :     f <- ms/(ssr/dfr)
586 :     P <- pf(f, df, dfr, lower.tail = FALSE)
587 :     table <- data.frame(df, ss, ms, dfr, f, P)
588 :     dimnames(table) <-
589 :     list(nmeffects,
590 :     c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
591 :     if ("(Intercept)" %in% nmeffects) table <- table[-1,]
592 :     attr(table, "heading") <- "Analysis of Variance Table"
593 :     class(table) <- c("anova", "data.frame")
594 :     table
595 : bates 446 }
596 : bates 316 })
597 : bates 446
598 :     setMethod("update", signature(object = "lmer"),
599 :     function(object, formula., ..., evaluate = TRUE)
600 :     {
601 :     call <- object@call
602 :     if (is.null(call))
603 :     stop("need an object with call component")
604 :     extras <- match.call(expand.dots = FALSE)$...
605 :     if (!missing(formula.))
606 :     call$formula <- update.formula(formula(object), formula.)
607 :     if (length(extras) > 0) {
608 :     existing <- !is.na(match(names(extras), names(call)))
609 :     for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
610 :     if (any(!existing)) {
611 :     call <- c(as.list(call), extras[!existing])
612 :     call <- as.call(call)
613 :     }
614 :     }
615 :     if (evaluate)
616 :     eval(call, parent.frame())
617 :     else call
618 :     })
619 :    
620 :    
621 :     setMethod("confint", signature(object = "lmer"),
622 : maechler 832 function (object, parm, level = 0.95, ...)
623 : bates 446 {
624 :     cf <- fixef(object)
625 :     pnames <- names(cf)
626 : maechler 832 if (missing(parm))
627 : bates 446 parm <- seq(along = pnames)
628 : maechler 832 else if (is.character(parm))
629 : bates 446 parm <- match(parm, pnames, nomatch = 0)
630 :     a <- (1 - level)/2
631 :     a <- c(a, 1 - a)
632 :     pct <- paste(round(100 * a, 1), "%")
633 :     ci <- array(NA, dim = c(length(parm), 2),
634 :     dimnames = list(pnames[parm], pct))
635 :     ses <- sqrt(diag(vcov(object)))[parm]
636 : bates 449 ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
637 : bates 446 ci
638 :     })
639 :    
640 : bates 755 setMethod("deviance", "mer",
641 : bates 449 function(object, REML = NULL, ...) {
642 :     .Call("lmer_factor", object, PACKAGE = "Matrix")
643 :     if (is.null(REML))
644 : bates 755 REML <- object@method == "REML"
645 : bates 449 object@deviance[[ifelse(REML, "REML", "ML")]]
646 :     })
647 : bates 446
648 : deepayan 721
649 : bates 769 setMethod("deviance", "lmer",
650 :     function(object, ...) -2 * c(object@logLik))
651 : deepayan 721
652 : bates 769
653 : bates 449 setMethod("chol", signature(x = "lmer"),
654 :     function(x, pivot = FALSE, LINPACK = pivot) {
655 :     x@status["factored"] <- FALSE # force a decomposition
656 :     .Call("lmer_factor", x, PACKAGE = "Matrix")
657 :     })
658 :    
659 :     setMethod("solve", signature(a = "lmer", b = "missing"),
660 :     function(a, b, ...)
661 : bates 562 .Call("lmer_invert", a, PACKAGE = "Matrix")
662 : bates 449 )
663 :    
664 :     setMethod("formula", "lmer", function(x, ...) x@call$formula)
665 :    
666 :     setMethod("vcov", signature(object = "lmer"),
667 : bates 755 function(object, REML = object@method == "REML", useScale = TRUE,...) {
668 : bates 449 sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
669 :     rr <- object@RXX
670 :     nms <- object@cnames[[".fixed"]]
671 :     dimnames(rr) <- list(nms, nms)
672 :     nr <- nrow(rr)
673 :     rr <- rr[-nr, -nr, drop = FALSE]
674 :     rr <- rr %*% t(rr)
675 :     if (useScale) {
676 :     rr = sc^2 * rr
677 :     }
678 :     rr
679 :     })
680 :    
681 : maechler 832 ## Extract the L matrix
682 : bates 550 setAs("lmer", "dtTMatrix",
683 :     function(from)
684 :     {
685 :     ## force a refactorization if the factors have been inverted
686 :     if (from@status["inverted"]) from@status["factored"] <- FALSE
687 :     .Call("lmer_factor", from, PACKAGE = "Matrix")
688 :     L <- lapply(from@L, as, "dgTMatrix")
689 :     nf <- length(from@D)
690 :     Gp <- from@Gp
691 :     nL <- Gp[nf + 1]
692 : bates 562 Li <- integer(0)
693 :     Lj <- integer(0)
694 :     Lx <- double(0)
695 : bates 550 for (i in 1:nf) {
696 :     for (j in 1:i) {
697 :     Lij <- L[[Lind(i, j)]]
698 : bates 562 Li <- c(Li, Lij@i + Gp[i])
699 :     Lj <- c(Lj, Lij@j + Gp[j])
700 :     Lx <- c(Lx, Lij@x)
701 : bates 550 }
702 :     }
703 : bates 562 new("dtTMatrix", Dim = as.integer(c(nL, nL)), i = Li, j = Lj, x = Lx,
704 : bates 550 uplo = "L", diag = "U")
705 :     })
706 : bates 562
707 :     ## Extract the ZZX matrix
708 :     setAs("lmer", "dsTMatrix",
709 :     function(from)
710 :     {
711 :     .Call("lmer_inflate", from, PACKAGE = "Matrix")
712 :     ZZpO <- lapply(from@ZZpO, as, "dgTMatrix")
713 :     ZZ <- lapply(from@ZtZ, as, "dgTMatrix")
714 :     nf <- length(ZZpO)
715 :     Gp <- from@Gp
716 :     nZ <- Gp[nf + 1]
717 :     Zi <- integer(0)
718 :     Zj <- integer(0)
719 :     Zx <- double(0)
720 :     for (i in 1:nf) {
721 :     ZZpOi <- ZZpO[[i]]
722 :     Zi <- c(Zi, ZZpOi@i + Gp[i])
723 :     Zj <- c(Zj, ZZpOi@j + Gp[i])
724 :     Zx <- c(Zx, ZZpOi@x)
725 :     if (i > 1) {
726 :     for (j in 1:(i-1)) {
727 :     ZZij <- ZZ[[Lind(i, j)]]
728 :     ## off-diagonal blocks are transposed
729 :     Zi <- c(Zi, ZZij@j + Gp[j])
730 :     Zj <- c(Zj, ZZij@i + Gp[i])
731 :     Zx <- c(Zx, ZZij@x)
732 :     }
733 :     }
734 :     }
735 :     new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx,
736 :     uplo = "U")
737 :     })
738 : bates 689
739 :     setMethod("fitted", signature(object = "lmer"),
740 : bates 691 function(object, ...)
741 :     napredict(attr(object@frame, "na.action"), object@fitted))
742 : bates 689
743 :     setMethod("residuals", signature(object = "lmer"),
744 : bates 691 function(object, ...)
745 :     naresid(attr(object@frame, "na.action"), object@residuals))
746 : bates 689
747 :     setMethod("resid", signature(object = "lmer"),
748 :     function(object, ...) do.call("residuals", c(list(object), list(...))))
749 :    
750 :     setMethod("coef", signature(object = "lmer"),
751 :     function(object, ...)
752 :     {
753 : bates 769 fef <- data.frame(rbind(object@fixed), check.names = FALSE)
754 : bates 689 ref <- as(ranef(object), "list")
755 :     names(ref) <- names(object@flist)
756 :     val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
757 :     for (i in seq(a = val)) {
758 :     refi <- ref[[i]]
759 :     row.names(val[[i]]) <- row.names(refi)
760 :     if (!all(names(refi) %in% names(fef)))
761 :     stop("unable to align random and fixed effects")
762 :     val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
763 :     }
764 :     new("lmer.coef", val)
765 :     })
766 :    
767 :     setMethod("plot", signature(x = "lmer.coef"),
768 :     function(x, y, ...)
769 :     {
770 : maechler 832 ## require("lattice", quietly = TRUE) -- now via Imports
771 :     varying <- unique(do.call("c",
772 :     lapply(x, function(el)
773 :     names(el)[sapply(el,
774 :     function(col)
775 :     any(col != col[1]))])))
776 :     gf <- do.call("rbind", lapply(x, "[", j = varying))
777 :     gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
778 :     switch(min(length(varying), 3),
779 :     qqmath(eval(substitute(~ x | .grp,
780 :     list(x = as.name(varying[1])))), gf, ...),
781 :     xyplot(eval(substitute(y ~ x | .grp,
782 :     list(y = as.name(varying[1]),
783 :     x = as.name(varying[2])))), gf, ...),
784 :     splom(~ gf | .grp, ...))
785 : bates 689 })
786 :    
787 :     setMethod("plot", signature(x = "lmer.ranef"),
788 :     function(x, y, ...)
789 :     {
790 : maechler 832 ## require("lattice", quietly = TRUE) -- now via Imports
791 :     lapply(x, function(x) {
792 :     cn <- lapply(colnames(x), as.name)
793 :     switch(min(ncol(x), 3),
794 :     qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
795 :     xyplot(eval(substitute(y ~ x,
796 :     list(y = cn[[1]],
797 :     x = cn[[2]]))), x, ...),
798 :     splom(~ x, ...))
799 :     })
800 : bates 689 })
801 :    
802 :     setMethod("with", signature(data = "lmer"),
803 : bates 690 function(data, expr, ...) {
804 : bates 691 dat <- eval(data@call$data)
805 :     if (!is.null(na.act <- attr(data@frame, "na.action")))
806 :     dat <- dat[-na.act, ]
807 :     lst <- c(list(. = data), data@flist, data@frame, dat)
808 :     eval(substitute(expr), lst[unique(names(lst))])
809 :     })
810 : bates 690
811 : bates 691 setMethod("terms", signature(x = "lmer"),
812 :     function(x, ...) x@terms)
813 : bates 767
814 :     setMethod("show", signature(object="VarCorr"),
815 :     function(object)
816 :     {
817 :     digits <- max(3, getOption("digits") - 2)
818 :     useScale <- length(object@useScale) > 0 && object@useScale[1]
819 :     sc <- ifelse(useScale, object@scale, 1.)
820 :     reStdDev <- c(lapply(object@reSumry,
821 :     function(x, sc)
822 :     sc*x@stdDev,
823 :     sc = sc), list(Residual = sc))
824 :     reLens <- unlist(c(lapply(reStdDev, length)))
825 :     reMat <- array('', c(sum(reLens), 4),
826 :     list(rep('', sum(reLens)),
827 :     c("Groups", "Name", "Variance", "Std.Dev.")))
828 :     reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
829 :     reMat[,2] <- c(unlist(lapply(reStdDev, names)), "")
830 :     reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
831 :     reMat[,4] <- format(unlist(reStdDev), digits = digits)
832 :     if (any(reLens > 1)) {
833 :     maxlen <- max(reLens)
834 :     corr <-
835 :     do.call("rbind",
836 :     lapply(object@reSumry,
837 :     function(x, maxlen) {
838 :     cc <- format(round(x, 3), nsmall = 3)
839 :     cc[!lower.tri(cc)] <- ""
840 :     nr <- dim(cc)[1]
841 :     if (nr >= maxlen) return(cc)
842 :     cbind(cc, matrix("", nr, maxlen-nr))
843 :     }, maxlen))
844 :     colnames(corr) <- c("Corr", rep("", maxlen - 1))
845 :     reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr))))
846 :     }
847 :     if (!useScale) reMat <- reMat[-nrow(reMat),]
848 :     print(reMat, quote = FALSE)
849 :     })
850 : bates 769
851 : bates 879 setMethod("mcmcsamp", signature(object = "lmer"),
852 :     function(object, n = 1, verbose = FALSE, saveb = FALSE,
853 : bates 824 trans = TRUE, ...)
854 : bates 820 {
855 : bates 879 if (object@family$family == "gaussian" &&
856 :     object@family$link == "identity") {
857 : bates 861 glmer <- FALSE
858 : bates 879 ans <- .Call("lmer_MCMCsamp", object, saveb, n, trans,
859 : bates 820 PACKAGE = "Matrix")
860 :     } else {
861 : bates 861 glmer <- TRUE
862 :     if (trans)
863 : bates 864 warning("trans option not currently allowed for generalized models")
864 : bates 861 trans <- FALSE
865 : bates 820 ## Check arguments
866 : bates 879 if (length(object@Omega) > 1 || object@nc[1] > 1)
867 : bates 820 stop("mcmcsamp currently defined for glmm models with only one variance component")
868 :     cv <- Matrix:::lmerControl()
869 :     if (verbose) cv$msVerbose <- 1
870 : bates 879 family <- object@family
871 :     frm <- object@frame
872 : bates 810
873 : bates 820 ## recreate model matrices
874 : bates 879 fixed.form <- Matrix:::nobars(object@call$formula)
875 : bates 820 if (inherits(fixed.form, "name")) # RHS is empty - use a constant
876 :     fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
877 :     glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE,
878 :     y = TRUE)
879 :     x <- glm.fit$x
880 :     y <- as.double(glm.fit$y)
881 : bates 879 bars <- Matrix:::findbars(object@call$formula[[3]])
882 : bates 820 random <-
883 :     lapply(bars,
884 :     function(x) list(model.matrix(eval(substitute(~term,
885 :     list(term=x[[2]]))),
886 :     frm),
887 :     eval(substitute(as.factor(fac)[,drop = TRUE],
888 :     list(fac = x[[3]])), frm)))
889 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
890 : bates 879 if (any(names(random) != names(object@flist)))
891 :     random <- random[names(object@flist)]
892 : bates 820 mmats <- c(lapply(random, "[[", 1),
893 :     .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
894 : bates 879 mer <- as(object, "mer")
895 : bates 781
896 : bates 820 ## establish the GS object and the ans matrix
897 : bates 879 eta <- glm.fit$linear.predictors # perhaps later change this to object@fitted?
898 : bates 820 wts <- glm.fit$prior.weights
899 :     wtssqr <- wts * wts
900 :     offset <- glm.fit$offset
901 :     if (is.null(offset)) offset <- numeric(length(eta))
902 :     off <- numeric(length(eta))
903 :     mu <- numeric(length(eta))
904 :     dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
905 :     linkinv <- quote(family$linkinv(eta))
906 :     mu.eta <- quote(family$mu.eta(eta))
907 :     variance <- quote(family$variance(mu))
908 :     LMEopt <- getAnywhere("LMEoptimize<-")
909 :     doLMEopt <- quote(LMEopt(x = mer, value = cv))
910 :     GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
911 : bates 879 fixed <- object@fixed
912 : bates 820 varc <- .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")
913 :     b <- .Call("lmer_ranef", mer, PACKAGE = "Matrix")
914 : bates 879 ans <- .Call("glmer_MCMCsamp", GSpt, b, fixed, varc, saveb, n,
915 : maechler 832 PACKAGE = "Matrix")
916 : bates 820 .Call("glmer_finalize", GSpt, PACKAGE = "Matrix");
917 :     }
918 : bates 879 gnms <- names(object@flist)
919 :     cnms <- object@cnames
920 :     ff <- fixef(object)
921 : bates 861 colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
922 : bates 842 unlist(lapply(seq(along = gnms),
923 :     function(i)
924 :     abbrvNms(gnms[i],cnms[[i]]))))
925 :     if (trans) {
926 :     ## parameter type: 0 => fixed effect, 1 => variance,
927 :     ## 2 => covariance
928 : bates 861 ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1,
929 : bates 842 unlist(lapply(seq(along = gnms),
930 :     function(i)
931 :     {
932 :     k <- length(cnms[[i]])
933 :     rep(1:2, c(k, (k*(k-1))/2))
934 :     })))
935 :     colnms[ptyp == 1] <-
936 :     paste("log(", colnms[ptyp == 1], ")", sep = "")
937 :     colnms[ptyp == 2] <-
938 :     paste("atanh(", colnms[ptyp == 2], ")", sep = "")
939 :     }
940 :     colnames(ans) <- colnms
941 : bates 820 ans
942 :     })
943 : bates 781
944 : bates 812 rWishart <- function(n, df, invScal)
945 :     .Call("Matrix_rWishart", n, df, invScal)
946 : bates 820
947 : bates 888
948 :     setMethod("model.matrix", signature(object = "lmer"),
949 :     function(object, ...)
950 : bates 878 {
951 : bates 879 frm <- object@frame
952 :     fixed.form <- Matrix:::nobars(object@call$formula)
953 : bates 878 if (inherits(fixed.form, "name")) # RHS is empty - use a constant
954 :     fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
955 : bates 888 glm.fit <- glm(eval(fixed.form), object@family, frm, x = TRUE, y = TRUE)
956 : bates 879 fxd <- unname(drop(glm.fit$x %*% fixef(object)))
957 : bates 888
958 :     ## Create the random effects model matrices
959 : bates 879 bars <- Matrix:::findbars(object@call$formula[[3]])
960 : bates 878 random <-
961 :     lapply(bars,
962 :     function(x) list(model.matrix(eval(substitute(~term,
963 :     list(term=x[[2]]))),
964 :     frm),
965 :     eval(substitute(as.factor(fac)[,drop = TRUE],
966 :     list(fac = x[[3]])), frm)))
967 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
968 :     ## re-order the random effects pairs if necessary
969 : bates 879 if (any(names(random) != names(object@flist)))
970 :     random <- random[names(object@flist)]
971 : bates 888 c(lapply(random, "[[", 1),
972 :     .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
973 :     })
974 :    
975 :     setMethod("simulate", signature(object = "lmer"),
976 : maechler 973 function(object, nsim = 1, seed = NULL, ...)
977 : bates 888 {
978 : maechler 973 if(!exists(".Random.seed", envir = .GlobalEnv))
979 :     runif(1) # initialize the RNG if necessary
980 :     if(is.null(seed))
981 :     RNGstate <- .Random.seed
982 :     else {
983 :     R.seed <- .Random.seed
984 :     set.seed(seed)
985 :     RNGstate <- structure(seed, kind = as.list(RNGkind()))
986 :     on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
987 :     }
988 : bates 888
989 :     family <- object@family
990 :     if (family$family != "gaussian" ||
991 : maechler 973 family$link != "identity")
992 : bates 888 stop("simulation of generalized linear mixed models not yet implemented")
993 :    
994 :     ## pieces we will need later
995 :     scale <- .Call("lmer_sigma", object, object@method == "REML",
996 :     PACKAGE = "Matrix")
997 :     mmats <- model.matrix(object)
998 :     ff <- fixef(object)
999 :    
1000 :     ###_FIXME: If the factor levels have been permuted, has the
1001 :     ### permutation been applied in the stored frame? Otherwise we
1002 :     ### need to check this.
1003 :    
1004 : bates 879 ## similate the linear predictors
1005 : bates 888 lpred <- .Call("lmer_simulate", as(object, "mer"), nsim,
1006 :     unname(drop(mmats[[length(mmats)]][,seq(a = ff)] %*% ff)),
1007 :     mmats, TRUE, PACKAGE = "Matrix")
1008 : bates 879 ## add per-observation noise term
1009 : bates 888 lpred <- as.data.frame(lpred + rnorm(prod(dim(lpred)), sd = scale))
1010 :    
1011 : maechler 973 ## save the seed
1012 :     attr(lpred, "seed") <- RNGstate
1013 : bates 888 lpred
1014 : bates 878 })
1015 : bates 879
1016 :     simulate2 <- function(object, n = 1, ...)
1017 :     {
1018 :     family <- object@family
1019 :     if (family$family != "gaussian" ||
1020 : maechler 973 family$link != "identity")
1021 : bates 879 stop("simulation of generalized linear mixed models not implemented yet")
1022 :    
1023 :     ## create the mean from the fixed effects
1024 :     frm <- object@frame
1025 :     fixed.form <- Matrix:::nobars(object@call$formula)
1026 :     if (inherits(fixed.form, "name")) # RHS is empty - use a constant
1027 :     fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
1028 :     glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE, y = TRUE)
1029 :     lpred <- matrix(glm.fit$x %*% fixef(object), nr = nrow(frm), nc = n)
1030 :    
1031 :     ## Create the random effects model matrices
1032 :     bars <- Matrix:::findbars(object@call$formula[[3]])
1033 :     random <-
1034 :     lapply(bars,
1035 :     function(x) list(model.matrix(eval(substitute(~term,
1036 :     list(term=x[[2]]))),
1037 :     frm),
1038 :     eval(substitute(as.factor(fac)[,drop = TRUE],
1039 :     list(fac = x[[3]])), frm)))
1040 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
1041 :     ## re-order the random effects pairs if necessary
1042 :     flist <- object@flist
1043 :     if (any(names(random) != names(flist)))
1044 :     random <- random[names(flist)]
1045 :     mmats <- lapply(random, "[[", 1)
1046 : maechler 973
1047 : bates 879 ## simulate the random effects
1048 :     scale <- .Call("lmer_sigma", object, object@method == "REML",
1049 :     PACKAGE = "Matrix")
1050 :     Omega <- object@Omega
1051 :     re <- lapply(seq(along = Omega),
1052 :     function(i) {
1053 :     om <- Omega[[i]]
1054 :     nr <- nrow(om)
1055 :     nlev <- length(levels(flist[[i]]))
1056 :     scale * array(solve(chol(new("dpoMatrix", Dim = dim(om),
1057 :     uplo = "U", x = c(om))),
1058 :     matrix(rnorm(nr * n * nlev),
1059 :     nr = nr))@x, c(nr, n, nlev))
1060 :     })
1061 :     ## apply the random effects
1062 :     for (j in seq(along = Omega)) {
1063 :     for (i in 1:nrow(lpred))
1064 :     lpred[i,] <- lpred[i,] + mmats[[j]][i,] %*% re[[j]][, , as.integer(flist[[j]])[i]]
1065 :     }
1066 :     ## add per-observation noise term
1067 :     lpred <- lpred + rnorm(prod(dim(lpred)), sd = scale)
1068 :     attr(lpred, "re") <- re
1069 :     lpred
1070 :     }
1071 :    
1072 :     refdist <- function(fm1, fm2, n, ...)
1073 :     {
1074 : bates 888 cv <- lmerControl()
1075 :     obs <- deviance(fm2) - deviance(fm1)
1076 : bates 879 newy <- simulate(fm2, n)
1077 : bates 888 mm1 <- model.matrix(fm1)
1078 :     mm2 <- model.matrix(fm2)
1079 : bates 879 ref <- numeric(n)
1080 : bates 888 mer1 <- as(fm1, "mer")
1081 :     mer2 <- as(fm2, "mer")
1082 : bates 879 for (j in 1:n) {
1083 : bates 888 .Call("lmer_update_y", mer2, newy[[j]], mm2, PACKAGE = "Matrix")
1084 :     LMEoptimize(mer2) <- cv
1085 :     .Call("lmer_update_y", mer1, newy[[j]], mm1, PACKAGE = "Matrix")
1086 :     LMEoptimize(mer1) <- cv
1087 :     ref[j] <- deviance(mer2) - deviance(mer1)
1088 : bates 879 }
1089 : bates 888 attr(ref, "observed") <- obs
1090 : bates 879 ref
1091 :     }
1092 : bates 940
1093 :     mer2 <-
1094 :     function(formula, data, family,
1095 :     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
1096 :     control = list(), start,
1097 :     subset, weights, na.action, offset,
1098 :     model = TRUE, x = FALSE, y = FALSE , ...)
1099 :     ## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
1100 :     {
1101 :     ## match and check parameters
1102 :     if (length(formula) < 3) stop("formula must be a two-sided formula")
1103 :     ## cv <- do.call("Matrix:::lmerControl", control)
1104 :     cv <- control
1105 : bates 978 cv$analyticGradient <- FALSE
1106 :     cv$msMaxIter <- as.integer(200)
1107 :     if (is.null(cv$msVerbose)) cv$msVerbose <- as.integer(1)
1108 : bates 981
1109 :     ## Must evaluate the model frame first and then fit the glm using
1110 :     ## that frame. Otherwise missing values in the grouping factors
1111 :     ## cause inconsistent numbers of observations.
1112 :    
1113 : bates 940 ## evaluate glm.fit, a generalized linear fit of fixed effects only
1114 :     mf <- match.call()
1115 :     m <- match(c("family", "data", "subset", "weights",
1116 :     "na.action", "offset"), names(mf), 0)
1117 : bates 981 mf <- fe <- mf[c(1, m)]
1118 : bates 940 frame.form <- subbars(formula) # substitute `+' for `|'
1119 :     fixed.form <- nobars(formula) # remove any terms with `|'
1120 :     if (inherits(fixed.form, "name")) # RHS is empty - use a constant
1121 :     fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
1122 :     environment(fixed.form) <- environment(frame.form) <- environment(formula)
1123 : bates 981
1124 :     ## evaluate a model frame for fixed and random effects
1125 :     mf$formula <- frame.form
1126 :     mf$family <- NULL
1127 :     mf$drop.unused.levels <- TRUE
1128 :     mf[[1]] <- as.name("model.frame")
1129 :     frm <- eval(mf, parent.frame())
1130 :    
1131 :     ## fit a glm model to the fixed formula
1132 :     fe$formula <- fixed.form
1133 :     fe$data <- frm
1134 :     fe$x <- fe$model <- fe$y <- TRUE
1135 :     fe[[1]] <- as.name("glm")
1136 :     glm.fit <- eval(fe, parent.frame())
1137 : bates 940 x <- glm.fit$x
1138 :     y <- as.double(glm.fit$y)
1139 :     family <- glm.fit$family
1140 : bates 981
1141 : bates 940 ## check for a linear mixed model
1142 :     lmm <- family$family == "gaussian" && family$link == "identity"
1143 :     if (lmm) { # linear mixed model
1144 :     method <- match.arg(method)
1145 :     if (method %in% c("PQL", "Laplace", "AGQ")) {
1146 :     warning(paste('Argument method = "', method,
1147 :     '" is not meaningful for a linear mixed model.\n',
1148 :     'Using method = "REML".\n', sep = ''))
1149 :     method <- "REML"
1150 :     }
1151 :     } else { # generalized linear mixed model
1152 :     if (missing(method)) method <- "PQL"
1153 :     else {
1154 :     method <- match.arg(method)
1155 :     if (method == "ML") method <- "PQL"
1156 :     if (method == "REML")
1157 :     warning('Argument method = "REML" is not meaningful ',
1158 :     'for a generalized linear mixed model.',
1159 :     '\nUsing method = "PQL".\n')
1160 :     }
1161 :     }
1162 : maechler 973
1163 : bates 940 ## grouping factors and model matrices for random effects
1164 :     bars <- findbars(formula[[3]])
1165 :     random <-
1166 :     lapply(bars, function(x)
1167 : bates 963 list(t(model.matrix(eval(substitute(~ expr,
1168 :     list(expr = x[[2]]))),
1169 :     frm)),
1170 : bates 940 eval(substitute(as.factor(fac)[,drop = TRUE],
1171 :     list(fac = x[[3]])), frm)))
1172 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
1173 : maechler 973
1174 : bates 940 ## order factor list by decreasing number of levels
1175 :     nlev <- sapply(random, function(x) length(levels(x[[2]])))
1176 :     if (any(diff(nlev) > 0)) {
1177 :     random <- random[rev(order(nlev))]
1178 :     }
1179 : maechler 973
1180 : bates 940 ## Create the model matrices and a mixed-effects representation (mer)
1181 : bates 963 mer <- .Call("mer2_create", random, x, y, method, PACKAGE = "Matrix")
1182 : bates 978 LMEoptimize(mer) <- cv
1183 : bates 940 mer
1184 :     }
1185 : bates 978
1186 :     ## Extract the permutation
1187 :     setAs("mer2", "pMatrix", function(from) .Call("mer2_pMatrix", from))
1188 :    
1189 :     ## Extract the L matrix
1190 :     setAs("mer2", "dtCMatrix", function(from) .Call("mer2_dtCMatrix", from))
1191 :    
1192 :     ## Optimization for mer2 objects
1193 :     setReplaceMethod("LMEoptimize", signature(x="mer2", value="list"),
1194 :     function(x, value)
1195 :     {
1196 :     if (value$msMaxIter < 1) return(x)
1197 :     nc <- x@nc
1198 : bates 984 constr <- unlist(lapply(nc[1:(length(nc) - 1)],
1199 : bates 978 function(k) 1:((k*(k+1))/2) <= k))
1200 :     fn <- function(pars)
1201 :     deviance(.Call("mer2_coefGets", x, pars, 2, PACKAGE = "Matrix"))
1202 :     gr <- NULL ## No gradient yet
1203 :     ## if (value$analyticGradient)
1204 :     ## function(pars) {
1205 :     ## if (!isTRUE(all.equal(pars,
1206 :     ## .Call("lmer_coef", x,
1207 :     ## 2, PACKAGE = "Matrix"))))
1208 :     ## .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
1209 :     ## .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
1210 :     ## }
1211 :     ## else NULL
1212 :     optimRes <- nlminb(.Call("mer2_coef", x, 2, PACKAGE = "Matrix"),
1213 :     fn, gr,
1214 :     lower = ifelse(constr, 5e-10, -Inf),
1215 :     control = list(iter.max = value$msMaxIter,
1216 :     trace = as.integer(value$msVerbose)))
1217 :     .Call("mer2_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
1218 :     if (optimRes$convergence != 0) {
1219 :     warning(paste("nlminb returned message",
1220 :     optimRes$message,"\n"))
1221 :     }
1222 :     return(x)
1223 :     })
1224 :    
1225 :     setMethod("deviance", "mer2",
1226 :     function(object, REML = NULL, ...) {
1227 :     .Call("mer2_factor", object, PACKAGE = "Matrix")
1228 :     if (is.null(REML))
1229 :     REML <- object@method == "REML"
1230 :     object@deviance[[ifelse(REML, "REML", "ML")]]
1231 :     })

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