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 776, Mon Jun 13 19:45:04 2005 UTC revision 777, Tue Jun 14 04:02:57 2005 UTC
# Line 201  Line 201 
201                         .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))                         .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))
202            ## set flag to skip fixed-effects in subsequent calls            ## set flag to skip fixed-effects in subsequent calls
203            mer@nc[length(mmats)] <- -mer@nc[length(mmats)]            mer@nc[length(mmats)] <- -mer@nc[length(mmats)]
204              ## indicator of constrained parameters
205            nAGQ <- 0            const <- c(rep(FALSE, length(fixInd)),
206                         unlist(lapply(mer@nc[seq(along = random)],
207            devLaplace <- function(pars)                                     function(k) 1:((k*(k+1))/2) <= k)
208                .Call("glmer_devAGQ", pars, GSpt, 1, PACKAGE = "Matrix")                              ))
   
209            devAGQ <- function(pars)            devAGQ <- function(pars)
210                .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")                .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
211    
212            if (method == "Laplace") {            nAGQ <- 1            # Laplacian approximation
213                nAGQ <- 1            loglik <- -devAGQ(PQLpars)/2
214                nc <- mer@nc            fxd <- PQLpars[fixInd]
               const <- c(rep(FALSE, length(fixInd)),  
                          unlist(lapply(nc[1:(length(nc) - 2)],  
                                        function(k) 1:((k*(k+1))/2) <= k)))  
215    
216                RV <- lapply(R.Version()[c("major", "minor")], as.numeric)            if (method %in% c("Laplace", "AGQ")) {
217                if (RV$major == 2 && RV$minor >= 2.0) {                if (exists("nlminb", mode = "function")) {
218                    optimRes <-                    optimRes <-
219                        nlminb(PQLpars, devLaplace,                        nlminb(PQLpars, devAGQ,
220                               lower = ifelse(const, 5e-10, -Inf),                               lower = ifelse(const, 5e-10, -Inf),
221                               control = list(trace = getOption("verbose"),                               control = list(trace = getOption("verbose"),
222                               iter.max = cv$msMaxIter))                               iter.max = cv$msMaxIter))
# Line 230  Line 226 
226                    loglik <- -optimRes$objective/2                    loglik <- -optimRes$objective/2
227                } else {                } else {
228                    optimRes <-                    optimRes <-
229                        optim(PQLpars, devLaplace, method = "L-BFGS-B",                        optim(PQLpars, devAGQ, method = "L-BFGS-B",
230                              lower = ifelse(const, 5e-10, -Inf),                              lower = ifelse(const, 5e-10, -Inf),
231                              control = list(trace = getOption("verbose"),                              control = list(trace = getOption("verbose"),
232                                             maxit = cv$msMaxIter))                                             maxit = cv$msMaxIter))
# Line 239  Line 235 
235                        warning("optim failed to converge")                        warning("optim failed to converge")
236                    loglik <- -optimRes$value                    loglik <- -optimRes$value
237                }                }
   
238                if (gVerb) {                if (gVerb) {
239                    cat(paste("convergence message", optimRes$message, "\n"))                    cat(paste("convergence message", optimRes$message, "\n"))
                   cat("Fixed effects:\n")  
                   print(optimRes$par[fixInd])  
                   cat("(box constrained) variance coefficients:\n")  
                   print(optimRes$par[-fixInd])  
240                }                }
241                fxd <- optpars[fixInd]                fxd[] <- optpars[fixInd]  ## preserve the names
               names(fxd) <- names(PQLpars)[fixInd]  
           } else {  
               loglik <- -devLaplace(PQLpars)/2  
               fxd <- PQLpars[fixInd]  
242            }            }
243    
244            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
245    
           ## reset flag to skip fixed-effects in subsequent calls  
           mer@nc[length(mmats)] <- -mer@nc[length(mmats)]  
   
246            attributes(loglik) <- attributes(logLik(mer))            attributes(loglik) <- attributes(logLik(mer))
247            new("lmer", mer, frame = frm, terms = glm.fit$terms,            new("lmer", mer, frame = frm, terms = glm.fit$terms,
248                assign = attr(glm.fit$x, "assign"), call = match.call(),                assign = attr(glm.fit$x, "assign"),
249                family = family, logLik = loglik, fixed = fxd)                call = match.call(), family = family,
250                  logLik = loglik, fixed = fxd)
251        })        })
252    
253  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
# Line 284  Line 269 
269                                   .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")                                   .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
270                               .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")                               .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
271                           }                           }
272                   RV <- lapply(R.Version()[c("major", "minor")], as.numeric)                   optimRes <-
273                   if (RV$major == 2 && RV$minor >= 2.0) {                       if (exists("nlminb", mode = "function"))
274                       optimRes <- nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),                           nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
275                                          fn, gr,                                          fn, gr,
276                                          lower = ifelse(constr, 5e-10, -Inf),                                          lower = ifelse(constr, 5e-10, -Inf),
277                                          control = list(iter.max = value$msMaxIter,                                          control = list(iter.max = value$msMaxIter,
278                                          trace = as.integer(value$msVerbose)))                                          trace = as.integer(value$msVerbose)))
279                   } else {                       else
280                       optimRes <- optim(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),                           optim(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
281                                         fn, gr, method = "L-BFGS-B",                                         fn, gr, method = "L-BFGS-B",
282                                         lower = ifelse(constr, 5e-10, -Inf),                                         lower = ifelse(constr, 5e-10, -Inf),
283                                         control = list(maxit = value$msMaxIter,                                         control = list(maxit = value$msMaxIter,
284                                         trace = as.integer(value$msVerbose)))                                         trace = as.integer(value$msVerbose)))
                  }  
285                   .Call("lmer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")                   .Call("lmer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
286                   if (optimRes$convergence != 0) {                   if (optimRes$convergence != 0) {
287                       warning(paste("optim returned message",optimRes$message,"\n"))                       warning(paste("optim or nlminb returned message",
288                                       optimRes$message,"\n"))
289                   }                   }
290                   return(x)                   return(x)
291               })               })
# Line 325  Line 310 
310            function(object, ...)            function(object, ...)
311                .Call("lmer_fixef", object, PACKAGE = "Matrix"))                .Call("lmer_fixef", object, PACKAGE = "Matrix"))
312    
   
313  setMethod("fixef", signature(object = "lmer"),  setMethod("fixef", signature(object = "lmer"),
314            function(object, ...) object@fixed)            function(object, ...) object@fixed)
315    
# Line 407  Line 391 
391                                 row.names = ""))                                 row.names = ""))
392                }                }
393                cat("Random effects:\n")                cat("Random effects:\n")
394                show(VarCorr(object))                show(VarCorr(object, useScale = useScale))
395                ngrps <- lapply(object@flist, function(x) length(levels(x)))                ngrps <- lapply(object@flist, function(x) length(levels(x)))
396                cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))                cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))
397                cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))                cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
# Line 464  Line 448 
448            function(object, ...)            function(object, ...)
449        {        {
450            nc <- object@nc[-seq(along = object@Omega)]            nc <- object@nc[-seq(along = object@Omega)]
451            p <- nc[1] - 1            p <- abs(nc[1]) - 1
452            n <- nc[2]            n <- nc[2]
453            rep(n-p, p)            rep(n-p, p)
454        })        })
# Line 600  Line 584 
584            ci            ci
585        })        })
586    
 ##setMethod("param", signature(object = "lmer"),  
 ##          function(object, unconst = FALSE, ...) {  
 ##              .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")  
 ##          })  
   
587  setMethod("deviance", "mer",  setMethod("deviance", "mer",
588            function(object, REML = NULL, ...) {            function(object, REML = NULL, ...) {
589                .Call("lmer_factor", object, PACKAGE = "Matrix")                .Call("lmer_factor", object, PACKAGE = "Matrix")
# Line 735  Line 714 
714  setMethod("plot", signature(x = "lmer.coef"),  setMethod("plot", signature(x = "lmer.coef"),
715            function(x, y, ...)            function(x, y, ...)
716        {        {
717              if (require("lattice", quietly = TRUE)) {
718            varying <- unique(do.call("c",            varying <- unique(do.call("c",
719                                      lapply(x, function(el)                                      lapply(x, function(el)
720                                             names(el)[sapply(el,                                             names(el)[sapply(el,
# Line 749  Line 729 
729                                          list(y = as.name(varying[1]),                                          list(y = as.name(varying[1]),
730                                               x = as.name(varying[2])))), gf, ...),                                               x = as.name(varying[2])))), gf, ...),
731                   splom(~ gf | .grp, ...))                   splom(~ gf | .grp, ...))
732              }
733        })        })
734    
735  setMethod("plot", signature(x = "lmer.ranef"),  setMethod("plot", signature(x = "lmer.ranef"),
736            function(x, y, ...)            function(x, y, ...)
737        {        {
738              if (require("lattice", quietly = TRUE))
739            lapply(x, function(x) {            lapply(x, function(x) {
740                cn <- lapply(colnames(x), as.name)                cn <- lapply(colnames(x), as.name)
741                switch(min(ncol(x), 3),                switch(min(ncol(x), 3),
742                       qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),                           qqmath(eval(substitute(~ x,
743                       xyplot(eval(substitute(y ~ x, list(y = cn[[1]], x = cn[[2]]))),                                                  list(x = cn[[1]]))),
744                                    x, ...),
745                             xyplot(eval(substitute(y ~ x,
746                                                    list(y = cn[[1]],
747                                                         x = cn[[2]]))),
748                              x, ...),                              x, ...),
749                       splom(~ x, ...))                       splom(~ x, ...))
750            })            })

Legend:
Removed from v.776  
changed lines
  Added in v.777

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