SCM

SCM Repository

[matrix] Annotation of /pkg/R/lmer.R
ViewVC logotype

Annotation of /pkg/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 446 - (view) (download)
Original Path: branches/trunk-lme4/R/lmer.R

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

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge