Annotation of /branches/trunk-lme4/R/lmer.R

Parent Directory | Revision Log

---

Revision 457 - (view) (download)

1 :	bates	446	contr.SAS <- function(n, contrasts = TRUE)
2 :			{
3 :	bates	449	contr.treatment(n, if (is.numeric(n) && length(n) == 1) n else length(n), contrasts)
4 :	bates	446	}
5 :
6 :	bates	435	lmerControl <- # Control parameters for lmer
7 :			function(maxIter = 50,
8 :			msMaxIter = 50,
9 :			tolerance = sqrt((.Machine$double.eps)),
10 :			niterEM = 20,
11 :			msTol = sqrt(.Machine$double.eps),
12 :			msVerbose = getOption("verbose"),
13 :			PQLmaxIt = 20,
14 :			.relStep = (.Machine$double.eps)^(1/3),
15 :			EMverbose = getOption("verbose"),
16 :			analyticGradient = TRUE,
17 :			analyticHessian=FALSE)
18 :			{
19 :			list(maxIter = maxIter,
20 :			msMaxIter = msMaxIter,
21 :			tolerance = tolerance,
22 :			niterEM = niterEM,
23 :			msTol = msTol,
24 :			msVerbose = msVerbose,
25 :			PQLmaxIt = PQLmaxIt,
26 :			.relStep = .relStep,
27 :			EMverbose=EMverbose,
28 :			analyticHessian=analyticHessian,
29 :			analyticGradient=analyticGradient)
30 :			}
31 :
32 :			setMethod("lmer", signature(formula = "formula"),
33 :			function(formula, data,
34 :			method = c("REML", "ML"),
35 :			control = list(),
36 :			subset, weights, na.action, offset,
37 :			model = TRUE, x = FALSE, y = FALSE, ...)
38 :			{
39 :			# match and check parameters
40 :	bates	449	REML <- match.arg(method) == "REML"
41 :	bates	435	controlvals <- do.call("lmerControl", control)
42 :	bates	449	controlvals$REML <- REML
43 :	bates	435	if (length(formula) < 3) stop("formula must be a two-sided formula")
44 :	bates	449
45 :			mf <- match.call() # create the model frame as frm
46 :	bates	435	m <- match(c("data", "subset", "weights", "na.action", "offset"),
47 :			names(mf), 0)
48 :			mf <- mf[c(1, m)]
49 :			mf[[1]] <- as.name("model.frame")
50 :			frame.form <- subbars(formula)
51 :			environment(frame.form) <- environment(formula)
52 :			mf$formula <- frame.form
53 :			mf$drop.unused.levels <- TRUE
54 :			frm <- eval(mf, parent.frame())
55 :	bates	449
56 :	bates	435	## grouping factors and model matrices for random effects
57 :			bars <- findbars(formula[[3]])
58 :			random <-
59 :			lapply(bars,
60 :			function(x) list(model.matrix(eval(substitute(~term,
61 :			list(term=x[[2]]))),
62 :			frm),
63 :	bates	452	eval(substitute(as.factor(fac)[,drop = TRUE],
64 :	bates	435	list(fac = x[[3]])), frm)))
65 :			names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
66 :	bates	449
67 :	bates	435	## order factor list by decreasing number of levels
68 :	bates	449	nlev <- sapply(random, function(x) length(levels(x[[2]])))
69 :	bates	452	if (any(diff(nlev) > 0)) {
70 :	bates	449	random <- random[rev(order(nlev))]
71 :	bates	435	}
72 :	bates	452	fixed.form <- nobars(formula)
73 :			if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula
74 :	bates	435	mmats <- c(lapply(random, "[[", 1),
75 :	bates	452	.fixed = list(cbind(model.matrix(fixed.form, frm),
76 :	bates	435	.response = model.response(frm))))
77 :			obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")
78 :	bates	446	obj@call <- match.call()
79 :	bates	449	obj@REML <- REML
80 :	bates	435	.Call("lmer_initial", obj, PACKAGE="Matrix")
81 :			.Call("lmer_ECMEsteps", obj,
82 :			controlvals$niterEM,
83 :			controlvals$REML,
84 :			controlvals$EMverbose,
85 :			PACKAGE = "Matrix")
86 :			LMEoptimize(obj) <- controlvals
87 :	bates	449	#fitted <- .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")
88 :			#residuals <- mmats$.Xy[,".response"] - fitted
89 :			#if (as.logical(x)[1]) x <- mmats else x = list()
90 :	bates	435	#rm(mmats)
91 :			obj
92 :			})
93 :
94 :	bates	413	setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),
95 :	bates	316	function(x, value)
96 :			{
97 :			if (value$msMaxIter < 1) return(x)
98 :			st <- ccoef(x) # starting values
99 :			nc <- x@nc
100 :			nc <- nc[1:(length(nc) - 2)]
101 :			constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
102 :			fn <- function(pars) {
103 :			ccoef(x) <- pars
104 :			deviance(x, REML = value$REML)
105 :			}
106 :	bates	380	gr <- if (value$analyticGradient)
107 :			function(pars) {
108 :	bates	457	if (!all.equal(pars, ccoef(x))) ccoef(x) <- pars
109 :			grad <- gradient(x, REML = value$REML, unconst = TRUE)
110 :	bates	380	grad[constr] <- -grad[constr]/pars[constr]
111 :			grad
112 :			} else NULL
113 :	bates	316	optimRes <- optim(st, fn, gr,
114 :			method = "L-BFGS-B",
115 :			lower = ifelse(constr, 1e-10, -Inf),
116 :	bates	362	control = list(maxit = value$msMaxIter,
117 :	bates	411	trace = as.integer(value$msVerbose)))
118 :	bates	316	if (optimRes$convergence != 0) {
119 :			warning(paste("optim returned message",optimRes$message,"\n"))
120 :			}
121 :	bates	411	ccoef(x) <- optimRes$par
122 :	bates	316	return(x)
123 :			})
124 :
125 :	bates	413	setMethod("ranef", signature(object = "lmer"),
126 :	bates	316	function(object, ...) {
127 :	bates	413	.Call("lmer_ranef", object, PACKAGE = "Matrix")
128 :	bates	316	})
129 :
130 :	bates	413	setMethod("fixef", signature(object = "lmer"),
131 :	bates	316	function(object, ...) {
132 :	bates	413	val = .Call("lmer_fixef", object, PACKAGE = "Matrix")
133 :	bates	316	names(val) = object@cnames[[".fixed"]]
134 :			val[-length(val)]
135 :			})
136 :
137 :	bates	413	setMethod("VarCorr", signature(x = "lmer"),
138 :	bates	316	function(x, REML = TRUE, useScale = TRUE, ...) {
139 :	bates	413	val = .Call("lmer_variances", x, PACKAGE = "Matrix")
140 :	bates	316	for (i in seq(along = val)) {
141 :			dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])
142 :			val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
143 :			}
144 :			new("VarCorr",
145 :	bates	449	scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),
146 :	bates	316	reSumry = val,
147 :			useScale = useScale)
148 :			})
149 :
150 :	bates	413	setMethod("gradient", signature(x = "lmer"),
151 :	bates	316	function(x, REML, unconst, ...)
152 :	bates	449	.Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix"))
153 :	bates	316
154 :	bates	449	setMethod("summary", signature(object = "lmer"),
155 :			function(object, ...)
156 :			new("summary.lmer", object, useScale = TRUE, showCorrelation = TRUE))
157 :	bates	316
158 :	bates	449	setMethod("show", signature(object = "lmer"),
159 :			function(object)
160 :			show(new("summary.lmer", object, useScale = TRUE, showCorrelation = FALSE))
161 :			)
162 :
163 :			setMethod("show", "summary.lmer",
164 :	bates	316	function(object) {
165 :			fcoef <- fixef(object)
166 :	bates	449	useScale <- object@useScale
167 :			corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"),
168 :	bates	316	"corrmatrix")
169 :			DF <- getFixDF(object)
170 :			coefs <- cbind(fcoef, corF@stdDev, DF)
171 :			nc <- object@nc
172 :			dimnames(coefs) <-
173 :			list(names(fcoef), c("Estimate", "Std. Error", "DF"))
174 :	bates	449	digits <- max(3, getOption("digits") - 2)
175 :			REML <- length(object@REML) > 0 && object@REML[1]
176 :			llik <- logLik(object)
177 :			dev <- object@deviance
178 :
179 :			rdig <- 5
180 :			cat("Linear mixed-effects model fit by ")
181 :			cat(ifelse(object@REML, "REML\n", "maximum likelihood\n") )
182 :			if (!is.null(object@call$formula)) {
183 :			cat("Formula:", deparse(object@call$formula),"\n")
184 :			}
185 :			if (!is.null(object@call$data)) {
186 :			cat(" Data:", deparse(object@call$data), "\n")
187 :			}
188 :			if (!is.null(object@call$subset)) {
189 :			cat(" Subset:",
190 :			deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
191 :			}
192 :			print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
193 :			logLik = c(llik),
194 :			MLdeviance = dev["ML"],
195 :			REMLdeviance = dev["REML"],
196 :			row.names = ""))
197 :			cat("Random effects:\n")
198 :			show(VarCorr(object))
199 :			ngrps <- lapply(object@flist, function(x) length(levels(x)))
200 :			cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))
201 :			cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
202 :			cat("\n")
203 :			if (!useScale)
204 :			cat("\nEstimated scale (compare to 1) ",
205 :			.Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"),
206 :			"\n")
207 :			if (nrow(coefs) > 0) {
208 :			if (useScale) {
209 :			stat <- coefs[,1]/coefs[,2]
210 :			pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
211 :			nms <- colnames(coefs)
212 :			coefs <- cbind(coefs, stat, pval)
213 :			colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")
214 :			} else {
215 :			coefs <- coefs[, 1:2, drop = FALSE]
216 :			stat <- coefs[,1]/coefs[,2]
217 :			pval <- 2*pnorm(abs(stat), lower = FALSE)
218 :			nms <- colnames(coefs)
219 :			coefs <- cbind(coefs, stat, pval)
220 :			colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
221 :			}
222 :			cat("\nFixed effects:\n")
223 :			printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
224 :			if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
225 :			rn <- rownames(coefs)
226 :			dimnames(corF) <- list(
227 :			abbreviate(rn, minlen=11),
228 :			abbreviate(rn, minlen=6))
229 :			if (!is.null(corF)) {
230 :			p <- NCOL(corF)
231 :			if (p > 1) {
232 :			cat("\nCorrelation of Fixed Effects:\n")
233 :			corF <- format(round(corF, 3), nsmall = 3)
234 :			corF[!lower.tri(corF)] <- ""
235 :			print(corF[-1, -p, drop=FALSE], quote = FALSE)
236 :			}
237 :			}
238 :			}
239 :			}
240 :			invisible(object)
241 :	bates	316	})
242 :
243 :			## calculates degrees of freedom for fixed effects Wald tests
244 :			## This is a placeholder. The answers are generally wrong. It will
245 :			## be very tricky to decide what a 'right' answer should be with
246 :			## crossed random effects.
247 :
248 :	bates	413	setMethod("getFixDF", signature(object="lmer"),
249 :	bates	316	function(object, ...)
250 :			{
251 :			nc <- object@nc[-seq(along = object@Omega)]
252 :			p <- nc[1] - 1
253 :			n <- nc[2]
254 :			rep(n-p, p)
255 :			})
256 :
257 :	bates	446	setMethod("logLik", signature(object="lmer"),
258 :			function(object, REML = object@REML, ...) {
259 :			val <- -deviance(object, REML = REML)/2
260 :			nc <- object@nc[-seq(a = object@Omega)]
261 :			attr(val, "nall") <- attr(val, "nobs") <- nc[2]
262 :			attr(val, "df") <- nc[1] + length(coef(object))
263 :			attr(val, "REML") <- REML
264 :			class(val) <- "logLik"
265 :			val
266 :			})
267 :
268 :			setMethod("anova", signature(object = "lmer"),
269 :			function(object, ...)
270 :			{
271 :			mCall <- match.call(expand.dots = TRUE)
272 :			dots <- list(...)
273 :			modp <- logical(0)
274 :			if (length(dots))
275 :			modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm")
276 :			if (any(modp)) { # multiple models - form table
277 :			opts <- dots[!modp]
278 :			mods <- c(list(object), dots[modp])
279 :			names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)
280 :			mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE), attr, "df"))]
281 :			calls <- lapply(mods, slot, "call")
282 :			data <- lapply(calls, "[[", "data")
283 :			if (any(data != data[[1]])) stop("all models must be fit to the same data object")
284 :			header <- paste("Data:", data[[1]])
285 :			subset <- lapply(calls, "[[", "subset")
286 :			if (any(subset != subset[[1]])) stop("all models must use the same subset")
287 :			if (!is.null(subset[[1]]))
288 :			header <-
289 :			c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
290 :			llks <- lapply(mods, logLik, REML = FALSE)
291 :			Df <- sapply(llks, attr, "df")
292 :			llk <- unlist(llks)
293 :			chisq <- 2 * pmax(0, c(NA, diff(llk)))
294 :			dfChisq <- c(NA, diff(Df))
295 :			val <- data.frame(Df = Df,
296 :			AIC = sapply(llks, AIC),
297 :			BIC = sapply(llks, BIC),
298 :			logLik = llk,
299 :			"Chisq" = chisq,
300 :			"Chi Df" = dfChisq,
301 :			"Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
302 :			check.names = FALSE)
303 :			class(val) <- c("anova", class(val))
304 :			attr(val, "heading") <-
305 :			c(header, "", "Models:",
306 :			paste(names(mods),
307 :			unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
308 :			sep = ": "),"")
309 :			return(val)
310 :			} else {
311 :			# mCall$terms <- object@terms
312 :			# mCall$assign <- object@assign
313 :	bates	316	foo <- object
314 :			foo@status["factored"] <- FALSE
315 :	bates	413	.Call("lmer_factor", foo, PACKAGE="Matrix")
316 :	bates	449	dfr <- getFixDF(foo)
317 :	bates	446	rcol <- ncol(foo@RXX)
318 :			ss <- foo@RXX[ , rcol]^2
319 :			ssr <- ss[[rcol]]
320 :	bates	316	ss <- ss[seq(along = dfr)]
321 :			names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
322 :			# FIXME: This only gives single degree of freedom tests
323 :			ms <- ss
324 :			df <- rep(1, length(ms))
325 :			f <- ms/(ssr/dfr)
326 :			P <- pf(f, df, dfr, lower.tail = FALSE)
327 :			table <- data.frame(df, ss, ms, dfr, f, P)
328 :			dimnames(table) <-
329 :			list(names(ss),
330 :			c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
331 :			if (any(match(names(ss), "(Intercept)", nomatch = 0)))
332 :			table <- table[-1,]
333 :			attr(table, "heading") <- "Analysis of Variance Table"
334 :			class(table) <- c("anova", "data.frame")
335 :			table
336 :	bates	446
337 :			}
338 :	bates	316	})
339 :	bates	446
340 :			setMethod("update", signature(object = "lmer"),
341 :			function(object, formula., ..., evaluate = TRUE)
342 :			{
343 :			call <- object@call
344 :			if (is.null(call))
345 :			stop("need an object with call component")
346 :			extras <- match.call(expand.dots = FALSE)$...
347 :			if (!missing(formula.))
348 :			call$formula <- update.formula(formula(object), formula.)
349 :			if (length(extras) > 0) {
350 :			existing <- !is.na(match(names(extras), names(call)))
351 :			for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
352 :			if (any(!existing)) {
353 :			call <- c(as.list(call), extras[!existing])
354 :			call <- as.call(call)
355 :			}
356 :			}
357 :			if (evaluate)
358 :			eval(call, parent.frame())
359 :			else call
360 :			})
361 :
362 :
363 :			setMethod("confint", signature(object = "lmer"),
364 :			function (object, parm, level = 0.95, ...)
365 :			{
366 :			cf <- fixef(object)
367 :			pnames <- names(cf)
368 :			if (missing(parm))
369 :			parm <- seq(along = pnames)
370 :			else if (is.character(parm))
371 :			parm <- match(parm, pnames, nomatch = 0)
372 :			a <- (1 - level)/2
373 :			a <- c(a, 1 - a)
374 :			pct <- paste(round(100 * a, 1), "%")
375 :			ci <- array(NA, dim = c(length(parm), 2),
376 :			dimnames = list(pnames[parm], pct))
377 :			ses <- sqrt(diag(vcov(object)))[parm]
378 :	bates	449	ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
379 :	bates	446	ci
380 :			})
381 :
382 :	bates	449	setReplaceMethod("coef", signature(object = "lmer", value = "numeric"),
383 :			function(object, unconst = FALSE, ..., value)
384 :			.Call("lmer_coefGets", object, as.double(value),
385 :			unconst, PACKAGE = "Matrix"))
386 :	bates	446
387 :	bates	449	setMethod("coef", signature(object = "lmer"),
388 :			function(object, unconst = FALSE, ...) {
389 :			.Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
390 :			})
391 :	bates	446
392 :	bates	449	setMethod("deviance", "lmer",
393 :			function(object, REML = NULL, ...) {
394 :			.Call("lmer_factor", object, PACKAGE = "Matrix")
395 :			if (is.null(REML))
396 :			REML <- if (length(oR <- object@REML)) oR else FALSE
397 :			object@deviance[[ifelse(REML, "REML", "ML")]]
398 :			})
399 :	bates	446
400 :	bates	449	setMethod("chol", signature(x = "lmer"),
401 :			function(x, pivot = FALSE, LINPACK = pivot) {
402 :			x@status["factored"] <- FALSE # force a decomposition
403 :			.Call("lmer_factor", x, PACKAGE = "Matrix")
404 :			})
405 :
406 :			setMethod("solve", signature(a = "lmer", b = "missing"),
407 :			function(a, b, ...)
408 :			.Call("lmer_invert", a)
409 :			)
410 :
411 :			setMethod("formula", "lmer", function(x, ...) x@call$formula)
412 :
413 :			setMethod("vcov", signature(object = "lmer"),
414 :			function(object, REML = object@REML, useScale = TRUE,...) {
415 :			## force an "lmer_invert"
416 :			sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
417 :			rr <- object@RXX
418 :			nms <- object@cnames[[".fixed"]]
419 :			dimnames(rr) <- list(nms, nms)
420 :			nr <- nrow(rr)
421 :			rr <- rr[-nr, -nr, drop = FALSE]
422 :			rr <- rr %*% t(rr)
423 :			if (useScale) {
424 :			rr = sc^2 * rr
425 :			}
426 :			rr
427 :			})
428 :

---

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