Diff of /pkg/R/lmer.R

Parent Directory | Revision Log | Patch

	revision 981, Tue Oct 11 14:12:19 2005 UTC	revision 1165, Sat Jan 14 23:41:55 2006 UTC
#	Line 2	Line 2
2
3	## Some utilities	## Some utilities
4
	## Return the index into the packed lower triangle
	Lind <- function(i,j) {
	if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))
	((i - 1) * i)/2 + j
	}
5	## Return the pairs of expressions separated by vertical bars	## Return the pairs of expressions separated by vertical bars
6	findbars <- function(term)	findbars <- function(term)
7	{	{
#	Line 23	Line 17
17	## that are separated by vertical bars	## that are separated by vertical bars
18	nobars <- function(term)	nobars <- function(term)
19	{	{
20	# FIXME: is the is.name in the condition redundant?	if (!('|' %in% all.names(term))) return(term)
	#   A name won't satisfy the first condition.
	if (!('|' %in% all.names(term)) || is.name(term)) return(term)
21	if (is.call(term) && term[[1]] == as.name('|')) return(NULL)	if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
22	if (length(term) == 2) {	if (length(term) == 2) {
23	nb <- nobars(term[[2]])	nb <- nobars(term[[2]])
#	Line 60	Line 52
52
53	## Expand an expression with colons to the sum of the lhs	## Expand an expression with colons to the sum of the lhs
54	## and the current expression	## and the current expression
55	colExpand <- function(term)	colExpand <- function(term)
56	{	{
57	if (is.name(term) || is.numeric(term)) return(term)	if (is.name(term) || is.numeric(term)) return(term)
#	Line 116	Line 107
107	analyticHessian = as.logical(analyticHessian))	analyticHessian = as.logical(analyticHessian))
108	}	}
109
110		setMethod("coef", signature(object = "lmer"),
111		function(object, ...)
112		{
113		fef <- data.frame(rbind(object@fixed), check.names = FALSE)
114		ref <- as(ranef(object), "list")
115		names(ref) <- names(object@flist)
116		val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
117		for (i in seq(a = val)) {
118		refi <- ref[[i]]
119		row.names(val[[i]]) <- row.names(refi)
120		if (!all(names(refi) %in% names(fef)))
121		stop("unable to align random and fixed effects")
122		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"),	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(), start,	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 , ...)
	## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
182	{	{
183	## match and check parameters	## match and check parameters
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)	cv <- do.call("lmerControl", control)
186	## evaluate glm.fit, a generalized linear fit of fixed effects only
187		## Must evaluate the model frame first and then fit the glm using
188		## that frame.  Otherwise missing values in the grouping factors
189		## cause inconsistent numbers of observations.
190	mf <- match.call()	mf <- match.call()
191	m <- match(c("family", "data", "subset", "weights",	m <- match(c("family", "data", "subset", "weights",
192	"na.action", "offset"), names(mf), 0)	"na.action", "offset"), names(mf), 0)
193	mf <- mf[c(1, m)]	mf <- fe <- mf[c(1, m)]
194	frame.form <- subbars(formula) # substitute `+' for `|'	frame.form <- subbars(formula) # substitute `+' for `|'
195	fixed.form <- nobars(formula)  # remove any terms with `|'	fixed.form <- nobars(formula)  # remove any terms with `|'
196	if (inherits(fixed.form, "name")) # RHS is empty - use a constant	if (inherits(fixed.form, "name")) # RHS is empty - use a constant
197	fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))	fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
198	environment(fixed.form) <- environment(frame.form) <- environment(formula)	environment(fixed.form) <- environment(frame.form) <- environment(formula)
199	mf$formula <- fixed.form
200	mf$x <- mf$model <- mf$y <- TRUE	## evaluate a model frame for fixed and random effects
201	mf[[1]] <- as.name("glm")	mf$formula <- frame.form
202	glm.fit <- eval(mf, parent.frame())	mf$family <- NULL
203		mf$drop.unused.levels <- TRUE
204		mf[[1]] <- as.name("model.frame")
205		frm <- eval(mf, parent.frame())
206
207		## 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	x <- glm.fit$x
215	y <- as.double(glm.fit$y)	y <- as.double(glm.fit$y)
216	family <- glm.fit$family	family <- glm.fit$family
217
218	## check for a linear mixed model	## check for a linear mixed model
219	lmm <- family$family == "gaussian" && family$link == "identity"	lmm <- family$family == "gaussian" && family$link == "identity"
220	if (lmm) { # linear mixed model	if (lmm) { # linear mixed model
#	Line 165	Line 236
236	'\nUsing method = "PQL".\n')	'\nUsing method = "PQL".\n')
237	}	}
238	}	}
239		## create factor list for the random effects
	## evaluate a model frame for fixed and random effects
	mf$formula <- frame.form
	mf$x <- mf$model <- mf$y <- mf$family <- NULL
	mf$drop.unused.levels <- TRUE
	mf[[1]] <- as.name("model.frame")
	frm <- eval(mf, parent.frame())
	## grouping factors and model matrices for 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, function(x)	fl <- lapply(bars,
243	list(model.matrix(eval(substitute(~ T, list(T = x[[2]]))),	function(x)
	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		fl <- fl[ord]
252	## Create the model matrices and a mixed-effects representation (mer)	}
253	mmats <- c(lapply(random, "[[", 1),	## create list of transposed model matrices for random effects
254	.fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))	Ztl <- lapply(bars, function(x)
255	mer <- .Call("lmer_create", lapply(random, "[[", 2),	t(model.matrix(eval(substitute(~ expr,
256	mmats, method, PACKAGE = "Matrix")	list(expr = x[[2]]))),
257	if (lmm) {                    ## linear mixed model	frm)))
258	if (missing(start)) .Call("lmer_initial", mer, PACKAGE="Matrix")	## Create the mixed-effects representation (mer) object
259	else .Call("lmer_set_initial", mer, start, PACKAGE = "Matrix")	mer <- .Call("mer_create", fl,
260	.Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")	.Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"),
261		x, y, method, sapply(Ztl, nrow),
262		c(lapply(Ztl, rownames), list(.fixed = colnames(x))),
263		!(family$family %in% c("binomial", "poisson")),
264		match.call(), family,
265		PACKAGE = "Matrix")
266		if (lmm) {
267		.Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
268	LMEoptimize(mer) <- cv	LMEoptimize(mer) <- cv
269	fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")	return(mer)
	return(new("lmer",
	mer,
	assign = attr(x, "assign"),
	call = match.call(),
	family = family, fitted = fits,
	fixed = fixef(mer),
	frame = if (model) frm else data.frame(),
	logLik = logLik(mer),
	residuals = unname(model.response(frm) - fits),
	terms = glm.fit$terms))
270	}	}
271
272	## The rest of the function applies to generalized linear mixed models	## The rest of the function applies to generalized linear mixed models
#	Line 219	Line 276
276	wtssqr <- wts * wts	wtssqr <- wts * wts
277	offset <- glm.fit$offset	offset <- glm.fit$offset
278	if (is.null(offset)) offset <- numeric(length(eta))	if (is.null(offset)) offset <- numeric(length(eta))
	mu <- numeric(length(eta))
	dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
279	linkinv <- quote(family$linkinv(eta))	linkinv <- quote(family$linkinv(eta))
280	mu.eta <- quote(family$mu.eta(eta))	mu.eta <- quote(family$mu.eta(eta))
281		mu <- family$linkinv(eta)
282	variance <- quote(family$variance(mu))	variance <- quote(family$variance(mu))
283		dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
284	LMEopt <- get("LMEoptimize<-")	LMEopt <- get("LMEoptimize<-")
285	doLMEopt <- quote(LMEopt(x = mer, value = cv))	doLMEopt <- quote(LMEopt(x = mer, value = cv))
286
287	GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")	GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
288	.Call("glmer_PQL", GSpt, PACKAGE = "Matrix")  # obtain PQL estimates	.Call("glmer_PQL", GSpt, PACKAGE = "Matrix")  # obtain PQL estimates
289	fixInd <- seq(ncol(x))	fixInd <- seq(ncol(x))
290	## pars[fixInd] == beta, pars[-fixInd] == theta	## pars[fixInd] == beta, pars[-fixInd] == theta
291	PQLpars <- c(fixef(mer),	PQLpars <- c(fixef(mer),
292	.Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))	.Call("mer_coef", mer, 2, PACKAGE = "Matrix"))
293	## set flag to skip fixed-effects in subsequent calls	.Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")
	mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
294	## indicator of constrained parameters	## indicator of constrained parameters
295	const <- c(rep(FALSE, length(fixInd)),	const <- c(rep(FALSE, length(fixInd)),
296	unlist(lapply(mer@nc[seq(along = random)],	unlist(lapply(mer@nc[seq(along = fl)],
297	function(k) 1:((k*(k+1))/2) <= k)	function(k) 1:((k*(k+1))/2) <= k)
298	))	))
299	devAGQ <- function(pars, n)	devLaplace <- function(pars)
300	.Call("glmer_devAGQ", pars, GSpt, n, PACKAGE = "Matrix")	.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)	deviance <- devAGQ(PQLpars, 1)
310
311	### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs	### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs
312	### AGQ for nc > 1 first.	### AGQ for nc > 1 first.
313	fxd <- PQLpars[fixInd]	fxd <- PQLpars[fixInd]
#	Line 265	Line 327
327	}	}
328	obj <- function(pars)	obj <- function(pars)
329	.Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")	.Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
	if (exists("nlminb", mode = "function")) {
330	optimRes <-	optimRes <-
331	nlminb(PQLpars, obj,	nlminb(PQLpars, obj,
332	lower = ifelse(const, 5e-10, -Inf),	lower = ifelse(const, 5e-10, -Inf),
#	Line 275	Line 336
336	if (optimRes$convergence != 0)	if (optimRes$convergence != 0)
337	warning("nlminb failed to converge")	warning("nlminb failed to converge")
338	deviance <- optimRes$objective	deviance <- optimRes$objective
339	} else {	if (gVerb)
	optimRes <-
	optim(PQLpars, obj, method = "L-BFGS-B",
	lower = ifelse(const, 5e-10, -Inf),
	control = list(trace = getOption("verbose"),
	maxit = cv$msMaxIter))
	optpars <- optimRes$par
	if (optimRes$convergence != 0)
	warning("optim failed to converge")
	deviance <- optimRes$value
	}
	if (gVerb) {
340	cat(paste("convergence message", optimRes$message, "\n"))	cat(paste("convergence message", optimRes$message, "\n"))
	}
341	fxd[] <- optpars[fixInd]  ## preserve the names	fxd[] <- optpars[fixInd]  ## preserve the names
342	.Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")	.Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
343	}	}
#	Line 302	Line 351
351	call = match.call(), family = family,	call = match.call(), family = family,
352	logLik = loglik, fixed = fxd)	logLik = loglik, fixed = fxd)
353	})	})
	## end{  "lmer . formula " }
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		## Optimization for mer objects
376	setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),	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)
380	nc <- x@nc	nc <- x@nc
381	constr <- unlist(lapply(nc[1:(length(nc) - 2)],	constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
	function(k) 1:((k*(k+1))/2) <= k))
382	fn <- function(pars)	fn <- function(pars)
383	deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))	deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
384	gr <-	gr <- if (value$analyticGradient)
	if (value$analyticGradient)
385	function(pars) {	function(pars) {
386	if (!isTRUE(all.equal(pars,	if (!isTRUE(all.equal(pars,
387	.Call("lmer_coef", x,	.Call("mer_coef", x,
388	2, PACKAGE = "Matrix"))))	2, PACKAGE = "Matrix"))))
389	.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")	.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")
390	.Call("lmer_gradient", x, 2, PACKAGE = "Matrix")	.Call("mer_gradient", x, 2, PACKAGE = "Matrix")
391	}	}
392	## else NULL	else NULL
393	optimRes <-	optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"),
	if (exists("nlminb", mode = "function"))
	nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
394	fn, gr,	fn, gr,
395	lower = ifelse(constr, 5e-10, -Inf),	lower = ifelse(constr, 5e-10, -Inf),
396	control = list(iter.max = value$msMaxIter,	control = list(iter.max = value$msMaxIter,
397	trace = as.integer(value$msVerbose)))	trace = as.integer(value$msVerbose)))
398	else	.Call("mer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
	optim(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
	fn, gr, method = "L-BFGS-B",
	lower = ifelse(constr, 5e-10, -Inf),
	control = list(maxit = value$msMaxIter,
	trace = as.integer(value$msVerbose)))
	.Call("lmer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
399	if (optimRes$convergence != 0) {	if (optimRes$convergence != 0) {
400	warning(paste("optim or nlminb returned message",	warning(paste("nlminb returned message",
401	optimRes$message,"\n"))	optimRes$message,"\n"))
402	}	}
403	return(x)	return(x)
404	})	})
405
406	setMethod("ranef", signature(object = "lmer"),	setMethod("deviance", signature(object = "mer"),
407	function(object, accumulate = FALSE, ...) {	function(object, ...) {
408	val <- new("lmer.ranef",	.Call("mer_factor", object, PACKAGE = "Matrix")
409	lapply(.Call("lmer_ranef", object, PACKAGE = "Matrix"),	object@deviance[[ifelse(object@method == "REML", "REML", "ML")]]
	data.frame, check.names = FALSE),
	varFac = object@bVar,
	stdErr = .Call("lmer_sigma", object,
	object@method == "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
410	})	})
411
412	setMethod("fixef", signature(object = "mer"),	setMethod("mcmcsamp", signature(object = "mer"),
413	function(object, ...)	function(object, n = 1, verbose = FALSE, saveb = FALSE,
414	.Call("lmer_fixef", object, PACKAGE = "Matrix"))	trans = TRUE, ...)
415		{
416		ans <- t(.Call("mer_MCMCsamp", object, saveb, n,
417		trans, PACKAGE = "Matrix"))
418		attr(ans, "mcpar") <- as.integer(c(1, n, 1))
419		class(ans) <- "mcmc"
420		glmer <- FALSE
421		gnms <- names(object@flist)
422		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("fixef", signature(object = "lmer"),	setMethod("simulate", signature(object = "mer"),
448	function(object, ...) object@fixed)	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	setMethod("VarCorr", signature(x = "lmer"),	family <- object@family
462	##FIXME - change this for reasonable defaults of useScale according to	if (family$family != "gaussian" ||
463	##the family slot.	family$link != "identity")
464	function(x, REML = TRUE, useScale = TRUE, ...) {	stop("simulation of generalized linear mixed models not yet implemented")
465	val <- .Call("lmer_variances", x, PACKAGE = "Matrix")	## similate the linear predictors
466	for (i in seq(along = val)) {	lpred <- .Call("mer_simulate", object, nsim, PACKAGE = "Matrix")
467	dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])	sc <- 1
468	val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")	if (object@useScale)
469	}	sc <- .Call("mer_sigma", object, object@method == "REML",
470	new("VarCorr",	PACKAGE = "Matrix")
471	scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),	## add fixed-effects contribution and per-observation noise term
472	reSumry = val,	lpred <- as.data.frame(lpred + drop(object@X %*% fixef(object)) +
473	useScale = useScale)	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, unconst, ...)
	.Call("lmer_gradient", x, unconst, PACKAGE = "Matrix"))
	setMethod("summary", signature(object = "lmer"),
	function(object, ...)
	new("summary.lmer", object,
	showCorrelation = TRUE,
	useScale = !((object@family)$family %in% c("binomial", "poisson"))))
	setMethod("show", signature(object = "lmer"),
	function(object)
	show(new("summary.lmer", object,
	showCorrelation = FALSE,
	useScale = !((object@family)$family %in% c("binomial", "poisson")))))
479
480	setMethod("show", "summary.lmer",	setMethod("show", "mer",
481	function(object) {	function(object) {
482	fcoef <- object@fixed	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 <- object@method == "REML"	REML <- object@method == "REML"
525	llik <- object@logLik	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"))) {	if (glz <- !(object@method %in% c("REML", "ML"))) {
#	Line 447	Line 559
559	row.names = ""))	row.names = ""))
560	}	}
561	cat("Random effects:\n")	cat("Random effects:\n")
562	show(VarCorr(object, useScale = useScale))	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, FALSE, 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 473	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("vcov", signature(object = "mer"),
607		function(object, REML = object@method == "REML",
608		useScale = object@useScale,...) {
609		sc <- if (object@useScale) {
610		.Call("mer_sigma", object, REML, PACKAGE = "Matrix")
611		} else { 1 }
612		rr <- as(sc^2 * tcrossprod(solve(object@RXX)), "dpoMatrix")
613		rr@factors$correlation <- as(rr, "correlation")
614		rr
615		})
616
617	## calculates degrees of freedom for fixed effects Wald tests	## calculates degrees of freedom for fixed effects Wald tests
618	## This is a placeholder.  The answers are generally wrong.  It will	## This is a placeholder.  The answers are generally wrong.  It will
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 <- abs(nc[1]) - 1
	n <- nc[2]
	rep(n-p, p)
626	})	})
627
628	setMethod("logLik", signature(object="mer"),	setMethod("logLik", signature(object="mer"),
629	function(object, REML = object@method == "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") <- abs(nc[1]) +	attr(val, "df") <- abs(devc[2]) +
634	length(.Call("lmer_coef", object, 0, PACKAGE = "Matrix"))	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("logLik", signature(object="lmer"),	setMethod("VarCorr", signature(x = "mer"),
641	function(object, ...) object@logLik)	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 = "lmer"),	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 566	Line 698
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 580	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 595	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("deviance", "mer",
	function(object, REML = NULL, ...) {
	.Call("lmer_factor", object, PACKAGE = "Matrix")
	if (is.null(REML))
	REML <- object@method == "REML"
	object@deviance[[ifelse(REML, "REML", "ML")]]
	})
	setMethod("deviance", "lmer",
	function(object, ...) -2 * c(object@logLik))
	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@method == "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
	})
	## Extract the L matrix
	setAs("lmer", "dtTMatrix",
	function(from)
	{
	## 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")
	})
740
741	## Extract the ZZX matrix	setMethod("formula", signature(x = "mer"),
742	setAs("lmer", "dsTMatrix",	function(x, ...)
743	function(from)	x@call$formula
744	{	)
	.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, ...)	function(object, ...)
748	napredict(attr(object@frame, "na.action"), object@fitted))	.NotYetImplemented()
749		)
750
751	setMethod("residuals", signature(object = "lmer"),	setMethod("resid", signature(object = "mer"),
752	function(object, ...)	function(object, ...)
753	naresid(attr(object@frame, "na.action"), object@residuals))	.NotYetImplemented()
754		)
	setMethod("resid", signature(object = "lmer"),
	function(object, ...) do.call("residuals", c(list(object), list(...))))
755
756	setMethod("coef", signature(object = "lmer"),	setMethod("summary", signature(object = "mer"),
757	function(object, ...)	function(object, ...)
758	{	.NotYetImplemented()
759	fef <- data.frame(rbind(object@fixed), 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)
	})
	setMethod("plot", signature(x = "lmer.coef"),
	function(x, y, ...)
	{
	## require("lattice", quietly = TRUE) -- now via Imports
	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, ...)
	{
	## require("lattice", quietly = TRUE) -- now via Imports
	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, ...))
	})
	})
	setMethod("with", signature(data = "lmer"),
	function(data, expr, ...) {
	dat <- eval(data@call$data)
	if (!is.null(na.act <- attr(data@frame, "na.action")))
	dat <- dat[-na.act, ]
	lst <- c(list(. = data), data@flist, data@frame, dat)
	eval(substitute(expr), lst[unique(names(lst))])
	})
	setMethod("terms", signature(x = "lmer"),
	function(x, ...) x@terms)
	setMethod("show", signature(object="VarCorr"),
	function(object)
	{
	digits <- max(3, getOption("digits") - 2)
	useScale <- length(object@useScale) > 0 && object@useScale[1]
	sc <- ifelse(useScale, object@scale,  1.)
	reStdDev <- c(lapply(object@reSumry,
	function(x, sc)
	sc*x@stdDev,
	sc = sc), list(Residual = sc))
	reLens <- unlist(c(lapply(reStdDev, length)))
	reMat <- array('', c(sum(reLens), 4),
	list(rep('', sum(reLens)),
	c("Groups", "Name", "Variance", "Std.Dev.")))
	reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
	reMat[,2] <- c(unlist(lapply(reStdDev, names)), "")
	reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
	reMat[,4] <- format(unlist(reStdDev), digits = digits)
	if (any(reLens > 1)) {
	maxlen <- max(reLens)
	corr <-
	do.call("rbind",
	lapply(object@reSumry,
	function(x, maxlen) {
	cc <- format(round(x, 3), nsmall = 3)
	cc[!lower.tri(cc)] <- ""
	nr <- dim(cc)[1]
	if (nr >= maxlen) return(cc)
	cbind(cc, matrix("", nr, maxlen-nr))
	}, maxlen))
	colnames(corr) <- c("Corr", rep("", maxlen - 1))
	reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr))))
	}
	if (!useScale) reMat <- reMat[-nrow(reMat),]
	print(reMat, quote = FALSE)
	})
	setMethod("mcmcsamp", signature(object = "lmer"),
	function(object, n = 1, verbose = FALSE, saveb = FALSE,
	trans = TRUE, ...)
	{
	if (object@family$family == "gaussian" &&
	object@family$link == "identity") {
	glmer <- FALSE
	ans <- .Call("lmer_MCMCsamp", object, saveb, n, trans,
	PACKAGE = "Matrix")
	} else {
	glmer <- TRUE
	if (trans)
	warning("trans option not currently allowed for generalized models")
	trans <- FALSE
	## Check arguments
	if (length(object@Omega) > 1 || object@nc[1] > 1)
	stop("mcmcsamp currently defined for glmm models with only one variance component")
	cv <- Matrix:::lmerControl()
	if (verbose) cv$msVerbose <- 1
	family <- object@family
	frm <- object@frame
	## recreate model matrices
	fixed.form <- Matrix:::nobars(object@call$formula)
	if (inherits(fixed.form, "name")) # RHS is empty - use a constant
	fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
	glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE,
	y = TRUE)
	x <- glm.fit$x
	y <- as.double(glm.fit$y)
	bars <- Matrix:::findbars(object@call$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]])))
	if (any(names(random) != names(object@flist)))
	random <- random[names(object@flist)]
	mmats <- c(lapply(random, "[[", 1),
	.fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
	mer <- as(object, "mer")
	## establish the GS object and the ans matrix
	eta <- glm.fit$linear.predictors # perhaps later change this to object@fitted?
	wts <- glm.fit$prior.weights
	wtssqr <- wts * wts
	offset <- glm.fit$offset
	if (is.null(offset)) offset <- numeric(length(eta))
	off <- numeric(length(eta))
	mu <- numeric(length(eta))
	dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
	linkinv <- quote(family$linkinv(eta))
	mu.eta <- quote(family$mu.eta(eta))
	variance <- quote(family$variance(mu))
	LMEopt <- getAnywhere("LMEoptimize<-")
	doLMEopt <- quote(LMEopt(x = mer, value = cv))
	GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
	fixed <- object@fixed
	varc <- .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")
	b <- .Call("lmer_ranef", mer, PACKAGE = "Matrix")
	ans <- .Call("glmer_MCMCsamp", GSpt, b, fixed, varc, saveb, n,
	PACKAGE = "Matrix")
	.Call("glmer_finalize", GSpt, PACKAGE = "Matrix");
	}
	gnms <- names(object@flist)
	cnms <- object@cnames
	ff <- fixef(object)
	colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
	unlist(lapply(seq(along = gnms),
	function(i)
	abbrvNms(gnms[i],cnms[[i]]))))
	if (trans) {
	## parameter type: 0 => fixed effect, 1 => variance,
	##                 2 => covariance
	ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1,
	unlist(lapply(seq(along = gnms),
	function(i)
	{
	k <- length(cnms[[i]])
	rep(1:2, c(k, (k*(k-1))/2))
	})))
	colnms[ptyp == 1] <-
	paste("log(", colnms[ptyp == 1], ")", sep = "")
	colnms[ptyp == 2] <-
	paste("atanh(", colnms[ptyp == 2], ")", sep = "")
	}
	colnames(ans) <- colnms
	ans
	})
	rWishart <- function(n, df, invScal)
	.Call("Matrix_rWishart", n, df, invScal)
760
761	setMethod("model.matrix", signature(object = "lmer"),	setMethod("update", signature(object = "mer"),
762	function(object, ...)	function(object, ...)
763	{	.NotYetImplemented()
764	frm <- object@frame	)
	fixed.form <- Matrix:::nobars(object@call$formula)
	if (inherits(fixed.form, "name")) # RHS is empty - use a constant
	fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
	glm.fit <- glm(eval(fixed.form), object@family, frm, x = TRUE, y = TRUE)
	fxd <- unname(drop(glm.fit$x %*% fixef(object)))
	## Create the random effects model matrices
	bars <- Matrix:::findbars(object@call$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]])))
	## re-order the random effects pairs if necessary
	if (any(names(random) != names(object@flist)))
	random <- random[names(object@flist)]
	c(lapply(random, "[[", 1),
	.fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
	})
	setMethod("simulate", signature(object = "lmer"),
	function(object, nsim = 1, seed = NULL, ...)
	{
	if(!exists(".Random.seed", envir = .GlobalEnv))
	runif(1)               # initialize the RNG if necessary
	if(is.null(seed))
	RNGstate <- .Random.seed
	else {
	R.seed <- .Random.seed
	set.seed(seed)
	RNGstate <- structure(seed, kind = as.list(RNGkind()))
	on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
	}
	family <- object@family
	if (family$family != "gaussian" ||
	family$link != "identity")
	stop("simulation of generalized linear mixed models not yet implemented")
	## pieces we will need later
	scale <- .Call("lmer_sigma", object, object@method == "REML",
	PACKAGE = "Matrix")
	mmats <- model.matrix(object)
	ff <- fixef(object)
	###_FIXME: If the factor levels have been permuted, has the
	###        permutation been applied in the stored frame?  Otherwise we
	###        need to check this.
	## similate the linear predictors
	lpred <- .Call("lmer_simulate", as(object, "mer"), nsim,
	unname(drop(mmats[[length(mmats)]][,seq(a = ff)] %*% ff)),
	mmats, TRUE, PACKAGE = "Matrix")
	## add per-observation noise term
	lpred <- as.data.frame(lpred + rnorm(prod(dim(lpred)), sd = scale))
765
	## save the seed
	attr(lpred, "seed") <- RNGstate
	lpred
	})
766
767	simulate2 <- function(object, n = 1, ...)	simss <- function(fm0, fma, nsim)
768	{	{
769	family <- object@family	ysim <- simulate(fm0, nsim)
770	if (family$family != "gaussian" ||	cv <- list(analyticGradient = FALSE, msMaxIter = 200:200,
771	family$link != "identity")	msVerbose = 0:0)
772	stop("simulation of generalized linear mixed models not implemented yet")	sapply(ysim, function(yy) {
773		.Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")
774	## create the mean from the fixed effects	LMEoptimize(fm0) <- cv
775	frm <- object@frame	.Call("mer_update_y", fma, yy, PACKAGE = "Matrix")
776	fixed.form <- Matrix:::nobars(object@call$formula)	LMEoptimize(fma) <- cv
777	if (inherits(fixed.form, "name")) # RHS is empty - use a constant	exp(c(H0 = fm0@devComp[["logryy2"]],
778	fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))	Ha = fma@devComp[["logryy2"]]))
	glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE, y = TRUE)
	lpred <- matrix(glm.fit$x %*% fixef(object), nr = nrow(frm), nc = n)
	## Create the random effects model matrices
	bars <- Matrix:::findbars(object@call$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]])))
	## re-order the random effects pairs if necessary
	flist <- object@flist
	if (any(names(random) != names(flist)))
	random <- random[names(flist)]
	mmats <- lapply(random, "[[", 1)
	## simulate the random effects
	scale <- .Call("lmer_sigma", object, object@method == "REML",
	PACKAGE = "Matrix")
	Omega <- object@Omega
	re <- lapply(seq(along = Omega),
	function(i) {
	om <- Omega[[i]]
	nr <- nrow(om)
	nlev <- length(levels(flist[[i]]))
	scale * array(solve(chol(new("dpoMatrix", Dim = dim(om),
	uplo = "U", x = c(om))),
	matrix(rnorm(nr * n * nlev),
	nr = nr))@x, c(nr, n, nlev))
779	})	})
	## apply the random effects
	for (j in seq(along = Omega)) {
	for (i in 1:nrow(lpred))
	lpred[i,] <- lpred[i,] + mmats[[j]][i,] %*% re[[j]][, , as.integer(flist[[j]])[i]]
	}
	## add per-observation noise term
	lpred <- lpred + rnorm(prod(dim(lpred)), sd = scale)
	attr(lpred, "re") <- re
	lpred
	}
	refdist <- function(fm1, fm2, n, ...)
	{
	cv <- lmerControl()
	obs <- deviance(fm2) - deviance(fm1)
	newy <- simulate(fm2, n)
	mm1 <- model.matrix(fm1)
	mm2 <- model.matrix(fm2)
	ref <- numeric(n)
	mer1 <- as(fm1, "mer")
	mer2 <- as(fm2, "mer")
	for (j in 1:n) {
	.Call("lmer_update_y", mer2, newy[[j]], mm2, PACKAGE = "Matrix")
	LMEoptimize(mer2) <- cv
	.Call("lmer_update_y", mer1, newy[[j]], mm1, PACKAGE = "Matrix")
	LMEoptimize(mer1) <- cv
	ref[j] <- deviance(mer2) - deviance(mer1)
780	}	}
	attr(ref, "observed") <- obs
	ref
	}
	mer2 <-
	function(formula, data, family,
	method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
	control = list(), start,
	subset, weights, na.action, offset,
	model = TRUE, x = FALSE, y = FALSE , ...)
	## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
	{
	## match and check parameters
	if (length(formula) < 3) stop("formula must be a two-sided formula")
	## cv <- do.call("Matrix:::lmerControl", control)
	cv <- control
	cv$analyticGradient <- FALSE
	cv$msMaxIter <- as.integer(200)
	if (is.null(cv$msVerbose)) cv$msVerbose <- as.integer(1)
	## Must evaluate the model frame first and then fit the glm using
	## that frame.  Otherwise missing values in the grouping factors
	## cause inconsistent numbers of observations.
	## evaluate glm.fit, a generalized linear fit of fixed effects only
	mf <- match.call()
	m <- match(c("family", "data", "subset", "weights",
	"na.action", "offset"), names(mf), 0)
	mf <- fe <- mf[c(1, m)]
	frame.form <- subbars(formula) # substitute `+' for `|'
	fixed.form <- nobars(formula)  # remove any terms with `|'
	if (inherits(fixed.form, "name")) # RHS is empty - use a constant
	fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
	environment(fixed.form) <- environment(frame.form) <- environment(formula)
	## evaluate a model frame for fixed and random effects
	mf$formula <- frame.form
	mf$family <- NULL
	mf$drop.unused.levels <- TRUE
	mf[[1]] <- as.name("model.frame")
	frm <- eval(mf, parent.frame())
	## fit a glm model to the fixed formula
	fe$formula <- fixed.form
	fe$data <- frm
	fe$x <- fe$model <- fe$y <- TRUE
	fe[[1]] <- as.name("glm")
	glm.fit <- eval(fe, parent.frame())
	x <- glm.fit$x
	y <- as.double(glm.fit$y)
	family <- glm.fit$family
	## check for a linear mixed model
	lmm <- family$family == "gaussian" && family$link == "identity"
	if (lmm) { # linear mixed model
	method <- match.arg(method)
	if (method %in% c("PQL", "Laplace", "AGQ")) {
	warning(paste('Argument method = "', method,
	'" is not meaningful for a linear mixed model.\n',
	'Using method = "REML".\n', sep = ''))
	method <- "REML"
	}
	} else { # generalized linear mixed model
	if (missing(method)) method <- "PQL"
	else {
	method <- match.arg(method)
	if (method == "ML") method <- "PQL"
	if (method == "REML")
	warning('Argument method = "REML" is not meaningful ',
	'for a generalized linear mixed model.',
	'\nUsing method = "PQL".\n')
	}
	}
	## grouping factors and model matrices for random effects
	bars <- findbars(formula[[3]])
	random <-
	lapply(bars, function(x)
	list(t(model.matrix(eval(substitute(~ expr,
	list(expr = 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))]
	}
	## Create the model matrices and a mixed-effects representation (mer)
	mer <- .Call("mer2_create", random, x, y, method, PACKAGE = "Matrix")
	LMEoptimize(mer) <- cv
	mer
	}
	## Extract the permutation
	setAs("mer2", "pMatrix", function(from) .Call("mer2_pMatrix", from))
	## Extract the L matrix
	setAs("mer2", "dtCMatrix", function(from) .Call("mer2_dtCMatrix", from))
	## Optimization for mer2 objects
	setReplaceMethod("LMEoptimize", signature(x="mer2", value="list"),
	function(x, value)
	{
	if (value$msMaxIter < 1) return(x)
	nc <- x@nc
	constr <- unlist(lapply(nc[1:(length(nc) - 2)],
	function(k) 1:((k*(k+1))/2) <= k))
	fn <- function(pars)
	deviance(.Call("mer2_coefGets", x, pars, 2, PACKAGE = "Matrix"))
	gr <- NULL  ## No gradient yet
	##                     if (value$analyticGradient)
	##                         function(pars) {
	##                             if (!isTRUE(all.equal(pars,
	##                                                   .Call("lmer_coef", x,
	##                                                         2, PACKAGE = "Matrix"))))
	##                                 .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
	##                             .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
	##                         }
	## else NULL
	optimRes <- nlminb(.Call("mer2_coef", x, 2, PACKAGE = "Matrix"),
	fn, gr,
	lower = ifelse(constr, 5e-10, -Inf),
	control = list(iter.max = value$msMaxIter,
	trace = as.integer(value$msVerbose)))
	.Call("mer2_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
	if (optimRes$convergence != 0) {
	warning(paste("nlminb returned message",
	optimRes$message,"\n"))
	}
	return(x)
	})
	setMethod("deviance", "mer2",
	function(object, REML = NULL, ...) {
	.Call("mer2_factor", object, PACKAGE = "Matrix")
	if (is.null(REML))
	REML <- object@method == "REML"
	object@deviance[[ifelse(REML, "REML", "ML")]]
	})

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

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

---

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