# SCM Repository

[matrix] Diff of /pkg/R/lmer.R
 [matrix] / pkg / R / lmer.R

# Diff of /pkg/R/lmer.R

revision 774, Fri Jun 10 18:40:54 2005 UTC revision 775, Mon Jun 13 01:22:28 2005 UTC
# Line 2  Line 2
2
3  ## Some utilities  ## Some utilities
4
5    ## Return the index into the packed lower triangle
6  Lind <- function(i,j) {  Lind <- function(i,j) {
7      if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))      if (i < j) stop(paste("Index i=", i,"must be >= index j=", j))
8      ((i - 1) * i)/2 + j      ((i - 1) * i)/2 + j
9  }  }
10
11  # Return the pairs of expressions separated by vertical bars  ## Return the pairs of expressions separated by vertical bars

12  findbars <- function(term)  findbars <- function(term)
13  {  {
14      if (is.name(term) || is.numeric(term)) return(NULL)      if (is.name(term) || is.numeric(term)) return(NULL)
# Line 19  Line 19
19      c(findbars(term[[2]]), findbars(term[[3]]))      c(findbars(term[[2]]), findbars(term[[3]]))
20  }  }
21
22  # Return the formula omitting the pairs of expressions separated by vertical bars  ## Return the formula omitting the pairs of expressions
23    ## that are separated by vertical bars
24  nobars <- function(term)  nobars <- function(term)
25  {  {
26      # FIXME: is the is.name in the condition redundant?      # FIXME: is the is.name in the condition redundant?
# Line 42  Line 42
42      term      term
43  }  }
44
45  # Substitute the '+' function for the '|' function  ## Substitute the '+' function for the '|' function

46  subbars <- function(term)  subbars <- function(term)
47  {  {
48      if (is.name(term) || is.numeric(term)) return(term)      if (is.name(term) || is.numeric(term)) return(term)
# Line 58  Line 57
57      term      term
58  }  }
59
60  lmerControl <-                            # Control parameters for lmer  ## Control parameters for lmer
61    lmerControl <-
62    function(maxIter = 50,    function(maxIter = 50,
63             msMaxIter = 200,             msMaxIter = 200,
64             tolerance = sqrt((.Machine\$double.eps)),             tolerance = sqrt((.Machine\$double.eps)),
# Line 70  Line 70
71             analyticHessian=FALSE)             analyticHessian=FALSE)
72  {  {
73      list(maxIter = maxIter,      list(maxIter = as.integer(maxIter),
74           msMaxIter = msMaxIter,           msMaxIter = as.integer(msMaxIter),
75           tolerance = tolerance,           tolerance = as.double(tolerance),
76           niterEM = niterEM,           niterEM = as.integer(niterEM),
77           msTol = msTol,           msTol = as.double(msTol),
78           msVerbose = msVerbose,           msVerbose = as.logical(msVerbose),
79           PQLmaxIt = PQLmaxIt,           PQLmaxIt = as.integer(PQLmaxIt),
80           EMverbose=EMverbose,           EMverbose = as.logical(EMverbose),
83  }  }
84
85  setMethod("lmer", signature(formula = "formula"),  setMethod("lmer", signature(formula = "formula"),
# Line 185  Line 185
185
186            dev.resids <- quote(family\$dev.resids(y, mu, wtssqr))            dev.resids <- quote(family\$dev.resids(y, mu, wtssqr))
189            mu.eta <- quote(family\$mu.eta(eta))            mu.eta <- quote(family\$mu.eta(eta))
190            variance <- quote(family\$variance(mu))            variance <- quote(family\$variance(mu))
191              LMEopt <- get("LMEoptimize<-")
192              doLMEopt <- quote(LMEopt(x = mer, value = cv))
193
194            mmo <- mmats            GSpt <- .Call("glmer_init", environment())
195            mmats[[1]][1,1] <- mmats[[1]][1,1]            .Call("glmer_PQL", GSpt)  # obtain PQL estimates
conv <- FALSE
firstIter <- TRUE
msMaxIter.orig <- cv\$msMaxIter

for (iter in seq(length = cv\$PQLmaxIt))
{
dmu.deta <- eval(mu.eta) # family\$mu.eta(eta)
## weights (note: wts is already square-rooted)
.Call("glmer_weight_matrix_list", mmo,
wts * dmu.deta / sqrt(eval(variance)), ## weights
eta - offset + (y - mu) / dmu.deta, ## working residual
mmats, PACKAGE="Matrix")
.Call("lmer_update_mm", mer, mmats, PACKAGE="Matrix")
if (firstIter) {
.Call("lmer_initial", mer, PACKAGE="Matrix")
if (gVerb) cat(" PQL iterations convergence criterion\n")
}
.Call("lmer_ECMEsteps", mer, cv\$niterEM, cv\$EMverbose,
PACKAGE = "Matrix")
LMEoptimize(mer) <- cv
eta <- offset + .Call("lmer_fitted", mer, mmo, TRUE,
PACKAGE = "Matrix")
crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))
if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))
## use this to determine convergence
if (crit < cv\$tolerance) {
conv <- TRUE
break
}
etaold[] <- eta
if (firstIter) {          # Change the number of EM and optimization
cv\$niterEM <- 2       # iterations for subsequent PQL iterations.
cv\$msMaxIter <- 10
firstIter <- FALSE
}
}
if (!conv) warning("IRLS iterations for glmm did not converge")
cv\$msMaxIter <- msMaxIter.orig
196
197            fixInd <- seq(ncol(x))            fixInd <- seq(ncol(x))
198            ## pars[fixInd] == beta, pars[-fixInd] == theta            ## pars[fixInd] == beta, pars[-fixInd] == theta
199            PQLpars <- c(fixef(mer),            PQLpars <- c(fixef(mer),
200                         .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))                         .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))
201            env <- environment()            ## set flag to skip fixed-effects in subsequent calls
202              mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
203
204              nAGQ <- 0
205              if (method == "Laplace") nAGQ <- 1
206
207            devLaplace <- function(pars)            devLaplace <- function(pars)
208                .Call("lmer_devLaplace", pars, cv\$tolerance, env, PACKAGE = "Matrix")                .Call("glmer_devLaplace", pars, GSpt, PACKAGE = "Matrix")
209
210            if (method == "Laplace") {            if (method == "Laplace") {
211                nc <- mer@nc                nc <- mer@nc
212                const <- c(rep(FALSE, length(fixInd)),                const <- c(rep(FALSE, length(fixInd)),
213                           unlist(lapply(nc[1:(length(nc) - 2)],                           unlist(lapply(nc[1:(length(nc) - 2)],
214                                         function(k) 1:((k*(k+1))/2) <= k)))                                         function(k) 1:((k*(k+1))/2) <= k)))
215                ## set flag to skip fixed-effects in subsequent mer computations
mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
216                RV <- lapply(R.Version()[c("major", "minor")], as.numeric)                RV <- lapply(R.Version()[c("major", "minor")], as.numeric)
217                if (RV\$major == 2 && RV\$minor >= 2.0) {                if (RV\$major == 2 && RV\$minor >= 2.0) {
218                    optimRes <-                    optimRes <-
# Line 277  Line 244
244                loglik <- -optimRes\$objective/2                loglik <- -optimRes\$objective/2
245                fxd <- optpars[fixInd]                fxd <- optpars[fixInd]
246                names(fxd) <- names(PQLpars)[fixInd]                names(fxd) <- names(PQLpars)[fixInd]
## reset flag to skip fixed-effects in mer computations
mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
247            } else {            } else {
248                loglik <- -devLaplace(PQLpars)/2                loglik <- -devLaplace(PQLpars)/2
249                fxd <- PQLpars[fixInd]                fxd <- PQLpars[fixInd]
250            }            }
251
252              ## reset flag to skip fixed-effects in subsequent calls
253              mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
254
255            attributes(loglik) <- attributes(logLik(mer))            attributes(loglik) <- attributes(logLik(mer))
256            new("lmer", mer, frame = frm, terms = glm.fit\$terms,            new("lmer", mer, frame = frm, terms = glm.fit\$terms,
257                assign = attr(glm.fit\$x, "assign"), call = match.call(),                assign = attr(glm.fit\$x, "assign"), call = match.call(),

Legend:
 Removed from v.774 changed lines Added in v.775