Diff of /pkg/R/lmer.R

Parent Directory | Revision Log | Patch

	branches/trunk-lme4/R/lmer.R revision 689, Sat Apr 2 14:41:46 2005 UTC	pkg/R/lmer.R revision 1165, Sat Jan 14 23:41:55 2006 UTC
#	Line 1	Line 1
1	## Methods for lmer and for the objects that it produces	# Methods for lmer and for the objects that it produces
2
3	## Some utilities	## Some utilities
4
5	contr.SAS <- function(n, contrasts = TRUE)	## Return the pairs of expressions separated by vertical bars
6	## Eliminate this function after R-2.1.0 is released	findbars <- function(term)
7	contr.treatment(n,	{
8	if (is.numeric(n) && length(n) == 1) n else length(n),	if (is.name(term) || is.numeric(term)) return(NULL)
9	contrasts)	if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
10		if (!is.call(term)) stop("term must be of class call")
11	Lind <- function(i,j) {	if (term[[1]] == as.name('|')) return(term)
12	if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))	if (length(term) == 2) return(findbars(term[[2]]))
13	((i - 1) * i)/2 + j	c(findbars(term[[2]]), findbars(term[[3]]))
14	}	}
15
16	Dhalf <- function(from) {	## Return the formula omitting the pairs of expressions
17	D <- from@D	## that are separated by vertical bars
18	nf <- length(D)	nobars <- function(term)
19	Gp <- from@Gp	{
20	res <- array(0, rep(Gp[nf+1],2))	if (!('|' %in% all.names(term))) return(term)
21	for (i in 1:nf) {	if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
22	DD <- D[[i]]	if (length(term) == 2) {
23	dd <- dim(DD)	nb <- nobars(term[[2]])
24	for (k in 1:dd[3]) {	if (is.null(nb)) return(NULL)
25	mm <- array(DD[ , , k], dd[1:2])	term[[2]] <- nb
26	base <- Gp[i] + (k - 1)*dd[1]	return(term)
27	res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)	}
28	}	nb2 <- nobars(term[[2]])
29	}	nb3 <- nobars(term[[3]])
30	res	if (is.null(nb2)) return(nb3)
31	}	if (is.null(nb3)) return(nb2)
32		term[[2]] <- nb2
33	lmerControl <-                            # Control parameters for lmer	term[[3]] <- nb3
34	function(maxIter = 50,	term
35	msMaxIter = 50,	}
36	tolerance = sqrt((.Machine$double.eps)),
37	niterEM = 20,	## Substitute the '+' function for the '|' function
38	msTol = sqrt(.Machine$double.eps),	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		if (is.call(term) && term[[1]] == as.name('|'))
47		term[[1]] <- as.name('+')
48		term[[2]] <- subbars(term[[2]])
49		term[[3]] <- subbars(term[[3]])
50		term
51		}
52
53		## 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		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		## Control parameters for lmer
84		lmerControl <-
85		function(maxIter = 200, # used in ../src/lmer.c only
86		tolerance = sqrt(.Machine$double.eps),# ditto
87		msMaxIter = 200,
88		## msTol = sqrt(.Machine$double.eps),
89		## FIXME:  should be able to pass tolerances to nlminb()
90	msVerbose = getOption("verbose"),	msVerbose = getOption("verbose"),
91	PQLmaxIt = 20,	niterEM = 15,
92	EMverbose = getOption("verbose"),	EMverbose = getOption("verbose"),
93		PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead
94	analyticGradient = TRUE,	analyticGradient = TRUE,
95	analyticHessian=FALSE)	analyticHessian = FALSE # unused _FIXME_
96		)
97		{
98		list(maxIter = as.integer(maxIter),
99		tolerance = as.double(tolerance),
100		msMaxIter = as.integer(msMaxIter),
101		## msTol = as.double(msTol),
102		msVerbose = as.integer(msVerbose),# "integer" on purpose
103		niterEM = as.integer(niterEM),
104		EMverbose = as.logical(EMverbose),
105		PQLmaxIt = as.integer(PQLmaxIt),
106		analyticGradient = as.logical(analyticGradient),
107		analyticHessian = as.logical(analyticHessian))
108		}
109
110		setMethod("coef", signature(object = "lmer"),
111		function(object, ...)
112	{	{
113	list(maxIter = maxIter,	fef <- data.frame(rbind(object@fixed), check.names = FALSE)
114	msMaxIter = msMaxIter,	ref <- as(ranef(object), "list")
115	tolerance = tolerance,	names(ref) <- names(object@flist)
116	niterEM = niterEM,	val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
117	msTol = msTol,	for (i in seq(a = val)) {
118	msVerbose = msVerbose,	refi <- ref[[i]]
119	PQLmaxIt = PQLmaxIt,	row.names(val[[i]]) <- row.names(refi)
120	EMverbose=EMverbose,	if (!all(names(refi) %in% names(fef)))
121	analyticHessian=analyticHessian,	stop("unable to align random and fixed effects")
122	analyticGradient=analyticGradient)	val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
123	}	}
124		new("lmer.coef", val)
125		})
126
127		## setMethod("plot", signature(x = "lmer.coef"),
128		##           function(x, y, ...)
129		##       {
130		##           varying <- unique(do.call("c",
131		##                                     lapply(x, function(el)
132		##                                            names(el)[sapply(el,
133		##                                                             function(col)
134		##                                                             any(col != col[1]))])))
135		##           gf <- do.call("rbind", lapply(x, "[", j = varying))
136		##           gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
137		##           switch(min(length(varying), 3),
138		##                  qqmath(eval(substitute(~ x | .grp,
139		##                                         list(x = as.name(varying[1])))), gf, ...),
140		##                  xyplot(eval(substitute(y ~ x | .grp,
141		##                                         list(y = as.name(varying[1]),
142		##                                              x = as.name(varying[2])))), gf, ...),
143		##                  splom(~ gf | .grp, ...))
144		##       })
145
146		## setMethod("plot", signature(x = "lmer.ranef"),
147		##           function(x, y, ...)
148		##       {
149		##           lapply(x, function(x) {
150		##               cn <- lapply(colnames(x), as.name)
151		##               switch(min(ncol(x), 3),
152		##                      qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
153		##                      xyplot(eval(substitute(y ~ x,
154		##                                             list(y = cn[[1]],
155		##                                                  x = cn[[2]]))), x, ...),
156		##                      splom(~ x, ...))
157		##           })
158		##       })
159
160		setMethod("with", signature(data = "lmer"),
161		function(data, expr, ...) {
162		dat <- eval(data@call$data)
163		if (!is.null(na.act <- attr(data@frame, "na.action")))
164		dat <- dat[-na.act, ]
165		lst <- c(list(. = data), data@flist, data@frame, dat)
166		eval(substitute(expr), lst[unique(names(lst))])
167		})
168
169		setMethod("terms", signature(x = "lmer"),
170		function(x, ...) x@terms)
171
172		rWishart <- function(n, df, invScal)
173		.Call("Matrix_rWishart", n, df, invScal, PACKAGE = "Matrix")
174
175
176	setMethod("lmer", signature(formula = "formula", family = "missing"),	setMethod("lmer", signature(formula = "formula"),
177	function(formula, data, family,	function(formula, data, family,
178	method = c("REML", "ML", "PQL", "Laplace", "AGQ"),	method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
179	control = list(),	control = list(), start,
180	subset, weights, na.action, offset,	subset, weights, na.action, offset,
181	model = TRUE, x = FALSE, y = FALSE, ...)	model = TRUE, x = FALSE, y = FALSE, ...)
182	{	{
183	# match and check parameters	## match and check parameters
	REML <- match.arg(method) == "REML"
	controlvals <- do.call("lmerControl", control)
	controlvals$REML <- REML
184	if (length(formula) < 3) stop("formula must be a two-sided formula")	if (length(formula) < 3) stop("formula must be a two-sided formula")
185		cv <- do.call("lmerControl", control)
186
187	mf <- match.call()           # create the model frame as frm	## Must evaluate the model frame first and then fit the glm using
188	m <- match(c("data", "subset", "weights", "na.action", "offset"),	## that frame.  Otherwise missing values in the grouping factors
189	names(mf), 0)	## cause inconsistent numbers of observations.
190	mf <- mf[c(1, m)]	mf <- match.call()
191	mf[[1]] <- as.name("model.frame")	m <- match(c("family", "data", "subset", "weights",
192	frame.form <- subbars(formula)	"na.action", "offset"), names(mf), 0)
193	environment(frame.form) <- environment(formula)	mf <- fe <- mf[c(1, m)]
194		frame.form <- subbars(formula) # substitute `+' for `|'
195		fixed.form <- nobars(formula)  # remove any terms with `|'
196		if (inherits(fixed.form, "name")) # RHS is empty - use a constant
197		fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
198		environment(fixed.form) <- environment(frame.form) <- environment(formula)
199
200		## evaluate a model frame for fixed and random effects
201	mf$formula <- frame.form	mf$formula <- frame.form
202		mf$family <- NULL
203	mf$drop.unused.levels <- TRUE	mf$drop.unused.levels <- TRUE
204		mf[[1]] <- as.name("model.frame")
205	frm <- eval(mf, parent.frame())	frm <- eval(mf, parent.frame())
206
207	## grouping factors and model matrices for random effects	## fit a glm model to the fixed formula
208		fe$formula <- fixed.form
209		fe$subset <- NULL             # subset has already been created in call to data.frame
210		fe$data <- frm
211		fe$x <- fe$model <- fe$y <- TRUE
212		fe[[1]] <- as.name("glm")
213		glm.fit <- eval(fe, parent.frame())
214		x <- glm.fit$x
215		y <- as.double(glm.fit$y)
216		family <- glm.fit$family
217
218		## check for a linear mixed model
219		lmm <- family$family == "gaussian" && family$link == "identity"
220		if (lmm) { # linear mixed model
221		method <- match.arg(method)
222		if (method %in% c("PQL", "Laplace", "AGQ")) {
223		warning(paste('Argument method = "', method,
224		'" is not meaningful for a linear mixed model.\n',
225		'Using method = "REML".\n', sep = ''))
226		method <- "REML"
227		}
228		} else { # generalized linear mixed model
229		if (missing(method)) method <- "PQL"
230		else {
231		method <- match.arg(method)
232		if (method == "ML") method <- "PQL"
233		if (method == "REML")
234		warning('Argument method = "REML" is not meaningful ',
235		'for a generalized linear mixed model.',
236		'\nUsing method = "PQL".\n')
237		}
238		}
239		## create factor list for the random effects
240	bars <- findbars(formula[[3]])	bars <- findbars(formula[[3]])
241	random <-	names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
242	lapply(bars,	fl <- lapply(bars,
243	function(x) list(model.matrix(eval(substitute(~term,	function(x)
	list(term=x[[2]]))),
	frm),
244	eval(substitute(as.factor(fac)[,drop = TRUE],	eval(substitute(as.factor(fac)[,drop = TRUE],
245	list(fac = x[[3]])), frm)))	list(fac = x[[3]])), frm))
	names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
246	## order factor list by decreasing number of levels	## order factor list by decreasing number of levels
247	nlev <- sapply(random, function(x) length(levels(x[[2]])))	nlev <- sapply(fl, function(x) length(levels(x)))
248	if (any(diff(nlev) > 0)) {	if (any(diff(nlev) > 0)) {
249	random <- random[rev(order(nlev))]	ord <- rev(order(nlev))
250	}	bars <- bars[ord]
251	fixed.form <- nobars(formula)	fl <- fl[ord]
252	if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula	}
253	Xmat <- model.matrix(fixed.form, frm)	## create list of transposed model matrices for random effects
254	mmats <- c(lapply(random, "[[", 1),	Ztl <- lapply(bars, function(x)
255	.fixed = list(cbind(Xmat, .response = model.response(frm))))	t(model.matrix(eval(substitute(~ expr,
256	## FIXME: Use Xfrm and Xmat to get the terms and assign	list(expr = x[[2]]))),
257	## slots, pass them to lmer_create, then destroy them	frm)))
258	obj <- .Call("lmer_create", lapply(random, "[[", 2),	## Create the mixed-effects representation (mer) object
259	mmats, PACKAGE = "Matrix")	mer <- .Call("mer_create", fl,
260	slot(obj, "terms") <- attr(model.frame(fixed.form, data), "terms")	.Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"),
261	slot(obj, "assign") <- attr(Xmat, "assign")	x, y, method, sapply(Ztl, nrow),
262	slot(obj, "call") <- match.call()	c(lapply(Ztl, rownames), list(.fixed = colnames(x))),
263	slot(obj, "REML") <- REML	!(family$family %in% c("binomial", "poisson")),
264	rm(Xmat)	match.call(), family,
	.Call("lmer_initial", obj, PACKAGE="Matrix")
	.Call("lmer_ECMEsteps", obj,
	controlvals$niterEM,
	controlvals$REML,
	controlvals$EMverbose,
265	PACKAGE = "Matrix")	PACKAGE = "Matrix")
266	LMEoptimize(obj) <- controlvals	if (lmm) {
267	slot(obj, "residuals") <-	.Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
268	unname(model.response(frm) -	LMEoptimize(mer) <- cv
269	(slot(obj, "fitted") <-	return(mer)
270	.Call("lmer_fitted", obj, mmats, TRUE, PACKAGE = "Matrix")))	}
271	obj
272	})	## The rest of the function applies to generalized linear mixed models
273		gVerb <- getOption("verbose")
274		eta <- glm.fit$linear.predictors
275		wts <- glm.fit$prior.weights
276		wtssqr <- wts * wts
277		offset <- glm.fit$offset
278		if (is.null(offset)) offset <- numeric(length(eta))
279		linkinv <- quote(family$linkinv(eta))
280		mu.eta <- quote(family$mu.eta(eta))
281		mu <- family$linkinv(eta)
282		variance <- quote(family$variance(mu))
283		dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
284		LMEopt <- get("LMEoptimize<-")
285		doLMEopt <- quote(LMEopt(x = mer, value = cv))
286
287		GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
288		.Call("glmer_PQL", GSpt, PACKAGE = "Matrix")  # obtain PQL estimates
289		fixInd <- seq(ncol(x))
290		## pars[fixInd] == beta, pars[-fixInd] == theta
291		PQLpars <- c(fixef(mer),
292		.Call("mer_coef", mer, 2, PACKAGE = "Matrix"))
293		.Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")
294		## indicator of constrained parameters
295		const <- c(rep(FALSE, length(fixInd)),
296		unlist(lapply(mer@nc[seq(along = fl)],
297		function(k) 1:((k*(k+1))/2) <= k)
298		))
299		devLaplace <- function(pars)
300		.Call("glmer_devLaplace", pars, GSpt, PACKAGE = "Matrix")
301
302		optimRes <-
303		nlminb(PQLpars, devLaplace,
304		lower = ifelse(const, 5e-10, -Inf),
305		control = list(trace = getOption("verbose"),
306		iter.max = cv$msMaxIter))
307		.Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
308		return(mer)
309		deviance <- devAGQ(PQLpars, 1)
310
311		### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs
312		### AGQ for nc > 1 first.
313		fxd <- PQLpars[fixInd]
314		loglik <- logLik(mer)
315
316		if (method %in% c("Laplace", "AGQ")) {
317		nAGQ <- 1
318		if (method == "AGQ") {    # determine nAGQ at PQL estimates
319		dev11 <- devAGQ(PQLpars, 11)
320		## FIXME: Should this be an absolute or a relative tolerance?
321		devTol <- sqrt(.Machine$double.eps) * abs(dev11)
322		for (nAGQ in c(9, 7, 5, 3, 1))
323		if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break
324		nAGQ <- nAGQ + 2
325		if (gVerb)
326		cat(paste("Using", nAGQ, "quadrature points per column\n"))
327		}
328		obj <- function(pars)
329		.Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
330		optimRes <-
331		nlminb(PQLpars, obj,
332		lower = ifelse(const, 5e-10, -Inf),
333		control = list(trace = getOption("verbose"),
334		iter.max = cv$msMaxIter))
335		optpars <- optimRes$par
336		if (optimRes$convergence != 0)
337		warning("nlminb failed to converge")
338		deviance <- optimRes$objective
339		if (gVerb)
340		cat(paste("convergence message", optimRes$message, "\n"))
341		fxd[] <- optpars[fixInd]  ## preserve the names
342		.Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
343		}
344
345		.Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
346		loglik[] <- -deviance/2
347		new("lmer", mer,
348		frame = if (model) frm else data.frame(),
349		terms = glm.fit$terms,
350		assign = attr(glm.fit$x, "assign"),
351		call = match.call(), family = family,
352		logLik = loglik, fixed = fxd)
353		})
354
355
356		## Extract the permutation
357		setAs("mer", "pMatrix", function(from)
358		.Call("mer_pMatrix", from, PACKAGE = "Matrix"))
359
360		## Extract the L matrix
361		setAs("mer", "dtCMatrix", function(from)
362		.Call("mer_dtCMatrix", from, PACKAGE = "Matrix"))
363
364		## Extract the fixed effects
365		setMethod("fixef", signature(object = "mer"),
366		function(object, ...)
367		.Call("mer_fixef", object, PACKAGE = "Matrix"))
368
369		## Extract the random effects
370		setMethod("ranef", signature(object = "mer"),
371		function(object, ...)
372		.Call("mer_ranef", object, PACKAGE = "Matrix")
373		)
374
375	setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),	## Optimization for mer objects
376		setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
377	function(x, value)	function(x, value)
378	{	{
379	if (value$msMaxIter < 1) return(x)	if (value$msMaxIter < 1) return(x)
	st <- ccoef(x)         # starting values
380	nc <- x@nc	nc <- x@nc
	nc <- nc[1:(length(nc) - 2)]
381	constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))	constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
382	fn <- function(pars) {	fn <- function(pars)
383	ccoef(x) <- pars	deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
	deviance(x, REML = value$REML)
	}
384	gr <- if (value$analyticGradient)	gr <- if (value$analyticGradient)
385	function(pars) {	function(pars) {
386	if (!identical(TRUE,all.equal(pars, ccoef(x)))) ccoef(x) <- pars	if (!isTRUE(all.equal(pars,
387	grad <- gradient(x, REML = value$REML, unconst = TRUE)	.Call("mer_coef", x,
388	grad[constr] <- -grad[constr]/pars[constr]	2, PACKAGE = "Matrix"))))
389	grad	.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")
390	} else NULL	.Call("mer_gradient", x, 2, PACKAGE = "Matrix")
391	optimRes <- optim(st, fn, gr,	}
392	method = "L-BFGS-B",	else NULL
393	lower = ifelse(constr, 1e-10, -Inf),	optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"),
394	control = list(maxit = value$msMaxIter,	fn, gr,
395		lower = ifelse(constr, 5e-10, -Inf),
396		control = list(iter.max = value$msMaxIter,
397	trace = as.integer(value$msVerbose)))	trace = as.integer(value$msVerbose)))
398		.Call("mer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
399	if (optimRes$convergence != 0) {	if (optimRes$convergence != 0) {
400	warning(paste("optim returned message",optimRes$message,"\n"))	warning(paste("nlminb returned message",
401		optimRes$message,"\n"))
402	}	}
	ccoef(x) <- optimRes$par
403	return(x)	return(x)
404	})	})
405
406	setMethod("ranef", signature(object = "lmer"),	setMethod("deviance", signature(object = "mer"),
	function(object, accumulate = FALSE, ...) {
	val <- new("lmer.ranef",
	lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),
	data.frame, check.names = FALSE),
	varFac = object@bVar,
	stdErr = .Call("lmer_sigma", object,
	object@REML, PACKAGE = "Matrix"))
	if (!accumulate || length(val@varFac) == 1) return(val)
	## check for nested factors
	L <- object@L
	if (any(sapply(seq(a = val), function(i) length(L[[Lind(i,i)]]@i))))
	error("Require nested grouping factors to accumulate random effects")
	val
	})
	setMethod("fixef", signature(object = "lmer"),
407	function(object, ...) {	function(object, ...) {
408	val <- .Call("lmer_fixef", object, PACKAGE = "Matrix")	.Call("mer_factor", object, PACKAGE = "Matrix")
409	val[-length(val)]	object@deviance[[ifelse(object@method == "REML", "REML", "ML")]]
410	})	})
411
412	setMethod("VarCorr", signature(x = "lmer"),	setMethod("mcmcsamp", signature(object = "mer"),
413	function(x, REML = TRUE, useScale = TRUE, ...) {	function(object, n = 1, verbose = FALSE, saveb = FALSE,
414	val <- .Call("lmer_variances", x, PACKAGE = "Matrix")	trans = TRUE, ...)
415	for (i in seq(along = val)) {	{
416	dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])	ans <- t(.Call("mer_MCMCsamp", object, saveb, n,
417	val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")	trans, PACKAGE = "Matrix"))
418	}	attr(ans, "mcpar") <- as.integer(c(1, n, 1))
419	new("VarCorr",	class(ans) <- "mcmc"
420	scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),	glmer <- FALSE
421	reSumry = val,	gnms <- names(object@flist)
422	useScale = useScale)	cnms <- object@cnames
423		ff <- fixef(object)
424		colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
425		unlist(lapply(seq(along = gnms),
426		function(i)
427		abbrvNms(gnms[i],cnms[[i]]))))
428		if (trans) {
429		## parameter type: 0 => fixed effect, 1 => variance,
430		##                 2 => covariance
431		ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1,
432		unlist(lapply(seq(along = gnms),
433		function(i)
434		{
435		k <- length(cnms[[i]])
436		rep(1:2, c(k, (k*(k-1))/2))
437		})))
438		colnms[ptyp == 1] <-
439		paste("log(", colnms[ptyp == 1], ")", sep = "")
440		colnms[ptyp == 2] <-
441		paste("atanh(", colnms[ptyp == 2], ")", sep = "")
442		}
443		colnames(ans) <- colnms
444		ans
445		})
446
447		setMethod("simulate", signature(object = "mer"),
448		function(object, nsim = 1, seed = NULL, ...)
449		{
450		if(!exists(".Random.seed", envir = .GlobalEnv))
451		runif(1)               # initialize the RNG if necessary
452		if(is.null(seed))
453		RNGstate <- .Random.seed
454		else {
455		R.seed <- .Random.seed
456		set.seed(seed)
457		RNGstate <- structure(seed, kind = as.list(RNGkind()))
458		on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
459		}
460
461		family <- object@family
462		if (family$family != "gaussian" ||
463		family$link != "identity")
464		stop("simulation of generalized linear mixed models not yet implemented")
465		## similate the linear predictors
466		lpred <- .Call("mer_simulate", object, nsim, PACKAGE = "Matrix")
467		sc <- 1
468		if (object@useScale)
469		sc <- .Call("mer_sigma", object, object@method == "REML",
470		PACKAGE = "Matrix")
471		## add fixed-effects contribution and per-observation noise term
472		lpred <- as.data.frame(lpred + drop(object@X %*% fixef(object)) +
473		rnorm(prod(dim(lpred)), sd = sc))
474		## save the seed
475		attr(lpred, "seed") <- RNGstate
476		lpred
477	})	})
478
	setMethod("gradient", signature(x = "lmer"),
	function(x, REML, unconst, ...)
	.Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix"))
	setMethod("summary", signature(object = "lmer"),
	function(object, ...)
	new("summary.lmer", object, useScale = TRUE, showCorrelation = TRUE))
	setMethod("show", signature(object = "lmer"),
	function(object)
	show(new("summary.lmer", object, useScale = TRUE,
	showCorrelation = FALSE))
	)
479
480	setMethod("show", "summary.lmer",	setMethod("show", "mer",
481	function(object) {	function(object) {
482	fcoef <- fixef(object)	vcShow <- function(varc, useScale)
483		{
484		digits <- max(3, getOption("digits") - 2)
485		sc <- attr(varc, "sc")
486		recorr <- lapply(varc, function(el) el@factors$correlation)
487		reStdDev <- c(lapply(recorr, slot, "sd"), list(Residual = sc))
488		reLens <- unlist(c(lapply(reStdDev, length)))
489		reMat <- array('', c(sum(reLens), 4),
490		list(rep('', sum(reLens)),
491		c("Groups", "Name", "Variance", "Std.Dev.")))
492		reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
493		reMat[,2] <- c(unlist(lapply(reStdDev, names)), "")
494		reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
495		reMat[,4] <- format(unlist(reStdDev), digits = digits)
496		if (any(reLens > 1)) {
497		maxlen <- max(reLens)
498		corr <-
499		do.call("rbind",
500		lapply(recorr,
501		function(x, maxlen) {
502		x <- as(x, "matrix")
503		cc <- format(round(x, 3), nsmall = 3)
504		cc[!lower.tri(cc)] <- ""
505		nr <- dim(cc)[1]
506		if (nr >= maxlen) return(cc)
507		cbind(cc, matrix("", nr, maxlen-nr))
508		}, maxlen))
509		colnames(corr) <- c("Corr", rep("", maxlen - 1))
510		reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr))))
511		}
512		if (!useScale) reMat <- reMat[-nrow(reMat),]
513		print(reMat, quote = FALSE)
514		}
515
516		fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")
517	useScale <- object@useScale	useScale <- object@useScale
518	corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"),	corF <- vcov(object)@factors$correlation
	"corrmatrix")
519	DF <- getFixDF(object)	DF <- getFixDF(object)
520	coefs <- cbind(fcoef, corF@stdDev, DF)	coefs <- cbind(fcoef, corF@sd, DF)
	nc <- object@nc
521	dimnames(coefs) <-	dimnames(coefs) <-
522	list(names(fcoef), c("Estimate", "Std. Error", "DF"))	list(names(fcoef), c("Estimate", "Std. Error", "DF"))
523	digits <- max(3, getOption("digits") - 2)	digits <- max(3, getOption("digits") - 2)
524	REML <- length(object@REML) > 0 && object@REML[1]	REML <- object@method == "REML"
525	llik <- logLik(object)	llik <- logLik(object, REML)
526	dev <- object@deviance	dev <- object@deviance
527		devc <- object@devComp
528
529	rdig <- 5	rdig <- 5
530		if (glz <- !(object@method %in% c("REML", "ML"))) {
531		cat(paste("Generalized linear mixed model fit using",
532		object@method, "\n"))
533		} else {
534	cat("Linear mixed-effects model fit by ")	cat("Linear mixed-effects model fit by ")
535	cat(ifelse(object@REML, "REML\n", "maximum likelihood\n") )	cat(if(REML) "REML\n" else "maximum likelihood\n")
536		}
537	if (!is.null(object@call$formula)) {	if (!is.null(object@call$formula)) {
538	cat("Formula:", deparse(object@call$formula),"\n")	cat("Formula:", deparse(object@call$formula),"\n")
539	}	}
#	Line 231	Line 544
544	cat(" Subset:",	cat(" Subset:",
545	deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")	deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
546	}	}
547		if (glz) {
548		cat(" Family: ", object@family$family, "(",
549		object@family$link, " link)\n", sep = "")
550		print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
551		logLik = c(llik),
552		deviance = -2*llik,
553		row.names = ""))
554		} else {
555	print(data.frame(AIC = AIC(llik), BIC = BIC(llik),	print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
556	logLik = c(llik),	logLik = c(llik),
557	MLdeviance = dev["ML"],	MLdeviance = dev["ML"],
558	REMLdeviance = dev["REML"],	REMLdeviance = dev["REML"],
559	row.names = ""))	row.names = ""))
560		}
561	cat("Random effects:\n")	cat("Random effects:\n")
562	show(VarCorr(object))	vcShow(VarCorr(object), useScale)
563	ngrps <- lapply(object@flist, function(x) length(levels(x)))	ngrps <- lapply(object@flist, function(x) length(levels(x)))
564	cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))	cat(sprintf("# of obs: %d, groups: ", devc[1]))
565	cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))	cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
566	cat("\n")	cat("\n")
567	if (!useScale)	if (!useScale)
568	cat("\nEstimated scale (compare to 1) ",	cat("\nEstimated scale (compare to 1) ",
569	.Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"),	.Call("mer_sigma", object, FALSE, PACKAGE = "Matrix"),
570	"\n")	"\n")
571	if (nrow(coefs) > 0) {	if (nrow(coefs) > 0) {
572	if (useScale) {	if (useScale) {
#	Line 263	Line 585
585	}	}
586	cat("\nFixed effects:\n")	cat("\nFixed effects:\n")
587	printCoefmat(coefs, tst.ind = 4, zap.ind = 3)	printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
	if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
588	rn <- rownames(coefs)	rn <- rownames(coefs)
	dimnames(corF) <- list(
	abbreviate(rn, minlen=11),
	abbreviate(rn, minlen=6))
589	if (!is.null(corF)) {	if (!is.null(corF)) {
590	p <- NCOL(corF)	p <- ncol(corF)
591	if (p > 1) {	if (p > 1) {
592	cat("\nCorrelation of Fixed Effects:\n")	cat("\nCorrelation of Fixed Effects:\n")
593	corF <- format(round(corF, 3), nsmall = 3)	corF <- matrix(format(round(corF@x, 3), nsmall = 3),
594		nc = p)
595		dimnames(corF) <- list(
596		abbreviate(rn, minlen=11),
597		abbreviate(rn, minlen=6))
598	corF[!lower.tri(corF)] <- ""	corF[!lower.tri(corF)] <- ""
599	print(corF[-1, -p, drop=FALSE], quote = FALSE)	print(corF[-1, -p, drop=FALSE], quote = FALSE)
600	}	}
601	}	}
602	}	}
	}
603	invisible(object)	invisible(object)
604	})	})
605
606	setMethod("lmer", signature(formula = "formula"),	setMethod("vcov", signature(object = "mer"),
607	function(formula, family, data,	function(object, REML = object@method == "REML",
608	method = c("REML", "ML", "PQL", "Laplace", "AGQ"),	useScale = object@useScale,...) {
609	control = list(),	sc <- if (object@useScale) {
610	subset, weights, na.action, offset,	.Call("mer_sigma", object, REML, PACKAGE = "Matrix")
611	model = TRUE, x = FALSE, y = FALSE, ...)	} else { 1 }
612	{	rr <- as(sc^2 * tcrossprod(solve(object@RXX)), "dpoMatrix")
613	gVerb <- getOption("verbose")	rr@factors$correlation <- as(rr, "correlation")
614	# match and check parameters	rr
	controlvals <- do.call("lmerControl", control)
	controlvals$REML <- FALSE
	if (length(formula) < 3) stop("formula must be a two-sided formula")
	## initial glm fit
	mf <- match.call()
	m <- match(c("family", "data", "subset", "weights",
	"na.action", "offset"),
	names(mf), 0)
	mf <- mf[c(1, m)]
	mf[[1]] <- as.name("glm")
	fixed.form <- nobars(formula)
	if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula
	environment(fixed.form) <- environment(formula)
	mf$formula <- fixed.form
	mf$x <- mf$model <- mf$y <- TRUE
	glm.fit <- eval(mf, parent.frame())
	family <- glm.fit$family
	## Note: offset is on the linear scale
	offset <- glm.fit$offset
	if (is.null(offset)) offset <- 0
	weights <- sqrt(abs(glm.fit$prior.weights))
	## initial 'fitted' values on linear scale
	etaold <- eta <- glm.fit$linear.predictors
	## evaluation of model frame
	mf$x <- mf$model <- mf$y <- mf$family <- NULL
	mf$drop.unused.levels <- TRUE
	this.form <- subbars(formula)
	environment(this.form) <- environment(formula)
	mf$formula <- this.form
	mf[[1]] <- as.name("model.frame")
	frm <- eval(mf, parent.frame())
	## grouping factors and model matrices for random effects
	bars <- findbars(formula[[3]])
	random <-
	lapply(bars,
	function(x) list(model.matrix(eval(substitute(~term,
	list(term=x[[2]]))),
	frm),
	eval(substitute(as.factor(fac)[,drop = TRUE],
	list(fac = x[[3]])), frm)))
	names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
	## order factor list by decreasing number of levels
	nlev <- sapply(random, function(x) length(levels(x[[2]])))
	if (any(diff(nlev) > 0)) {
	random <- random[rev(order(nlev))]
	}
	mmats <- c(lapply(random, "[[", 1),
	.fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
	## FIXME: Use Xfrm and Xmat to get the terms and assign
	## slots, pass these to lmer_create, then destroy Xfrm, Xmat, etc.
	obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")
	obj@terms <- attr(glm.fit$model, "terms")
	obj@assign <- attr(glm.fit$x, "assign")
	obj@call <- match.call()
	obj@REML <- FALSE
	rm(glm.fit)
	.Call("lmer_initial", obj, PACKAGE="Matrix")
	mmats.unadjusted <- mmats
	mmats[[1]][1,1] <- mmats[[1]][1,1]
	conv <- FALSE
	firstIter <- TRUE
	msMaxIter.orig <- controlvals$msMaxIter
	responseIndex <- ncol(mmats$.fixed)
	for (iter in seq(length = controlvals$PQLmaxIt))
	{
	mu <- family$linkinv(eta)
	dmu.deta <- family$mu.eta(eta)
	## weights (note: weights is already square-rooted)
	w <- weights * dmu.deta / sqrt(family$variance(mu))
	## adjusted response (should be comparable to X \beta, not including offset
	z <- eta - offset + (mmats.unadjusted$.fixed[, responseIndex] - mu) / dmu.deta
	.Call("nlme_weight_matrix_list",
	mmats.unadjusted, w, z, mmats, PACKAGE="Matrix")
	.Call("lmer_update_mm", obj, mmats, PACKAGE="Matrix")
	if (firstIter) {
	.Call("lmer_initial", obj, PACKAGE="Matrix")
	if (gVerb) cat(" PQL iterations convergence criterion\n")
	}
	.Call("lmer_ECMEsteps", obj,
	controlvals$niterEM,
	FALSE,
	controlvals$EMverbose,
	PACKAGE = "Matrix")
	LMEoptimize(obj) <- controlvals
	eta[] <- offset + ## FIXME: should the offset be here ?
	.Call("lmer_fitted", obj,
	mmats.unadjusted, TRUE, PACKAGE = "Matrix")
	crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))
	if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))
	## use this to determine convergence
	if (crit < controlvals$tolerance) {
	conv <- TRUE
	break
	}
	etaold[] <- eta
	## Changing number of iterations on second and
	## subsequent iterations.
	if (firstIter)
	{
	controlvals$niterEM <- 2
	controlvals$msMaxIter <- 10
	firstIter <- FALSE
	}
	}
	if (!conv) warning("IRLS iterations for glmm did not converge")
	obj
615	})	})
616
617	## calculates degrees of freedom for fixed effects Wald tests	## calculates degrees of freedom for fixed effects Wald tests
#	Line 410	Line 619
619	## be very tricky to decide what a 'right' answer should be with	## be very tricky to decide what a 'right' answer should be with
620	## crossed random effects.	## crossed random effects.
621
622	setMethod("getFixDF", signature(object="lmer"),	setMethod("getFixDF", signature(object="mer"),
623	function(object, ...)	function(object, ...) {
624	{	devc <- object@devComp
625	nc <- object@nc[-seq(along = object@Omega)]	rep(as.integer(devc[1]- devc[2]), devc[2])
	p <- nc[1] - 1
	n <- nc[2]
	rep(n-p, p)
626	})	})
627
628	setMethod("logLik", signature(object="lmer"),	setMethod("logLik", signature(object="mer"),
629	function(object, REML = object@REML, ...) {	function(object, REML = object@method == "REML", ...) {
630	val <- -deviance(object, REML = REML)/2	val <- -deviance(object, REML = REML)/2
631	nc <- object@nc[-seq(a = object@Omega)]	devc <- as.integer(object@devComp[1:2])
632	attr(val, "nall") <- attr(val, "nobs") <- nc[2]	attr(val, "nall") <- attr(val, "nobs") <- devc[1]
633	attr(val, "df") <- nc[1] + length(ccoef(object))	attr(val, "df") <- abs(devc[2]) +
634		length(.Call("mer_coef", object, 0, PACKAGE = "Matrix"))
635	attr(val, "REML") <- REML	attr(val, "REML") <- REML
636	class(val) <- "logLik"	class(val) <- "logLik"
637	val	val
638	})	})
639
640	setMethod("anova", signature(object = "lmer"),	setMethod("VarCorr", signature(x = "mer"),
641		function(x, REML = x@method == "REML", useScale = x@useScale, ...)
642		{
643		sc <- 1
644		if (useScale)
645		sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix")
646		sc2 <- sc * sc
647		ans <- lapply(x@Omega, function(el) {
648		el <- as(sc2 * solve(el), "dpoMatrix")
649		el@factors$correlation <- as(el, "correlation")
650		el
651		})
652		attr(ans, "sc") <- sc
653		ans
654		})
655
656		setMethod("anova", signature(object = "mer"),
657	function(object, ...)	function(object, ...)
658	{	{
659	mCall <- match.call(expand.dots = TRUE)	mCall <- match.call(expand.dots = TRUE)
660	dots <- list(...)	dots <- list(...)
661	modp <- logical(0)	modp <- logical(0)
662	if (length(dots))	if (length(dots))
663	modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm")	modp <- sapply(dots, inherits, "mer") | sapply(dots, inherits, "lm")
664	if (any(modp)) {              # multiple models - form table	if (any(modp)) {              # multiple models - form table
665	opts <- dots[!modp]	opts <- dots[!modp]
666	mods <- c(list(object), dots[modp])	mods <- c(list(object), dots[modp])
#	Line 467	Line 690
690	check.names = FALSE)	check.names = FALSE)
691	class(val) <- c("anova", class(val))	class(val) <- c("anova", class(val))
692	attr(val, "heading") <-	attr(val, "heading") <-
693	c(header, "", "Models:",	c(header, "Models:",
694	paste(names(mods),	paste(names(mods),
695	unlist(lapply(lapply(calls, "[[", "formula"), deparse)),	unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
696	sep = ": "),"")	sep = ": "))
697	return(val)	return(val)
698	} else {	} else {
699	foo <- object	foo <- object
700	foo@status["factored"] <- FALSE	foo@status["factored"] <- FALSE
701	.Call("lmer_factor", foo, PACKAGE="Matrix")	.Call("mer_factor", foo, PACKAGE="Matrix")
702	dfr <- getFixDF(foo)	dfr <- getFixDF(foo)
703	rcol <- ncol(foo@RXX)	ss <- foo@rXy^2
704	ss <- foo@RXX[ , rcol]^2	ssr <- exp(foo@devComp["logryy2"])
	ssr <- ss[[rcol]]
705	ss <- ss[seq(along = dfr)]	ss <- ss[seq(along = dfr)]
706	names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]	names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
707	asgn <- foo@assign	asgn <- foo@assign
#	Line 489	Line 711
711	nmeffects <- c("(Intercept)", nmeffects)	nmeffects <- c("(Intercept)", nmeffects)
712	ss <- unlist(lapply(split(ss, asgn), sum))	ss <- unlist(lapply(split(ss, asgn), sum))
713	df <- unlist(lapply(split(asgn,  asgn), length))	df <- unlist(lapply(split(asgn,  asgn), length))
714	dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))	#dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1]))
715	ms <- ss/df	ms <- ss/df
716	f <- ms/(ssr/dfr)	#f <- ms/(ssr/dfr)
717	P <- pf(f, df, dfr, lower.tail = FALSE)	#P <- pf(f, df, dfr, lower.tail = FALSE)
718	table <- data.frame(df, ss, ms, dfr, f, P)	#table <- data.frame(df, ss, ms, dfr, f, P)
719		table <- data.frame(df, ss, ms)
720	dimnames(table) <-	dimnames(table) <-
721	list(nmeffects,	list(nmeffects,
722	c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))	#                       c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
723		c("Df", "Sum Sq", "Mean Sq"))
724	if ("(Intercept)" %in% nmeffects) table <- table[-1,]	if ("(Intercept)" %in% nmeffects) table <- table[-1,]
725	attr(table, "heading") <- "Analysis of Variance Table"	attr(table, "heading") <- "Analysis of Variance Table"
726	class(table) <- c("anova", "data.frame")	class(table) <- c("anova", "data.frame")
#	Line 504	Line 728
728	}	}
729	})	})
730
731	setMethod("update", signature(object = "lmer"),	setMethod("confint", signature(object = "mer"),
	function(object, formula., ..., evaluate = TRUE)
	{
	call <- object@call
	if (is.null(call))
	stop("need an object with call component")
	extras <- match.call(expand.dots = FALSE)$...
	if (!missing(formula.))
	call$formula <- update.formula(formula(object), formula.)
	if (length(extras) > 0) {
	existing <- !is.na(match(names(extras), names(call)))
	for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
	if (any(!existing)) {
	call <- c(as.list(call), extras[!existing])
	call <- as.call(call)
	}
	}
	if (evaluate)
	eval(call, parent.frame())
	else call
	})
	setMethod("confint", signature(object = "lmer"),
732	function (object, parm, level = 0.95, ...)	function (object, parm, level = 0.95, ...)
733	{	.NotYetImplemented()
	cf <- fixef(object)
	pnames <- names(cf)
	if (missing(parm))
	parm <- seq(along = pnames)
	else if (is.character(parm))
	parm <- match(parm, pnames, nomatch = 0)
	a <- (1 - level)/2
	a <- c(a, 1 - a)
	pct <- paste(round(100 * a, 1), "%")
	ci <- array(NA, dim = c(length(parm), 2),
	dimnames = list(pnames[parm], pct))
	ses <- sqrt(diag(vcov(object)))[parm]
	ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
	ci
	})
	setMethod("param", signature(object = "lmer"),
	function(object, unconst = FALSE, ...) {
	.Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
	})
	setMethod("deviance", "lmer",
	function(object, REML = NULL, ...) {
	.Call("lmer_factor", object, PACKAGE = "Matrix")
	if (is.null(REML))
	REML <- if (length(oR <- object@REML)) oR else FALSE
	object@deviance[[ifelse(REML, "REML", "ML")]]
	})
	setMethod("chol", signature(x = "lmer"),
	function(x, pivot = FALSE, LINPACK = pivot) {
	x@status["factored"] <- FALSE # force a decomposition
	.Call("lmer_factor", x, PACKAGE = "Matrix")
	})
	setMethod("solve", signature(a = "lmer", b = "missing"),
	function(a, b, ...)
	.Call("lmer_invert", a, PACKAGE = "Matrix")
734	)	)
735
736	setMethod("formula", "lmer", function(x, ...) x@call$formula)	setMethod("fitted", signature(object = "mer"),
737		function(object, ...)
738	setMethod("vcov", signature(object = "lmer"),	.Call("mer_fitted", object, PACKAGE = "Matrix")
739	function(object, REML = object@REML, useScale = TRUE,...) {	)
	sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
	rr <- object@RXX
	nms <- object@cnames[[".fixed"]]
	dimnames(rr) <- list(nms, nms)
	nr <- nrow(rr)
	rr <- rr[-nr, -nr, drop = FALSE]
	rr <- rr %*% t(rr)
	if (useScale) {
	rr = sc^2 * rr
	}
	rr
	})
740
741	## Extract the L matrix	setMethod("formula", signature(x = "mer"),
742	setAs("lmer", "dtTMatrix",	function(x, ...)
743	function(from)	x@call$formula
744	{	)
	## force a refactorization if the factors have been inverted
	if (from@status["inverted"]) from@status["factored"] <- FALSE
	.Call("lmer_factor", from, PACKAGE = "Matrix")
	L <- lapply(from@L, as, "dgTMatrix")
	nf <- length(from@D)
	Gp <- from@Gp
	nL <- Gp[nf + 1]
	Li <- integer(0)
	Lj <- integer(0)
	Lx <- double(0)
	for (i in 1:nf) {
	for (j in 1:i) {
	Lij <- L[[Lind(i, j)]]
	Li <- c(Li, Lij@i + Gp[i])
	Lj <- c(Lj, Lij@j + Gp[j])
	Lx <- c(Lx, Lij@x)
	}
	}
	new("dtTMatrix", Dim = as.integer(c(nL, nL)), i = Li, j = Lj, x = Lx,
	uplo = "L", diag = "U")
	})
	## Extract the ZZX matrix
	setAs("lmer", "dsTMatrix",
	function(from)
	{
	.Call("lmer_inflate", from, PACKAGE = "Matrix")
	ZZpO <- lapply(from@ZZpO, as, "dgTMatrix")
	ZZ <- lapply(from@ZtZ, as, "dgTMatrix")
	nf <- length(ZZpO)
	Gp <- from@Gp
	nZ <- Gp[nf + 1]
	Zi <- integer(0)
	Zj <- integer(0)
	Zx <- double(0)
	for (i in 1:nf) {
	ZZpOi <- ZZpO[[i]]
	Zi <- c(Zi, ZZpOi@i + Gp[i])
	Zj <- c(Zj, ZZpOi@j + Gp[i])
	Zx <- c(Zx, ZZpOi@x)
	if (i > 1) {
	for (j in 1:(i-1)) {
	ZZij <- ZZ[[Lind(i, j)]]
	## off-diagonal blocks are transposed
	Zi <- c(Zi, ZZij@j + Gp[j])
	Zj <- c(Zj, ZZij@i + Gp[i])
	Zx <- c(Zx, ZZij@x)
	}
	}
	}
	new("dsTMatrix", Dim = as.integer(c(nZ, nZ)), i = Zi, j = Zj, x = Zx,
	uplo = "U")
	})
745
746	setMethod("fitted", signature(object = "lmer"),	setMethod("residuals", signature(object = "mer"),
747	function(object, ...) object@fitted)	function(object, ...)
748		.NotYetImplemented()
749		)
750
751	setMethod("residuals", signature(object = "lmer"),	setMethod("resid", signature(object = "mer"),
752	function(object, ...) object@residuals)	function(object, ...)
753		.NotYetImplemented()
754		)
755
756	setMethod("resid", signature(object = "lmer"),	setMethod("summary", signature(object = "mer"),
757	function(object, ...) do.call("residuals", c(list(object), list(...))))	function(object, ...)
758		.NotYetImplemented()
759		)
760
761	setMethod("coef", signature(object = "lmer"),	setMethod("update", signature(object = "mer"),
762	function(object, ...)	function(object, ...)
763	{	.NotYetImplemented()
764	fef <- data.frame(rbind(fixef(object)), check.names = FALSE)	)
	ref <- as(ranef(object), "list")
	names(ref) <- names(object@flist)
	val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
	for (i in seq(a = val)) {
	refi <- ref[[i]]
	row.names(val[[i]]) <- row.names(refi)
	if (!all(names(refi) %in% names(fef)))
	stop("unable to align random and fixed effects")
	val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi
	}
	new("lmer.coef", val)
	})
765
	setMethod("plot", signature(x = "lmer.coef"),
	function(x, y, ...)
	{
	varying <- unique(do.call("c",
	lapply(x, function(el)
	names(el)[sapply(el,
	function(col)
	any(col != col[1]))])))
	gf <- do.call("rbind", lapply(x, "[", j = varying))
	gf$.grp <- factor(rep(names(x), sapply(x, nrow)))
	switch(min(length(varying), 3),
	qqmath(eval(substitute(~ x | .grp,
	list(x = as.name(varying[1])))), gf, ...),
	xyplot(eval(substitute(y ~ x | .grp,
	list(y = as.name(varying[1]),
	x = as.name(varying[2])))), gf, ...),
	splom(~ gf | .grp, ...))
	})
	setMethod("plot", signature(x = "lmer.ranef"),
	function(x, y, ...)
	{
	lapply(x, function(x) {
	cn <- lapply(colnames(x), as.name)
	switch(min(ncol(x), 3),
	qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
	xyplot(eval(substitute(y ~ x, list(y = cn[[1]], x = cn[[2]]))),
	x, ...),
	splom(~ x, ...))
	})
	})
766
767	setMethod("with", signature(data = "lmer"),	simss <- function(fm0, fma, nsim)
768	function(data, expr, ...)	{
769	eval(substitute(substitute(expr, list(. = quote(data))),	ysim <- simulate(fm0, nsim)
770	append(data@flist, eval(data@call$data)),	cv <- list(analyticGradient = FALSE, msMaxIter = 200:200,
771	enclos = parent.frame())	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		}

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.689   changed lines   Added in v.1165

---

root@r-forge.r-project.org	ViewVC Help
Powered by ViewVC 1.0.0