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 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 
70             analyticGradient = TRUE,             analyticGradient = TRUE,
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),
81           analyticHessian=analyticHessian,           analyticGradient = as.logical(analyticGradient),
82           analyticGradient=analyticGradient)           analyticHessian = as.logical(analyticHessian))
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))
187            linkinv <- quote(family$linkinv(eta))            linkinv <- quote(family$linkinv(eta))
188              mu <- eval(linkinv)
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))  
           {  
               mu <- eval(linkinv) # family$linkinv(eta)  
               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

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