| 1 |
|
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 |
setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"), |
setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"), |
| 87 |
function(x, value) |
function(x, value) |
| 88 |
{ |
{ |
Parent Directory
| 
Revision Log
| 
Patch
