Diff of /pkg/R/lmer.R

Parent Directory | Revision Log | Patch

	revision 831, Mon Aug 8 08:24:10 2005 UTC	revision 832, Mon Aug 8 08:33:55 2005 UTC
#	Line 51	Line 51
51	return(term)	return(term)
52	}	}
53	stopifnot(length(term) == 3)	stopifnot(length(term) == 3)
54	if (is.call(term) && term[[1]] == as.name('|')) term[[1]] <- as.name('+')	if (is.call(term) && term[[1]] == as.name('|'))
55		term[[1]] <- as.name('+')
56	term[[2]] <- subbars(term[[2]])	term[[2]] <- subbars(term[[2]])
57	term[[3]] <- subbars(term[[3]])	term[[3]] <- subbars(term[[3]])
58	term	term
#	Line 71	Line 72
72
73	## Control parameters for lmer	## Control parameters for lmer
74	lmerControl <-	lmerControl <-
75	function(maxIter = 200,	function(maxIter = 200, # used in ../src/lmer.c only
76		tolerance = sqrt(.Machine$double.eps),# dito
77	msMaxIter = 200,	msMaxIter = 200,
78	tolerance = sqrt((.Machine$double.eps)),	## msTol = sqrt(.Machine$double.eps),
79	niterEM = 15,	## FIXME:  should be able to pass tolerances to nlminb()
	msTol = sqrt(.Machine$double.eps),
80	msVerbose = getOption("verbose"),	msVerbose = getOption("verbose"),
81	PQLmaxIt = 30,	niterEM = 15,
82	EMverbose = getOption("verbose"),	EMverbose = getOption("verbose"),
83		PQLmaxIt = 30,# unused _FIXME_
84	analyticGradient = TRUE,	analyticGradient = TRUE,
85	analyticHessian = FALSE)	analyticHessian = FALSE # unused _FIXME_
86		)
87	{	{
88	list(maxIter = as.integer(maxIter),	list(maxIter = as.integer(maxIter),
	msMaxIter = as.integer(msMaxIter),
89	tolerance = as.double(tolerance),	tolerance = as.double(tolerance),
90		msMaxIter = as.integer(msMaxIter),
91		## msTol = as.double(msTol),
92		msVerbose = as.integer(msVerbose),# "integer" on purpose
93	niterEM = as.integer(niterEM),	niterEM = as.integer(niterEM),
	msTol = as.double(msTol),
	msVerbose = as.logical(msVerbose),
	PQLmaxIt = as.integer(PQLmaxIt),
94	EMverbose = as.logical(EMverbose),	EMverbose = as.logical(EMverbose),
95		PQLmaxIt = as.integer(PQLmaxIt),
96	analyticGradient = as.logical(analyticGradient),	analyticGradient = as.logical(analyticGradient),
97	analyticHessian = as.logical(analyticHessian))	analyticHessian = as.logical(analyticHessian))
98	}	}
#	Line 100	Line 103
103	control = list(),	control = list(),
104	subset, weights, na.action, offset,	subset, weights, na.action, offset,
105	model = TRUE, x = FALSE, y = FALSE, ...)	model = TRUE, x = FALSE, y = FALSE, ...)
106	{                                 ## match and check parameters	## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
107		{
108		## match and check parameters
109	if (length(formula) < 3) stop("formula must be a two-sided formula")	if (length(formula) < 3) stop("formula must be a two-sided formula")
110	cv <- do.call("lmerControl", control)	cv <- do.call("lmerControl", control)
	if (lmm <- missing(family)) { # 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
	method <- if (missing(method)) "PQL" else match.arg(method)
	if (method == "ML") method <- "PQL"
	if (method == "REML")
	warning(paste('Argument method = "REML" is not meaningful',
	'for a generalized linear mixed model.',
	'\nUsing method = "PQL".\n'))
	}
111	## evaluate glm.fit, a generalized linear fit of fixed effects only	## evaluate glm.fit, a generalized linear fit of fixed effects only
112	mf <- match.call()	mf <- match.call()
113	m <- match(c("family", "data", "subset", "weights",	m <- match(c("family", "data", "subset", "weights",
#	Line 136	Line 124
124	glm.fit <- eval(mf, parent.frame())	glm.fit <- eval(mf, parent.frame())
125	x <- glm.fit$x	x <- glm.fit$x
126	y <- as.double(glm.fit$y)	y <- as.double(glm.fit$y)
127		lmm <- missing(family) # => linear mixed model
128	family <- glm.fit$family	family <- glm.fit$family
129
130		if (lmm) { # linear mixed model
131		method <- match.arg(method)
132		if (method %in% c("PQL", "Laplace", "AGQ")) {
133		warning(paste('Argument method = "', method,
134		'" is not meaningful for a linear mixed model.\n',
135		'Using method = "REML".\n', sep = ''))
136		method <- "REML"
137		}
138		} else { # generalized linear mixed model
139		if (missing(method)) method <- "PQL"
140		else {
141		method <- match.arg(method)
142
143		## <FIXME>: only
144		## if( !(family == "gaussian" && link == "identity") ) {
145		if (method == "ML") method <- "PQL"
146		if (method == "REML")
147		warning('Argument method = "REML" is not meaningful ',
148		'for a generalized linear mixed model.',
149		'\nUsing method = "PQL".\n')
150		## } </FIXME>
151		}
152		}
153
154	## evaluate a model frame for fixed and random effects	## evaluate a model frame for fixed and random effects
155	mf$formula <- frame.form	mf$formula <- frame.form
156	mf$x <- mf$model <- mf$y <- mf$family <- NULL	mf$x <- mf$model <- mf$y <- mf$family <- NULL
#	Line 148	Line 161
161	## grouping factors and model matrices for random effects	## grouping factors and model matrices for random effects
162	bars <- findbars(formula[[3]])	bars <- findbars(formula[[3]])
163	random <-	random <-
164	lapply(bars,	lapply(bars, function(x)
165	function(x) list(model.matrix(eval(substitute(~term,	list(model.matrix(eval(substitute(~ T, list(T = x[[2]]))),
	list(term=x[[2]]))),
166	frm),	frm),
167	eval(substitute(as.factor(fac)[,drop = TRUE],	eval(substitute(as.factor(fac)[,drop = TRUE],
168	list(fac = x[[3]])), frm)))	list(fac = x[[3]])), frm)))
#	Line 267	Line 279
279
280	.Call("glmer_finalize", GSpt, PACKAGE = "Matrix")	.Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
281	loglik[] <- -deviance/2	loglik[] <- -deviance/2
282	new("lmer", mer, frame = frm, terms = glm.fit$terms,	new("lmer", mer,
283		frame = if (model) frm else data.frame(),
284		terms = glm.fit$terms,
285	assign = attr(glm.fit$x, "assign"),	assign = attr(glm.fit$x, "assign"),
286	call = match.call(), family = family,	call = match.call(), family = family,
287	logLik = loglik, fixed = fxd)	logLik = loglik, fixed = fxd)
288	})	})
289		## end{  "lmer . formula " }
290
291	setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),	setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
292	function(x, value)	function(x, value)
#	Line 282	Line 297
297	function(k) 1:((k*(k+1))/2) <= k))	function(k) 1:((k*(k+1))/2) <= k))
298	fn <- function(pars)	fn <- function(pars)
299	deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))	deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
	gr <- NULL
	if (value$analyticGradient)
300	gr <-	gr <-
301		if (value$analyticGradient)
302	function(pars) {	function(pars) {
303	if (!isTRUE(all.equal(pars,	if (!isTRUE(all.equal(pars,
304	.Call("lmer_coef", x,	.Call("lmer_coef", x,
#	Line 292	Line 306
306	.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")	.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
307	.Call("lmer_gradient", x, 2, PACKAGE = "Matrix")	.Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
308	}	}
309		## else NULL
310	optimRes <-	optimRes <-
311	if (exists("nlminb", mode = "function"))	if (exists("nlminb", mode = "function"))
312	nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),	nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
#	Line 737	Line 752
752	setMethod("plot", signature(x = "lmer.coef"),	setMethod("plot", signature(x = "lmer.coef"),
753	function(x, y, ...)	function(x, y, ...)
754	{	{
755	if (require("lattice", quietly = TRUE)) {	## require("lattice", quietly = TRUE) -- now via Imports
756	varying <- unique(do.call("c",	varying <- unique(do.call("c",
757	lapply(x, function(el)	lapply(x, function(el)
758	names(el)[sapply(el,	names(el)[sapply(el,
#	Line 752	Line 767
767	list(y = as.name(varying[1]),	list(y = as.name(varying[1]),
768	x = as.name(varying[2])))), gf, ...),	x = as.name(varying[2])))), gf, ...),
769	splom(~ gf | .grp, ...))	splom(~ gf | .grp, ...))
	}
770	})	})
771
772	setMethod("plot", signature(x = "lmer.ranef"),	setMethod("plot", signature(x = "lmer.ranef"),
773	function(x, y, ...)	function(x, y, ...)
774	{	{
775	if (require("lattice", quietly = TRUE))	## require("lattice", quietly = TRUE) -- now via Imports
776	lapply(x, function(x) {	lapply(x, function(x) {
777	cn <- lapply(colnames(x), as.name)	cn <- lapply(colnames(x), as.name)
778	switch(min(ncol(x), 3),	switch(min(ncol(x), 3),
779	qqmath(eval(substitute(~ x,	qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
	list(x = cn[[1]]))),
	x, ...),
780	xyplot(eval(substitute(y ~ x,	xyplot(eval(substitute(y ~ x,
781	list(y = cn[[1]],	list(y = cn[[1]],
782	x = cn[[2]]))),	x = cn[[2]]))), x, ...),
	x, ...),
783	splom(~ x, ...))	splom(~ x, ...))
784	})	})
785	})	})

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.831   changed lines   Added in v.832

---

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