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 435 - (view) (download)
Original Path: branches/trunk-lme4/R/lmer.R

1 : bates 435 lmerControl <- # Control parameters for lmer
2 :     function(maxIter = 50,
3 :     msMaxIter = 50,
4 :     tolerance = sqrt((.Machine$double.eps)),
5 :     niterEM = 20,
6 :     msTol = sqrt(.Machine$double.eps),
7 :     msVerbose = getOption("verbose"),
8 :     PQLmaxIt = 20,
9 :     .relStep = (.Machine$double.eps)^(1/3),
10 :     EMverbose = getOption("verbose"),
11 :     analyticGradient = TRUE,
12 :     analyticHessian=FALSE)
13 :     {
14 :     list(maxIter = maxIter,
15 :     msMaxIter = msMaxIter,
16 :     tolerance = tolerance,
17 :     niterEM = niterEM,
18 :     msTol = msTol,
19 :     msVerbose = msVerbose,
20 :     PQLmaxIt = PQLmaxIt,
21 :     .relStep = .relStep,
22 :     EMverbose=EMverbose,
23 :     analyticHessian=analyticHessian,
24 :     analyticGradient=analyticGradient)
25 :     }
26 :    
27 :     setMethod("lmer", signature(formula = "formula"),
28 :     function(formula, data,
29 :     method = c("REML", "ML"),
30 :     control = list(),
31 :     subset, weights, na.action, offset,
32 :     model = TRUE, x = FALSE, y = FALSE, ...)
33 :     {
34 :     # match and check parameters
35 :     method <- match.arg(method)
36 :     controlvals <- do.call("lmerControl", control)
37 :     controlvals$REML <- method == "REML"
38 :     if (length(formula) < 3) stop("formula must be a two-sided formula")
39 :     # create the model frame as frm
40 :     mf <- match.call()
41 :     m <- match(c("data", "subset", "weights", "na.action", "offset"),
42 :     names(mf), 0)
43 :     mf <- mf[c(1, m)]
44 :     mf[[1]] <- as.name("model.frame")
45 :     frame.form <- subbars(formula)
46 :     environment(frame.form) <- environment(formula)
47 :     mf$formula <- frame.form
48 :     mf$drop.unused.levels <- TRUE
49 :     frm <- eval(mf, parent.frame())
50 :    
51 :     ## grouping factors and model matrices for random effects
52 :     bars <- findbars(formula[[3]])
53 :     random <-
54 :     lapply(bars,
55 :     function(x) list(model.matrix(eval(substitute(~term,
56 :     list(term=x[[2]]))),
57 :     frm),
58 :     eval(substitute(as.factor(fac),
59 :     list(fac = x[[3]])), frm)))
60 :     names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
61 :    
62 :     ## order factor list by decreasing number of levels
63 :     ford <- rev(order(sapply(random, function(x) length(levels(x[[2]])))))
64 :     if (any(ford != seq(a = random))) { # re-order both facs and random
65 :     random <- random[ford]
66 :     }
67 :     mmats <- c(lapply(random, "[[", 1),
68 :     .fixed = list(cbind(model.matrix(nobars(formula), frm),
69 :     .response = model.response(frm))))
70 :     obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")
71 :     .Call("lmer_initial", obj, PACKAGE="Matrix")
72 :     .Call("lmer_ECMEsteps", obj,
73 :     controlvals$niterEM,
74 :     controlvals$REML,
75 :     controlvals$EMverbose,
76 :     PACKAGE = "Matrix")
77 :     LMEoptimize(obj) <- controlvals
78 :     obj@call <- match.call()
79 :     #fitted = .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")
80 :     #residuals = mmats$.Xy[,".response"] - fitted
81 :     #if (as.logical(x)[1]) x = mmats else x = list()
82 :     #rm(mmats)
83 :     obj
84 :     })
85 :    
86 : bates 413 setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),
87 : bates 316 function(x, value)
88 :     {
89 :     if (value$msMaxIter < 1) return(x)
90 :     st <- ccoef(x) # starting values
91 :     nc <- x@nc
92 :     nc <- nc[1:(length(nc) - 2)]
93 :     constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
94 :     fn <- function(pars) {
95 :     ccoef(x) <- pars
96 :     deviance(x, REML = value$REML)
97 :     }
98 : bates 380 gr <- if (value$analyticGradient)
99 :     function(pars) {
100 :     ccoef(x) <- pars
101 :     grad <- lme4:::gradient(x, REML = value$REML, unconst = TRUE)
102 :     grad[constr] <- -grad[constr]/pars[constr]
103 :     grad
104 :     } else NULL
105 : bates 316 optimRes <- optim(st, fn, gr,
106 :     method = "L-BFGS-B",
107 :     lower = ifelse(constr, 1e-10, -Inf),
108 : bates 362 control = list(maxit = value$msMaxIter,
109 : bates 411 trace = as.integer(value$msVerbose)))
110 : bates 316 if (optimRes$convergence != 0) {
111 :     warning(paste("optim returned message",optimRes$message,"\n"))
112 :     }
113 : bates 411 ccoef(x) <- optimRes$par
114 : bates 316 return(x)
115 :     })
116 :    
117 : bates 413 setMethod("deviance", signature(object = "lmer"),
118 : bates 316 function(object, REML = FALSE, ...) {
119 : bates 413 .Call("lmer_factor", object, PACKAGE = "Matrix")
120 : bates 316 object@deviance[ifelse(REML, 2, 1)]
121 :     })
122 :    
123 : bates 413 setMethod("ranef", signature(object = "lmer"),
124 : bates 316 function(object, ...) {
125 : bates 413 .Call("lmer_ranef", object, PACKAGE = "Matrix")
126 : bates 316 })
127 :    
128 : bates 413 setMethod("fixef", signature(object = "lmer"),
129 : bates 316 function(object, ...) {
130 : bates 413 val = .Call("lmer_fixef", object, PACKAGE = "Matrix")
131 : bates 316 names(val) = object@cnames[[".fixed"]]
132 :     val[-length(val)]
133 :     })
134 :    
135 : bates 413 setMethod("vcov", signature(object = "lmer"),
136 : bates 316 function(object, REML = TRUE, useScale = TRUE,...) {
137 : bates 413 ## force an "lmer_invert"
138 :     sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
139 : bates 316 rr <- object@RXX
140 :     nms <- object@cnames[[".fixed"]]
141 :     dimnames(rr) <- list(nms, nms)
142 :     nr <- nrow(rr)
143 :     rr <- rr[-nr, -nr, drop = FALSE]
144 :     rr <- rr %*% t(rr)
145 :     if (useScale) {
146 :     rr = sc^2 * rr
147 :     }
148 :     rr
149 :     })
150 :    
151 : bates 413 setMethod("VarCorr", signature(x = "lmer"),
152 : bates 316 function(x, REML = TRUE, useScale = TRUE, ...) {
153 : bates 413 val = .Call("lmer_variances", x, PACKAGE = "Matrix")
154 : bates 316 for (i in seq(along = val)) {
155 :     dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])
156 :     val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
157 :     }
158 :     new("VarCorr",
159 : bates 413 scale = .Call("lmer_sigma", x, REML),
160 : bates 316 reSumry = val,
161 :     useScale = useScale)
162 :     })
163 :    
164 : bates 413 setMethod("gradient", signature(x = "lmer"),
165 : bates 316 function(x, REML, unconst, ...)
166 : bates 413 .Call("lmer_gradient", x, REML, unconst))
167 : bates 316
168 : bates 413 setMethod("summary", "lmer",
169 : bates 316 function(object, REML = TRUE, useScale = TRUE, ...) {
170 :     fcoef <- fixef(object)
171 :     corF <- as(as(vcov(object, REML, useScale), "pdmatrix"),
172 :     "corrmatrix")
173 :     DF <- getFixDF(object)
174 :     coefs <- cbind(fcoef, corF@stdDev, DF)
175 :     nc <- object@nc
176 :     dimnames(coefs) <-
177 :     list(names(fcoef), c("Estimate", "Std. Error", "DF"))
178 :     new("summary.ssclme",
179 :     coefficients = as.matrix(coefs),
180 : bates 413 scale = .Call("lmer_sigma", object, REML),
181 : bates 316 denomDF = as.integer(DF),
182 :     REML = REML,
183 : bates 423 ngrps = unlist(lapply(object@flist,
184 :     function(x) length(levels(x)))),
185 : bates 316 nobs = nc[length(nc)],
186 :     corFixed = corF,
187 :     VarCorr = VarCorr(object, REML, useScale),
188 :     useScale = useScale,
189 :     showCorrelation = FALSE)
190 :     })
191 :    
192 : bates 413 setMethod("show", "lmer",
193 : bates 316 function(object) {
194 :     fcoef <- fixef(object)
195 :     corF <- as(as(vcov(object, REML = TRUE, useScale = TRUE), "pdmatrix"),
196 :     "corrmatrix")
197 :     DF <- getFixDF(object)
198 :     coefs <- cbind(fcoef, corF@stdDev, DF)
199 :     nc <- object@nc
200 :     dimnames(coefs) <-
201 :     list(names(fcoef), c("Estimate", "Std. Error", "DF"))
202 :     new("summary.ssclme",
203 :     coefficients = as.matrix(coefs),
204 : bates 413 scale = .Call("lmer_sigma", object, REML = TRUE),
205 : bates 316 denomDF = as.integer(DF),
206 :     REML = TRUE,
207 : bates 423 ngrps = unlist(lapply(object@flist,
208 :     function(x) length(levels(x)))),
209 : bates 316 nobs = nc[length(nc)],
210 :     corFixed = corF,
211 :     VarCorr = VarCorr(object, REML = TRUE, useScale = TRUE),
212 :     useScale = TRUE,
213 :     showCorrelation = FALSE)
214 :     })
215 :    
216 :     ## calculates degrees of freedom for fixed effects Wald tests
217 :     ## This is a placeholder. The answers are generally wrong. It will
218 :     ## be very tricky to decide what a 'right' answer should be with
219 :     ## crossed random effects.
220 :    
221 : bates 413 setMethod("getFixDF", signature(object="lmer"),
222 : bates 316 function(object, ...)
223 :     {
224 :     nc <- object@nc[-seq(along = object@Omega)]
225 :     p <- nc[1] - 1
226 :     n <- nc[2]
227 :     rep(n-p, p)
228 :     })
229 :    
230 : bates 413 setMethod("anova", signature(object="lmer"),
231 : bates 316 function(object, ...)
232 :     {
233 :     foo <- object
234 :     foo@status["factored"] <- FALSE
235 : bates 413 .Call("lmer_factor", foo, PACKAGE="Matrix")
236 : bates 316 dfr <- getFixDF(foo)
237 :     ss <- foo@RXX[ , ".response"]^2
238 :     ssr <- ss[[".response"]]
239 :     ss <- ss[seq(along = dfr)]
240 :     names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
241 :     # FIXME: This only gives single degree of freedom tests
242 :     ms <- ss
243 :     df <- rep(1, length(ms))
244 :     f <- ms/(ssr/dfr)
245 :     P <- pf(f, df, dfr, lower.tail = FALSE)
246 :     table <- data.frame(df, ss, ms, dfr, f, P)
247 :     dimnames(table) <-
248 :     list(names(ss),
249 :     c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
250 :     if (any(match(names(ss), "(Intercept)", nomatch = 0)))
251 :     table <- table[-1,]
252 :     attr(table, "heading") <- "Analysis of Variance Table"
253 :     class(table) <- c("anova", "data.frame")
254 :     table
255 :     })

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