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

Parent Directory | Revision Log | Patch

	revision 446, Wed Jan 19 19:08:01 2005 UTC	revision 449, Fri Jan 21 22:01:21 2005 UTC
#	Line 1	Line 1
1	contr.SAS <- function(n, contrasts = TRUE)	contr.SAS <- function(n, contrasts = TRUE)
2	{	{
3	if (is.numeric(n) && length(n) == 1) contr.treatment(n, n, contrasts)	contr.treatment(n, if (is.numeric(n) && length(n) == 1) n else length(n), contrasts)
	else contr.treatment(n, length(n), contrasts)
4	}	}
5
6	lmerControl <-                            # Control parameters for lmer	lmerControl <-                            # Control parameters for lmer
#	Line 38	Line 37
37	model = TRUE, x = FALSE, y = FALSE, ...)	model = TRUE, x = FALSE, y = FALSE, ...)
38	{	{
39	# match and check parameters	# match and check parameters
40	method <- match.arg(method)	REML <- match.arg(method) == "REML"
41	controlvals <- do.call("lmerControl", control)	controlvals <- do.call("lmerControl", control)
42	controlvals$REML <- method == "REML"	controlvals$REML <- REML
43	if (length(formula) < 3) stop("formula must be a two-sided formula")	if (length(formula) < 3) stop("formula must be a two-sided formula")
44	# create the model frame as frm
45	mf <- match.call()	mf <- match.call()           # create the model frame as frm
46	m <- match(c("data", "subset", "weights", "na.action", "offset"),	m <- match(c("data", "subset", "weights", "na.action", "offset"),
47	names(mf), 0)	names(mf), 0)
48	mf <- mf[c(1, m)]	mf <- mf[c(1, m)]
#	Line 66	Line 65
65	names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))	names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
66
67	## order factor list by decreasing number of levels	## order factor list by decreasing number of levels
68	ford <- rev(order(sapply(random, function(x) length(levels(x[[2]])))))	nlev <- sapply(random, function(x) length(levels(x[[2]])))
69	if (any(ford != seq(a = random))) { # re-order both facs and random	if (any(diff(nlev) < 0)) {
70	random <- random[ford]	random <- random[rev(order(nlev))]
71	}	}
72	mmats <- c(lapply(random, "[[", 1),	mmats <- c(lapply(random, "[[", 1),
73	.fixed = list(cbind(model.matrix(nobars(formula), frm),	.fixed = list(cbind(model.matrix(nobars(formula), frm),
74	.response = model.response(frm))))	.response = model.response(frm))))
75	obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")	obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")
76	obj@call <- match.call()	obj@call <- match.call()
77	obj@REML <- method == "REML"	obj@REML <- REML
78	.Call("lmer_initial", obj, PACKAGE="Matrix")	.Call("lmer_initial", obj, PACKAGE="Matrix")
79	.Call("lmer_ECMEsteps", obj,	.Call("lmer_ECMEsteps", obj,
80	controlvals$niterEM,	controlvals$niterEM,
#	Line 83	Line 82
82	controlvals$EMverbose,	controlvals$EMverbose,
83	PACKAGE = "Matrix")	PACKAGE = "Matrix")
84	LMEoptimize(obj) <- controlvals	LMEoptimize(obj) <- controlvals
85	#fitted = .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")	#fitted <- .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")
86	#residuals = mmats$.Xy[,".response"] - fitted	#residuals <- mmats$.Xy[,".response"] - fitted
87	#if (as.logical(x)[1]) x = mmats else x = list()	#if (as.logical(x)[1]) x <- mmats else x = list()
88	#rm(mmats)	#rm(mmats)
89	obj	obj
90	})	})
#	Line 121	Line 120
120	return(x)	return(x)
121	})	})
122
	setMethod("deviance", signature(object = "lmer"),
	function(object, REML = FALSE, ...) {
	.Call("lmer_factor", object, PACKAGE = "Matrix")
	object@deviance[ifelse(REML, 2, 1)]
	})
123	setMethod("ranef", signature(object = "lmer"),	setMethod("ranef", signature(object = "lmer"),
124	function(object, ...) {	function(object, ...) {
125	.Call("lmer_ranef", object, PACKAGE = "Matrix")	.Call("lmer_ranef", object, PACKAGE = "Matrix")
#	Line 139	Line 132
132	val[-length(val)]	val[-length(val)]
133	})	})
134
	setMethod("vcov", signature(object = "lmer"),
	function(object, REML = TRUE, useScale = TRUE,...) {
	## force an "lmer_invert"
	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
	})
135	setMethod("VarCorr", signature(x = "lmer"),	setMethod("VarCorr", signature(x = "lmer"),
136	function(x, REML = TRUE, useScale = TRUE, ...) {	function(x, REML = TRUE, useScale = TRUE, ...) {
137	val = .Call("lmer_variances", x, PACKAGE = "Matrix")	val = .Call("lmer_variances", x, PACKAGE = "Matrix")
#	Line 163	Line 140
140	val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")	val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
141	}	}
142	new("VarCorr",	new("VarCorr",
143	scale = .Call("lmer_sigma", x, REML),	scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),
144	reSumry = val,	reSumry = val,
145	useScale = useScale)	useScale = useScale)
146	})	})
147
148	setMethod("gradient", signature(x = "lmer"),	setMethod("gradient", signature(x = "lmer"),
149	function(x, REML, unconst, ...)	function(x, REML, unconst, ...)
150	.Call("lmer_gradient", x, REML, unconst))	.Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix"))
151
152	setMethod("summary", "lmer",	setMethod("summary", signature(object = "lmer"),
153	function(object, REML = TRUE, useScale = TRUE, ...) {	function(object, ...)
154	fcoef <- fixef(object)	new("summary.lmer", object, useScale = TRUE, showCorrelation = TRUE))
155	corF <- as(as(vcov(object, REML, useScale), "pdmatrix"),
156	"corrmatrix")	setMethod("show", signature(object = "lmer"),
157	DF <- getFixDF(object)	function(object)
158	coefs <- cbind(fcoef, corF@stdDev, DF)	show(new("summary.lmer", object, useScale = TRUE, showCorrelation = FALSE))
159	nc <- object@nc	)
	dimnames(coefs) <-
	list(names(fcoef), c("Estimate", "Std. Error", "DF"))
	new("summary.lmer",
	coefficients = as.matrix(coefs),
	scale = .Call("lmer_sigma", object, REML),
	denomDF = as.integer(DF),
	REML = REML,
	ngrps = unlist(lapply(object@flist,
	function(x) length(levels(x)))),
	nobs = nc[length(nc)],
	corFixed = corF,
	VarCorr = VarCorr(object, REML, useScale),
	useScale = useScale,
	showCorrelation = FALSE)
	})
160
161	setMethod("show", "lmer",	setMethod("show", "summary.lmer",
162	function(object) {	function(object) {
163	fcoef <- fixef(object)	fcoef <- fixef(object)
164	corF <- as(as(vcov(object, REML = TRUE, useScale = TRUE), "pdmatrix"),	useScale <- object@useScale
165		corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"),
166	"corrmatrix")	"corrmatrix")
167	DF <- getFixDF(object)	DF <- getFixDF(object)
168	coefs <- cbind(fcoef, corF@stdDev, DF)	coefs <- cbind(fcoef, corF@stdDev, DF)
169	nc <- object@nc	nc <- object@nc
170	dimnames(coefs) <-	dimnames(coefs) <-
171	list(names(fcoef), c("Estimate", "Std. Error", "DF"))	list(names(fcoef), c("Estimate", "Std. Error", "DF"))
172	new("summary.ssclme",	digits <- max(3, getOption("digits") - 2)
173	coefficients = as.matrix(coefs),	REML <- length(object@REML) > 0 && object@REML[1]
174	scale = .Call("lmer_sigma", object, REML = TRUE),	llik <- logLik(object)
175	denomDF = as.integer(DF),	dev <- object@deviance
176	REML = TRUE,
177	ngrps = unlist(lapply(object@flist,	rdig <- 5
178	function(x) length(levels(x)))),	cat("Linear mixed-effects model fit by ")
179	nobs = nc[length(nc)],	cat(ifelse(object@REML, "REML\n", "maximum likelihood\n") )
180	corFixed = corF,	if (!is.null(object@call$formula)) {
181	VarCorr = VarCorr(object, REML = TRUE, useScale = TRUE),	cat("Formula:", deparse(object@call$formula),"\n")
182	useScale = TRUE,	}
183	showCorrelation = FALSE)	if (!is.null(object@call$data)) {
184		cat("   Data:", deparse(object@call$data), "\n")
185		}
186		if (!is.null(object@call$subset)) {
187		cat(" Subset:",
188		deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
189		}
190		print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
191		logLik = c(llik),
192		MLdeviance = dev["ML"],
193		REMLdeviance = dev["REML"],
194		row.names = ""))
195		cat("Random effects:\n")
196		show(VarCorr(object))
197		ngrps <- lapply(object@flist, function(x) length(levels(x)))
198		cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))
199		cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
200		cat("\n")
201		if (!useScale)
202		cat("\nEstimated scale (compare to 1) ",
203		.Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"),
204		"\n")
205		if (nrow(coefs) > 0) {
206		if (useScale) {
207		stat <- coefs[,1]/coefs[,2]
208		pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
209		nms <- colnames(coefs)
210		coefs <- cbind(coefs, stat, pval)
211		colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")
212		} else {
213		coefs <- coefs[, 1:2, drop = FALSE]
214		stat <- coefs[,1]/coefs[,2]
215		pval <- 2*pnorm(abs(stat), lower = FALSE)
216		nms <- colnames(coefs)
217		coefs <- cbind(coefs, stat, pval)
218		colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
219		}
220		cat("\nFixed effects:\n")
221		printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
222		if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
223		rn <- rownames(coefs)
224		dimnames(corF) <- list(
225		abbreviate(rn, minlen=11),
226		abbreviate(rn, minlen=6))
227		if (!is.null(corF)) {
228		p <- NCOL(corF)
229		if (p > 1) {
230		cat("\nCorrelation of Fixed Effects:\n")
231		corF <- format(round(corF, 3), nsmall = 3)
232		corF[!lower.tri(corF)] <- ""
233		print(corF[-1, -p, drop=FALSE], quote = FALSE)
234		}
235		}
236		}
237		}
238		invisible(object)
239	})	})
240
241	## calculates degrees of freedom for fixed effects Wald tests	## calculates degrees of freedom for fixed effects Wald tests
#	Line 234	Line 252
252	rep(n-p, p)	rep(n-p, p)
253	})	})
254
	setMethod("fitted", signature(object="lmer"),
	function(object, ...)
	{
	object@fitted
	})
	setMethod("residuals", signature(object="lmer"),
	function(object, ...) object@residuals )
255	setMethod("logLik", signature(object="lmer"),	setMethod("logLik", signature(object="lmer"),
256	function(object, REML = object@REML, ...) {	function(object, REML = object@REML, ...) {
257	val <- -deviance(object, REML = REML)/2	val <- -deviance(object, REML = REML)/2
#	Line 254	Line 263
263	val	val
264	})	})
265
	# setMethod("summary", signature(object="lme"),
	#           function(object, ...) {
	#               llik <- logLik(object)
	#               resd <- residuals(object, type="pearson")
	#               if (length(resd) > 5) {
	#                   resd <- quantile(resd)
	#                   names(resd) <- c("Min","Q1","Med","Q3","Max")
	#               }
	#               new("summary.lme",
	#                   call = object@call,
	#                   logLik = llik,
	#                   re = summary(object@rep, REML = object@REML,
	#                                useScale = TRUE),
	#                   residuals = resd)
	#           })
	# setMethod("show", signature(object = "summary.lme"),
	#           function(object)
	#       {
	#           rdig <- 5
	#           cat("Linear mixed-effects model fit by ")
	#           cat(ifelse(object@re@REML, "REML\n", "maximum likelihood\n") )
	#           if (!is.null(object@call$formula)) {
	#               cat("Fixed:", deparse(object@call$formula),"\n")
	#           }
	#           if (!is.null(object@call$data)) {
	#               cat(" Data:", deparse(object@call$data), "\n")
	#           }
	#           if (!is.null(object@call$subset)) {
	#               cat(" Subset:",
	#                   deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
	#           }
	#           llik <- object@logLik
	#           print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
	#                            logLik = c(object@logLik), row.names = ""))
	#           cat("\n")
	#           object@re@useScale <- TRUE
	#           object@re@showCorrelation <- TRUE
	#           show(object@re)
	#           invisible(object)
	#       })
266	setMethod("anova", signature(object = "lmer"),	setMethod("anova", signature(object = "lmer"),
267	function(object, ...)	function(object, ...)
268	{	{
#	Line 344	Line 311
311	foo <- object	foo <- object
312	foo@status["factored"] <- FALSE	foo@status["factored"] <- FALSE
313	.Call("lmer_factor", foo, PACKAGE="Matrix")	.Call("lmer_factor", foo, PACKAGE="Matrix")
314	dfr <- lme4:::getFixDF(foo)	dfr <- getFixDF(foo)
315	rcol <- ncol(foo@RXX)	rcol <- ncol(foo@RXX)
316	ss <- foo@RXX[ , rcol]^2	ss <- foo@RXX[ , rcol]^2
317	ssr <- ss[[rcol]]	ssr <- ss[[rcol]]
#	Line 368	Line 335
335	}	}
336	})	})
337
	setMethod("formula", "lmer", function(x, ...) x@call$formula)
	setMethod("plot", signature(x = "lmer"),
	function(x, y, ...)
	cat("plot method for lmer not yet implemented\n"))
338	setMethod("update", signature(object = "lmer"),	setMethod("update", signature(object = "lmer"),
339	function(object, formula., ..., evaluate = TRUE)	function(object, formula., ..., evaluate = TRUE)
340	{	{
#	Line 412	Line 373
373	ci <- array(NA, dim = c(length(parm), 2),	ci <- array(NA, dim = c(length(parm), 2),
374	dimnames = list(pnames[parm], pct))	dimnames = list(pnames[parm], pct))
375	ses <- sqrt(diag(vcov(object)))[parm]	ses <- sqrt(diag(vcov(object)))[parm]
376	ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object@rep)[parm], qt))	ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
377	ci	ci
378	})	})
379
380		setReplaceMethod("coef", signature(object = "lmer", value = "numeric"),
381		function(object, unconst = FALSE, ..., value)
382		.Call("lmer_coefGets", object, as.double(value),
383		unconst, PACKAGE = "Matrix"))
384
385		setMethod("coef", signature(object = "lmer"),
386		function(object, unconst = FALSE, ...) {
387		.Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
388		})
389
390		setMethod("deviance", "lmer",
391		function(object, REML = NULL, ...) {
392		.Call("lmer_factor", object, PACKAGE = "Matrix")
393		if (is.null(REML))
394		REML <- if (length(oR <- object@REML)) oR else FALSE
395		object@deviance[[ifelse(REML, "REML", "ML")]]
396		})
397
398		setMethod("chol", signature(x = "lmer"),
399		function(x, pivot = FALSE, LINPACK = pivot) {
400		x@status["factored"] <- FALSE # force a decomposition
401		.Call("lmer_factor", x, PACKAGE = "Matrix")
402		})
403
404		setMethod("solve", signature(a = "lmer", b = "missing"),
405		function(a, b, ...)
406		.Call("lmer_invert", a)
407		)
408
409		setMethod("formula", "lmer", function(x, ...) x@call$formula)
410
411		setMethod("vcov", signature(object = "lmer"),
412		function(object, REML = object@REML, useScale = TRUE,...) {
413		## force an "lmer_invert"
414		sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
415		rr <- object@RXX
416		nms <- object@cnames[[".fixed"]]
417		dimnames(rr) <- list(nms, nms)
418		nr <- nrow(rr)
419		rr <- rr[-nr, -nr, drop = FALSE]
420		rr <- rr %*% t(rr)
421		if (useScale) {
422		rr = sc^2 * rr
423		}
424		rr
425		})
426

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.446   changed lines   Added in v.449

---

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