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 751, Tue May 24 17:49:33 2005 UTC revision 752, Sun May 29 02:14:58 2005 UTC
# Line 28  Line 28 
28    function(maxIter = 50,    function(maxIter = 50,
29             msMaxIter = 50,             msMaxIter = 50,
30             tolerance = sqrt((.Machine$double.eps)),             tolerance = sqrt((.Machine$double.eps)),
31             niterEM = 20,             niterEM = 15,
32             msTol = sqrt(.Machine$double.eps),             msTol = sqrt(.Machine$double.eps),
33             msVerbose = getOption("verbose"),             msVerbose = getOption("verbose"),
34             PQLmaxIt = 20,             PQLmaxIt = 20,
# Line 126  Line 126 
126                   function(x, value)                   function(x, value)
127               {               {
128                   if (value$msMaxIter < 1) return(x)                   if (value$msMaxIter < 1) return(x)
                  st <- ccoef(x)         # starting values  
129                   nc <- x@nc                   nc <- x@nc
130                   nc <- nc[1:(length(nc) - 2)]                   nc <- nc[1:(length(nc) - 2)]
131                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
132                   fn <- function(pars) {                   fn <- function(pars)
133                       .Call("lmer_coefGets", x, pars, FALSE)                       deviance(.Call("lmer_coefGets", x, pars,
134                       #ccoef(x) <- pars                                      2, PACKAGE = "Matrix"),
135                       deviance(x, REML = value$REML)                                REML = value$REML)
                  }  
136                   gr <- if (value$analyticGradient)                   gr <- if (value$analyticGradient)
137                       function(pars) {                       function(pars) {
138                           if (!isTRUE(all.equal(pars, ccoef(x))))                           if (!isTRUE(all.equal(pars,
139                               .Call("lmer_coefGets", x, pars, FALSE)                                                 .Call("lmer_coef", x, 2))))
140                               #ccoef(x) <- pars                               .Call("lmer_coefGets", x, pars, 2)
141                           #grad <- gradient(x, REML = value$REML, unconst = TRUE)                           .Call("lmer_gradient", x, value$REML, 2)
                          gradient(x, REML = value$REML, unconst = FALSE)  
                          #grad[constr] <- -grad[constr]/pars[constr]  
                          #grad  
142                       } else NULL                       } else NULL
143                   if (require("port", quietly = TRUE)) {                   RV <- lapply(R.Version()[c("major", "minor")], as.numeric)
144                       optimRes <- portOptim(.Call("lmer_coef", x, FALSE),                   if ((RV$major > 1 && RV$minor >= 2.0) ||
145                         require("port", quietly = TRUE)) {
146                         optimRes <- nlminb(.Call("lmer_coef", x, 2),
147                                             fn, gr,                                             fn, gr,
148                                             lower = ifelse(constr, 1e-10, -Inf),                                             lower = ifelse(constr, 1e-10, -Inf),
149                                             control = list(maxit = value$msMaxIter,                                          control = list(iter.max = value$msMaxIter,
150                                             trace = as.integer(value$msVerbose)))                                             trace = as.integer(value$msVerbose)))
                      .Call("lmer_coefGets", x, optimRes$par, FALSE)  
151                   } else {                   } else {
152                       optimRes <- optim(ccoef(x), fn, gr,                       optimRes <- optim(.Call("lmer_coef", x, 2), fn, gr,
153                                         method = "L-BFGS-B",                                         method = "L-BFGS-B",
154                                         lower = ifelse(constr, 1e-10, -Inf),                                         lower = ifelse(constr, 1e-10, -Inf),
155                                         control = list(maxit = value$msMaxIter,                                         control = list(maxit = value$msMaxIter,
156                                         trace = as.integer(value$msVerbose)))                                         trace = as.integer(value$msVerbose)))
                      ccoef(x) <- optimRes$par  
157                   }                   }
158                     .Call("lmer_coefGets", x, optimRes$par, 2)
159                   if (optimRes$convergence != 0) {                   if (optimRes$convergence != 0) {
160                       warning(paste("optim returned message",optimRes$message,"\n"))                       warning(paste("optim returned message",optimRes$message,"\n"))
161                   }                   }
# Line 687  Line 683 
683    
684            if (method == "Laplace")            if (method == "Laplace")
685            {            {
686                if (require("port")) {                RV <- lapply(R.Version()[c("major", "minor")], as.numeric)
687                  if ((RV$major > 1 && RV$minor >= 2.0) ||
688                      require("port", quietly = TRUE)) {
689                    optimRes <-                    optimRes <-
690                        portOptim(fn = devLaplace,                        nlminb(c(fixef(obj),
691                                  par =                                 .Call("lmer_coef", obj, TRUE, PACKAGE = "Matrix")),
692                                  c(fixef(obj),                               devLaplace,
                                   .Call("lmer_coef",  
                                         obj,  
                                         TRUE,  
                                         PACKAGE = "Matrix")),  
                         ## WAS: coef(obj, unconst = TRUE)),  
                         method = "BFGS", #hessian = TRUE,  
693                          control = list(trace = getOption("verbose"),                          control = list(trace = getOption("verbose"),
694                          reltol = controlvals$msTol,                               iter.max = controlvals$msMaxIter))
                         maxit = controlvals$msMaxIter))  
695                    optpars <- optimRes$par                    optpars <- optimRes$par
696                } else {                } else {
697                    optimRes <-                    optimRes <-
# Line 800  Line 791 
791                val <- -deviance(object, REML = REML)/2                val <- -deviance(object, REML = REML)/2
792                nc <- object@nc[-seq(a = object@Omega)]                nc <- object@nc[-seq(a = object@Omega)]
793                attr(val, "nall") <- attr(val, "nobs") <- nc[2]                attr(val, "nall") <- attr(val, "nobs") <- nc[2]
794                attr(val, "df") <- nc[1] + length(ccoef(object))                attr(val, "df") <-
795                      nc[1] + length(.Call("lmer_coef", object, 2))
796                attr(val, "REML") <- REML                attr(val, "REML") <- REML
797                class(val) <- "logLik"                class(val) <- "logLik"
798                val                val
# Line 811  Line 803 
803                val <- object@glmmll                val <- object@glmmll
804                nc <- object@nc[-seq(a = object@Omega)]                nc <- object@nc[-seq(a = object@Omega)]
805                attr(val, "nall") <- attr(val, "nobs") <- nc[2]                attr(val, "nall") <- attr(val, "nobs") <- nc[2]
806                attr(val, "df") <- nc[1] + length(ccoef(object))                attr(val, "df") <-
807                ## attr(val, "REML") <- REML                    nc[1] + length(.Call("lmer_coef", object, 2))
808                class(val) <- "logLik"                class(val) <- "logLik"
809                val                val
810            })            })

Legend:
Removed from v.751  
changed lines
  Added in v.752

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