SCM

SCM Repository

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

Diff of /pkg/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 773, Thu Jun 9 18:30:33 2005 UTC revision 774, Fri Jun 10 18:40:54 2005 UTC
# Line 179  Line 179 
179            gVerb <- getOption("verbose")            gVerb <- getOption("verbose")
180            etaold <- eta <- glm.fit$linear.predictors            etaold <- eta <- glm.fit$linear.predictors
181            wts <- glm.fit$prior.weights            wts <- glm.fit$prior.weights
182              wtssqr <- wts * wts
183            offset <- glm.fit$offset            offset <- glm.fit$offset
184            if (is.null(offset)) offset <- numeric(length(eta))            if (is.null(offset)) offset <- numeric(length(eta))
185    
186            dev.resids <- quote(family$dev.resids(y, mu, wts*wts))            dev.resids <- quote(family$dev.resids(y, mu, wtssqr))
187            linkinv <- quote(family$linkinv(eta))            linkinv <- quote(family$linkinv(eta))
188            mu.eta <- quote(family$mu.eta(eta))            mu.eta <- quote(family$mu.eta(eta))
189            variance <- quote(family$variance(mu))            variance <- quote(family$variance(mu))
# Line 229  Line 230 
230            if (!conv) warning("IRLS iterations for glmm did not converge")            if (!conv) warning("IRLS iterations for glmm did not converge")
231            cv$msMaxIter <- msMaxIter.orig            cv$msMaxIter <- msMaxIter.orig
232    
233            rdobj <- .Call("lmer_collapse", mer, PACKAGE = "Matrix")            fixInd <- seq(ncol(x))
234            unwtd <- mmo            ## pars[fixInd] == beta, pars[-fixInd] == theta
235            unwtd$.fixed <- as.matrix(y)            PQLpars <- c(fixef(mer),
236            wtd <- unwtd                         .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))
           wtd[[1]][1,1] <- wtd[[1]][1,1]  
           fixInd <- seq(along = fixef(mer))  
237            env <- environment()            env <- environment()
238            bhat <- function(pars) { #pars[fixInd] == beta, pars[-fixInd] == theta  
239  #              if (cv$analyticHessian) {            devLaplace <- function(pars)
240                    .Call("glmer_bhat_iterate", pars, cv$tolerance,                .Call("lmer_devLaplace", pars, cv$tolerance, env, PACKAGE = "Matrix")
                         env, PACKAGE = "Matrix")  
 ##              } else {  
 ##                  etaold <- off <- eta <- drop(glm.fit$x %*% pars[fixInd]) + offset  
 ##                  .Call("lmer_coefGets", rdobj, as.double(pars[-fixInd]),  
 ##                        2, PACKAGE = "Matrix")  
 ##                  niter <- 20  
 ##                  conv <- FALSE  
 ##                  for (iter in seq(length = niter)) {  
 ##                      mu <- eval(linkinv)  
 ##                      dmu.deta <- eval(mu.eta)  
 ##                      .Call("glmer_weight_matrix_list", unwtd,  
 ##                            wts * dmu.deta / sqrt(eval(variance)),  
 ##                            eta - off + (y - mu)/dmu.deta, wtd,  
 ##                            PACKAGE="Matrix")  
 ##                      .Call("lmer_update_mm", rdobj, wtd,  
 ##                            PACKAGE="Matrix")  
 ##                      eta <- off + .Call("lmer_fitted", rdobj,  
 ##                                         unwtd, TRUE, PACKAGE = "Matrix")  
 ##                      if (max(abs(eta - etaold)) <  
 ##                          (0.1 + max(abs(eta))) * cv$tolerance) {  
 ##                          conv <- TRUE  
 ##                          break  
 ##                      }  
 ##                      etaold <- eta  
 ##                  }  
 ##                  if (!conv) warning("iterations for bhat did not converge")  
 ##              }  
 ##              invisible(eval(linkinv))  
           }  
   
           devLaplace <- function(pars) {  
               ## FIXME: This actually returns half the deviance.  
               bhat(pars)  
               dev <- sum(family$dev.resids(y, mu, wts*wts))/2  
               lap <- .Call("glmer_Laplace_devComp", rdobj, PACKAGE = "Matrix")  
               if (gVerb) cat(sprintf("  %#11g %#11g\n", dev, lap))  
               dev - lap  
           }  
241    
242            if (method == "Laplace") {            if (method == "Laplace") {
243                nc <- mer@nc                nc <- mer@nc
244                const <- c(rep(FALSE, length(fixef(mer))),                const <- c(rep(FALSE, length(fixInd)),
245                           unlist(lapply(nc[1:(length(nc) - 2)],                           unlist(lapply(nc[1:(length(nc) - 2)],
246                                         function(k) 1:((k*(k+1))/2) <= k)))                                         function(k) 1:((k*(k+1))/2) <= k)))
247                  ## set flag to skip fixed-effects in subsequent mer computations
248                  mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
249                RV <- lapply(R.Version()[c("major", "minor")], as.numeric)                RV <- lapply(R.Version()[c("major", "minor")], as.numeric)
250                if (RV$major == 2 && RV$minor >= 2.0) {                if (RV$major == 2 && RV$minor >= 2.0) {
251                    optimRes <-                    optimRes <-
252                        nlminb(c(fixef(mer),                        nlminb(PQLpars, devLaplace,
                                .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")),  
                              devLaplace,  
253                               lower = ifelse(const, 5e-10, -Inf),                               lower = ifelse(const, 5e-10, -Inf),
254                               control = list(trace = getOption("verbose"),                               control = list(trace = getOption("verbose"),
255                               iter.max = cv$msMaxIter))                               iter.max = cv$msMaxIter))
# Line 297  Line 258 
258                        warning("nlminb failed to converge")                        warning("nlminb failed to converge")
259                } else {                } else {
260                    optimRes <-                    optimRes <-
261                        optim(c(fixef(mer),                        optim(PQLpars, devLaplace, method = "L-BFGS-B",
                               .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")),  
                             devLaplace, method = "L-BFGS-B",  
262                              lower = ifelse(const, 5e-10, -Inf),                              lower = ifelse(const, 5e-10, -Inf),
263                              control = list(trace = getOption("verbose"),                              control = list(trace = getOption("verbose"),
264                                   reltol = cv$msTol, maxit = cv$msMaxIter))                                   reltol = cv$msTol, maxit = cv$msMaxIter))
# Line 308  Line 267 
267                        warning("optim failed to converge")                        warning("optim failed to converge")
268                }                }
269    
270                if (getOption("verbose")) {                if (gVerb) {
271                    cat(paste("convergence message", optimRes$message, "\n"))                    cat(paste("convergence message", optimRes$message, "\n"))
272                    cat("Fixed effects:\n")                    cat("Fixed effects:\n")
273                    print(optimRes$par[fixInd])                    print(optimRes$par[fixInd])
274                    cat("(box constrained) variance coefficients:\n")                    cat("(box constrained) variance coefficients:\n")
275                    print(optimRes$par[-fixInd])                    print(optimRes$par[-fixInd])
276                }                }
277                  loglik <- -optimRes$objective/2
278                  fxd <- optpars[fixInd]
279                  names(fxd) <- names(PQLpars)[fixInd]
280                  ## reset flag to skip fixed-effects in mer computations
281                  mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
282            } else {            } else {
283                optpars <- c(fixef(mer),                loglik <- -devLaplace(PQLpars)/2
284                             .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))                fxd <- PQLpars[fixInd]
285            }            }
286    
           loglik <- -devLaplace(optpars)  
287            attributes(loglik) <- attributes(logLik(mer))            attributes(loglik) <- attributes(logLik(mer))
           ff <- optpars[fixInd]  
           names(ff) <- names(fixef(mer))  
           .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")  
288            new("lmer", mer, frame = frm, terms = glm.fit$terms,            new("lmer", mer, frame = frm, terms = glm.fit$terms,
289                assign = attr(glm.fit$x, "assign"), call = match.call(),                assign = attr(glm.fit$x, "assign"), call = match.call(),
290                family = family, logLik = loglik, fixed = ff)                family = family, logLik = loglik, fixed = fxd)
291        })        })
292    
293  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
# Line 387  Line 347 
347            })            })
348    
349  setMethod("fixef", signature(object = "mer"),  setMethod("fixef", signature(object = "mer"),
350            function(object, ...) {            function(object, ...)
351                val <- .Call("lmer_fixef", object, PACKAGE = "Matrix")                .Call("lmer_fixef", object, PACKAGE = "Matrix"))
352                val[-length(val)]  
           })  
353    
354  setMethod("fixef", signature(object = "lmer"),  setMethod("fixef", signature(object = "lmer"),
355            function(object, ...) object@fixed)            function(object, ...) object@fixed)

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

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