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

branches/trunk-lme4/R/lmer.R revision 767, Tue Jun 7 23:05:00 2005 UTC pkg/R/lmer.R revision 769, Wed Jun 8 23:47:45 2005 UTC
# Line 7  Line 7 
7      ((i - 1) * i)/2 + j      ((i - 1) * i)/2 + j
8  }  }
9    
10  ##Dhalf <- function(from) {  # Return the pairs of expressions separated by vertical bars
11  ##    D <- from@D  
12  ##    nf <- length(D)  findbars <- function(term)
13  ##    Gp <- from@Gp  {
14  ##    res <- array(0, rep(Gp[nf+1],2))      if (is.name(term) || is.numeric(term)) return(NULL)
15  ##    for (i in 1:nf) {      if (term[[1]] == as.name("(")) return(findbars(term[[2]]))
16  ##        DD <- D[[i]]      if (!is.call(term)) stop("term must be of class call")
17  ##        dd <- dim(DD)      if (term[[1]] == as.name('|')) return(term)
18  ##        for (k in 1:dd[3]) {      if (length(term) == 2) return(findbars(term[[2]]))
19  ##            mm <- array(DD[ , , k], dd[1:2])      c(findbars(term[[2]]), findbars(term[[3]]))
20  ##            base <- Gp[i] + (k - 1)*dd[1]  }
21  ##            res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)  
22  ##        }  # Return the formula omitting the pairs of expressions separated by vertical bars
23  ##    }  
24  ##    res  nobars <- function(term)
25  ##}  {
26        # FIXME: is the is.name in the condition redundant?
27        #   A name won't satisfy the first condition.
28        if (!('|' %in% all.names(term)) || is.name(term)) return(term)
29        if (is.call(term) && term[[1]] == as.name('|')) return(NULL)
30        if (length(term) == 2) {
31            nb <- nobars(term[[2]])
32            if (is.null(nb)) return(NULL)
33            term[[2]] <- nb
34            return(term)
35        }
36        nb2 <- nobars(term[[2]])
37        nb3 <- nobars(term[[3]])
38        if (is.null(nb2)) return(nb3)
39        if (is.null(nb3)) return(nb2)
40        term[[2]] <- nb2
41        term[[3]] <- nb3
42        term
43    }
44    
45    # Substitute the '+' function for the '|' function
46    
47    subbars <- function(term)
48    {
49        if (is.name(term) || is.numeric(term)) return(term)
50        if (length(term) == 2) {
51            term[[2]] <- subbars(term[[2]])
52            return(term)
53        }
54        stopifnot(length(term) == 3)
55        if (is.call(term) && term[[1]] == as.name('|')) term[[1]] <- as.name('+')
56        term[[2]] <- subbars(term[[2]])
57        term[[3]] <- subbars(term[[3]])
58        term
59    }
60    
61  lmerControl <-                            # Control parameters for lmer  lmerControl <-                            # Control parameters for lmer
62    function(maxIter = 50,    function(maxIter = 50,
63             msMaxIter = 50,             msMaxIter = 200,
64             tolerance = sqrt((.Machine$double.eps)),             tolerance = sqrt((.Machine$double.eps)),
65             niterEM = 15,             niterEM = 15,
66             msTol = sqrt(.Machine$double.eps),             msTol = sqrt(.Machine$double.eps),
67             msVerbose = getOption("verbose"),             msVerbose = getOption("verbose"),
68             PQLmaxIt = 20,             PQLmaxIt = 30,
69             EMverbose = getOption("verbose"),             EMverbose = getOption("verbose"),
70             analyticGradient = TRUE,             analyticGradient = TRUE,
71             analyticHessian=FALSE)             analyticHessian=FALSE)
# Line 93  Line 127 
127            family <- glm.fit$family            family <- glm.fit$family
128            x <- glm.fit$x            x <- glm.fit$x
129            y <- as.double(glm.fit$y)            y <- as.double(glm.fit$y)
130              family <- glm.fit$family
131    
132            ## evaluate a model frame for fixed and random effects            ## evaluate a model frame for fixed and random effects
133            mf$formula <- frame.form            mf$formula <- frame.form
# Line 129  Line 164 
164                LMEoptimize(mer) <- cv                LMEoptimize(mer) <- cv
165                fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")                fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")
166                return(new("lmer",                return(new("lmer",
167                           mer, frame = if (frame) frm else data.frame(),                           mer,
                          terms = glm.fit$terms,  
168                           assign = attr(x, "assign"),                           assign = attr(x, "assign"),
169                           call = match.call(),                           call = match.call(),
170                             family = family, fitted = fits,
171                             fixed = fixef(mer),
172                             frame = if (model) frm else data.frame(),
173                             logLik = logLik(mer),
174                           residuals = unname(model.response(frm) - fits),                           residuals = unname(model.response(frm) - fits),
175                           fitted = fits))                           terms = glm.fit$terms))
176            }            }
177    
178            ## The rest of the function applies to generalized linear mixed models            ## The rest of the function applies to generalized linear mixed models
# Line 144  Line 182 
182            offset <- glm.fit$offset            offset <- glm.fit$offset
183            if (is.null(offset)) offset <- numeric(length(eta))            if (is.null(offset)) offset <- numeric(length(eta))
184    
185            dev.resids <- quote(family$dev.resids(y, mu, wts))            dev.resids <- quote(family$dev.resids(y, mu, wts*wts))
186            linkinv <- quote(family$linkinv(eta))            linkinv <- quote(family$linkinv(eta))
187            mu.eta <- quote(family$mu.eta(eta))            mu.eta <- quote(family$mu.eta(eta))
188            variance <- quote(family$variance(mu))            variance <- quote(family$variance(mu))
# Line 199  Line 237 
237            fixInd <- seq(along = fixef(mer))            fixInd <- seq(along = fixef(mer))
238            env <- environment()            env <- environment()
239            bhat <- function(pars) { #pars[fixInd] == beta, pars[-fixInd] == theta            bhat <- function(pars) { #pars[fixInd] == beta, pars[-fixInd] == theta
240                if (cv$analyticHessian) {  #              if (cv$analyticHessian) {
241                    .Call("glmer_bhat_iterate", pars, cv$tolerance,                    .Call("glmer_bhat_iterate", pars, cv$tolerance,
242                          env, PACKAGE = "Matrix")                          env, PACKAGE = "Matrix")
243                } else {  ##              } else {
244                    etaold <- off <- eta <- drop(glm.fit$x %*% pars[fixInd]) + offset  ##                  etaold <- off <- eta <- drop(glm.fit$x %*% pars[fixInd]) + offset
245                    .Call("lmer_coefGets", rdobj, as.double(pars[-fixInd]),  ##                  .Call("lmer_coefGets", rdobj, as.double(pars[-fixInd]),
246                          2, PACKAGE = "Matrix")  ##                        2, PACKAGE = "Matrix")
247                    niter <- 20  ##                  niter <- 20
248                    conv <- FALSE  ##                  conv <- FALSE
249                    for (iter in seq(length = niter)) {  ##                  for (iter in seq(length = niter)) {
250                        mu <- eval(linkinv)  ##                      mu <- eval(linkinv)
251                        dmu.deta <- eval(mu.eta)  ##                      dmu.deta <- eval(mu.eta)
252                        .Call("glmer_weight_matrix_list", unwtd,  ##                      .Call("glmer_weight_matrix_list", unwtd,
253                              wts * dmu.deta / sqrt(eval(variance)),  ##                            wts * dmu.deta / sqrt(eval(variance)),
254                              eta - off + (y - mu)/dmu.deta, wtd,  ##                            eta - off + (y - mu)/dmu.deta, wtd,
255                              PACKAGE="Matrix")  ##                            PACKAGE="Matrix")
256                        .Call("lmer_update_mm", rdobj, wtd,  ##                      .Call("lmer_update_mm", rdobj, wtd,
257                              PACKAGE="Matrix")  ##                            PACKAGE="Matrix")
258                        eta <- off + .Call("lmer_fitted", rdobj,  ##                      eta <- off + .Call("lmer_fitted", rdobj,
259                                           unwtd, TRUE, PACKAGE = "Matrix")  ##                                         unwtd, TRUE, PACKAGE = "Matrix")
260                        if (max(abs(eta - etaold)) <  ##                      if (max(abs(eta - etaold)) <
261                            (0.1 + max(abs(eta))) * cv$tolerance) {  ##                          (0.1 + max(abs(eta))) * cv$tolerance) {
262                            conv <- TRUE  ##                          conv <- TRUE
263                            break  ##                          break
264                        }  ##                      }
265                        etaold <- eta  ##                      etaold <- eta
266                    }  ##                  }
267                    if (!conv) warning("iterations for bhat did not converge")  ##                  if (!conv) warning("iterations for bhat did not converge")
268                }  ##              }
269                ## bhat doesn't really need to return anything, we  ##              invisible(eval(linkinv))
               ## just want the side-effect of modifying rdobj  
               ## In particular, we are interested in  
               ## ranef(rdobj) and rdobj@bVar (?). But  
               ## the mu-scale response will be useful for log-lik  
               ## calculations later, so return them anyway  
   
               invisible(eval(linkinv))  
270            }            }
271    
   
           ## function that calculates log likelihood (the thing that  
           ## needs to get evaluated at each Gauss-Hermite location)  
   
           ## log scale ? worry about details later, get the pieces in place  
   
           ## this is for the Laplace approximation only. GH is more  
           ## complicated  
   
272            devLaplace <- function(pars) {            devLaplace <- function(pars) {
273                mu <- bhat(pars)                mu <- bhat(pars)
   
               ## GLM family log likelihood (upto constant ?)(log scale)  
               ## FIXME: need to adjust for sigma^2 for appropriate models (not trivial!)  
   
               ## Keep everything on (log) likelihood scale  
   
               ## log lik from observations given fixed and random effects  
               ## get deviance, then multiply by -1/2 (since deviance = -2 log lik)  
274                ## FIXME: This actually returns half the deviance.                ## FIXME: This actually returns half the deviance.
275                dev <- sum(family$dev.resids(y = glm.fit$y,                dev <- sum(eval(dev.resids))/2
                                            mu = mu, wt = wts^2))/2  
276                lap <- .Call("glmer_Laplace_devComp", rdobj, PACKAGE = "Matrix")                lap <- .Call("glmer_Laplace_devComp", rdobj, PACKAGE = "Matrix")
277                if (gVerb) cat(sprintf("  %#11g %#11g\n", dev, lap))                if (gVerb) cat(sprintf("  %#11g %#11g\n", dev, lap))
278                dev - lap                dev - lap
# Line 299  Line 312 
312                    cat(paste("optim convergence code",                    cat(paste("optim convergence code",
313                              optimRes$convergence, "\n"))                              optimRes$convergence, "\n"))
314                    cat("Fixed effects:\n")                    cat("Fixed effects:\n")
                   #print(fixef(mer))  ## no longer correct values  
315                    print(optimRes$par[fixInd])                    print(optimRes$par[fixInd])
316                    cat("(box constrained) variance coefficients:\n")                    cat("(box constrained) variance coefficients:\n")
                   #print(.Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))  
317                    print(optimRes$par[-fixInd])                    print(optimRes$par[-fixInd])
318                }                }
319            } else {            } else {
# Line 311  Line 322 
322            }            }
323    
324            loglik <- -devLaplace(optpars)            loglik <- -devLaplace(optpars)
325              attributes(loglik) <- attributes(logLik(mer))
326            ff <- optpars[fixInd]            ff <- optpars[fixInd]
327            names(ff) <- names(fixef(mer))            names(ff) <- names(fixef(mer))
328            .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")            .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
329            new("glmer", new("lmer", mer, frame = frm, terms = glm.fit$terms,            new("lmer", mer, frame = frm, terms = glm.fit$terms,
330                assign = attr(glm.fit$x, "assign"), call = match.call()),                assign = attr(glm.fit$x, "assign"), call = match.call(),
331                family = family, glmmll = loglik, fixed = ff)                family = family, logLik = loglik, fixed = ff)
332        })        })
333    
334  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
# Line 381  Line 393 
393                val[-length(val)]                val[-length(val)]
394            })            })
395    
396  setMethod("fixef", signature(object = "glmer"),  setMethod("fixef", signature(object = "lmer"),
397            function(object, ...) {            function(object, ...) object@fixed)
               object@fixed  
           })  
398    
399  setMethod("VarCorr", signature(x = "lmer"),  setMethod("VarCorr", signature(x = "lmer"),
400            function(x, REML = TRUE, useScale = TRUE, ...) {            function(x, REML = TRUE, useScale = TRUE, ...) {
# Line 405  Line 415 
415    
416  setMethod("summary", signature(object = "lmer"),  setMethod("summary", signature(object = "lmer"),
417            function(object, ...)            function(object, ...)
418            new("summary.lmer", object, useScale = TRUE,            new("summary.lmer", object,
               showCorrelation = TRUE,  
               method = object@method,  
               family = gaussian(),  
               logLik = logLik(object),  
               fixed = fixef(object)))  
   
 ## FIXME: glmm-s with scale not handled  
 setMethod("summary", signature(object = "glmer"),  
           function(object, ...)  
           new("summary.lmer", object, useScale = FALSE,  
419                showCorrelation = TRUE,                showCorrelation = TRUE,
420                method = object@method,                useScale = !((object@family)$family %in% c("binomial", "poisson"))))
               family = object@family,  
               logLik = logLik(object),  
               fixed = fixef(object)))  
421    
422  setMethod("show", signature(object = "lmer"),  setMethod("show", signature(object = "lmer"),
423            function(object)            function(object)
424            show(new("summary.lmer", object, useScale = TRUE,            show(new("summary.lmer", object,
425                     showCorrelation = FALSE,                     showCorrelation = FALSE,
426                     method = object@method,                     useScale = !((object@family)$family %in% c("binomial", "poisson")))))
                    family = gaussian(),  
                    logLik = logLik(object),  
                    fixed = fixef(object)))  
           )  
   
 setMethod("show", signature(object = "glmer"),  
           function(object)  
           show(new("summary.lmer", object, useScale = FALSE,  
                    showCorrelation = FALSE,  
                    method = object@method,  
                    family = object@family,  
                    logLik = logLik(object),  
                    fixed = fixef(object)))  
           )  
427    
428  setMethod("show", "summary.lmer",  setMethod("show", "summary.lmer",
429            function(object) {            function(object) {
# Line 565  Line 548 
548                val                val
549            })            })
550    
551  setMethod("logLik", signature(object="glmer"),  setMethod("logLik", signature(object="lmer"),
552            function(object, ...) {            function(object, ...) object@logLik)
               val <- object@glmmll  
               nc <- object@nc[-seq(a = object@Omega)]  
               attr(val, "nall") <- attr(val, "nobs") <- nc[2]  
               attr(val, "df") <- nc[1] +  
                   length(.Call("lmer_coef", object, 0, PACKAGE = "Matrix"))  
               class(val) <- "logLik"  
               val  
           })  
553    
554  setMethod("anova", signature(object = "lmer"),  setMethod("anova", signature(object = "lmer"),
555            function(object, ...)            function(object, ...)
# Line 692  Line 667 
667            ci            ci
668        })        })
669    
670  setMethod("param", signature(object = "lmer"),  ##setMethod("param", signature(object = "lmer"),
671            function(object, unconst = FALSE, ...) {  ##          function(object, unconst = FALSE, ...) {
672                .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")  ##              .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
673            })  ##          })
674    
675  setMethod("deviance", "mer",  setMethod("deviance", "mer",
676            function(object, REML = NULL, ...) {            function(object, REML = NULL, ...) {
# Line 706  Line 681 
681            })            })
682    
683    
684  setMethod("deviance", "glmer",  setMethod("deviance", "lmer",
685            function(object, ...) {            function(object, ...) -2 * c(object@logLik))
686                -2 * object@glmmll  
           })  
687    
688  setMethod("chol", signature(x = "lmer"),  setMethod("chol", signature(x = "lmer"),
689            function(x, pivot = FALSE, LINPACK = pivot) {            function(x, pivot = FALSE, LINPACK = pivot) {
# Line 811  Line 785 
785  setMethod("coef", signature(object = "lmer"),  setMethod("coef", signature(object = "lmer"),
786            function(object, ...)            function(object, ...)
787        {        {
788            fef <- data.frame(rbind(fixef(object)), check.names = FALSE)            fef <- data.frame(rbind(object@fixed), check.names = FALSE)
789            ref <- as(ranef(object), "list")            ref <- as(ranef(object), "list")
790            names(ref) <- names(object@flist)            names(ref) <- names(object@flist)
791            val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])            val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])
# Line 905  Line 879 
879            if (!useScale) reMat <- reMat[-nrow(reMat),]            if (!useScale) reMat <- reMat[-nrow(reMat),]
880            print(reMat, quote = FALSE)            print(reMat, quote = FALSE)
881        })        })
882    

Legend:
Removed from v.767  
changed lines
  Added in v.769

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