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 1165, Sat Jan 14 23:41:55 2006 UTC revision 1175, Tue Jan 17 23:41:36 2006 UTC
# Line 91  Line 91 
91             niterEM = 15,             niterEM = 15,
92             EMverbose = getOption("verbose"),             EMverbose = getOption("verbose"),
93             PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead             PQLmaxIt = 30,# FIXME: unused; PQL currently uses 'maxIter' instead
94             analyticGradient = TRUE,             gradient = TRUE,
95             analyticHessian = FALSE # unused _FIXME_             Hessian = FALSE # unused _FIXME_
96             )             )
97  {  {
98      list(maxIter = as.integer(maxIter),      list(maxIter = as.integer(maxIter),
# Line 103  Line 103 
103           niterEM = as.integer(niterEM),           niterEM = as.integer(niterEM),
104           EMverbose = as.logical(EMverbose),           EMverbose = as.logical(EMverbose),
105           PQLmaxIt = as.integer(PQLmaxIt),           PQLmaxIt = as.integer(PQLmaxIt),
106           analyticGradient = as.logical(analyticGradient),           gradient = as.logical(gradient),
107           analyticHessian = as.logical(analyticHessian))           Hessian = as.logical(Hessian))
108  }  }
109    
110  setMethod("coef", signature(object = "lmer"),  rWishart <- function(n, df, invScal)
111        .Call("Matrix_rWishart", n, df, invScal, PACKAGE = "Matrix")
112    
113    setMethod("coef", signature(object = "mer"),
114            function(object, ...)            function(object, ...)
115        {        {
116            fef <- data.frame(rbind(object@fixed), check.names = FALSE)             fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
117            ref <- as(ranef(object), "list")             ref <- ranef(object)
118            names(ref) <- names(object@flist)             val <- lapply(ref, function(x) fef[rep(1, nrow(x)),,drop = FALSE])
           val <- lapply(ref, function(x) fef[rep(1, nrow(x)),])  
119            for (i in seq(a = val)) {            for (i in seq(a = val)) {
120                refi <- ref[[i]]                refi <- ref[[i]]
121                row.names(val[[i]]) <- row.names(refi)                row.names(val[[i]]) <- row.names(refi)
122                if (!all(names(refi) %in% names(fef)))                 nmsi <- colnames(refi)
123                   if (!all(nmsi %in% names(fef)))
124                    stop("unable to align random and fixed effects")                    stop("unable to align random and fixed effects")
125                val[[i]][ , names(refi)] <- val[[i]][ , names(refi)] + refi                 for (nm in nmsi) val[[i]][[nm]] <- val[[i]][[nm]] + refi[,nm]
126            }            }
127            new("lmer.coef", val)             val
128        })        })
129    
130  ## setMethod("plot", signature(x = "lmer.coef"),  ## setMethod("plot", signature(x = "lmer.coef"),
# Line 143  Line 146 
146  ##                  splom(~ gf | .grp, ...))  ##                  splom(~ gf | .grp, ...))
147  ##       })  ##       })
148    
149  ## setMethod("plot", signature(x = "lmer.ranef"),  setMethod("plot", signature(x = "lmer.ranef"),
150  ##           function(x, y, ...)            function(x, y, ...)
151  ##       {        {
152  ##           lapply(x, function(x) {            lapply(x, function(x) {
153  ##               cn <- lapply(colnames(x), as.name)                cn <- lapply(colnames(x), as.name)
154  ##               switch(min(ncol(x), 3),                switch(min(ncol(x), 3),
155  ##                      qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),                       qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...),
156  ##                      xyplot(eval(substitute(y ~ x,                       xyplot(eval(substitute(y ~ x,
157  ##                                             list(y = cn[[1]],                                              list(y = cn[[1]],
158  ##                                                  x = cn[[2]]))), x, ...),                                                   x = cn[[2]]))), x, ...),
159  ##                      splom(~ x, ...))                       splom(~ x, ...))
160  ##           })            })
161  ##       })        })
162    
163  setMethod("with", signature(data = "lmer"),  setMethod("with", signature(data = "lmer"),
164            function(data, expr, ...) {            function(data, expr, ...) {
# Line 169  Line 172 
172  setMethod("terms", signature(x = "lmer"),  setMethod("terms", signature(x = "lmer"),
173            function(x, ...) x@terms)            function(x, ...) x@terms)
174    
 rWishart <- function(n, df, invScal)  
     .Call("Matrix_rWishart", n, df, invScal, PACKAGE = "Matrix")  
   
175    
176  setMethod("lmer", signature(formula = "formula"),  setMethod("lmer", signature(formula = "formula"),
177            function(formula, data, family,            function(formula, data, family = gaussian,
178                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
179                     control = list(), start,                     control = list(), start,
180                     subset, weights, na.action, offset,                     subset, weights, na.action, offset, contrasts = NULL,
181                     model = TRUE, x = FALSE, y = FALSE , ...)                     model = TRUE, x = TRUE, y = TRUE , ...)
182        {        {
183            ## match and check parameters            ## match and check parameters
184            if (length(formula) < 3) stop("formula must be a two-sided formula")            if (length(formula) < 3) stop("formula must be a two-sided formula")
# Line 188  Line 188 
188            ## that frame.  Otherwise missing values in the grouping factors            ## that frame.  Otherwise missing values in the grouping factors
189            ## cause inconsistent numbers of observations.            ## cause inconsistent numbers of observations.
190            mf <- match.call()            mf <- match.call()
191            m <- match(c("family", "data", "subset", "weights",            m <- match(c("data", "subset", "weights",
192                         "na.action", "offset"), names(mf), 0)                         "na.action", "offset"), names(mf), 0)
193            mf <- fe <- mf[c(1, m)]            mf <- mf[c(1, m)]
194            frame.form <- subbars(formula) # substitute `+' for `|'            frame.form <- subbars(formula) # substitute `+' for `|'
195            fixed.form <- nobars(formula)  # remove any terms with `|'            fixed.form <- nobars(formula)  # remove any terms with `|'
196            if (inherits(fixed.form, "name")) # RHS is empty - use a constant            if (inherits(fixed.form, "name")) # RHS is empty - use a constant
# Line 199  Line 199 
199    
200            ## evaluate a model frame for fixed and random effects            ## evaluate a model frame for fixed and random effects
201            mf$formula <- frame.form            mf$formula <- frame.form
           mf$family <- NULL  
202            mf$drop.unused.levels <- TRUE            mf$drop.unused.levels <- TRUE
203            mf[[1]] <- as.name("model.frame")            mf[[1]] <- as.name("model.frame")
204            frm <- eval(mf, parent.frame())            fe <- mf
205              mf <- eval(mf, parent.frame())
206    
207            ## fit a glm model to the fixed formula            ## get the terms for the fixed-effects only
208            fe$formula <- fixed.form            fe$formula <- fixed.form
209            fe$subset <- NULL             # subset has already been created in call to data.frame            fe <- eval(fe, parent.frame())
210            fe$data <- frm            mt <- attr(fe, "terms")   # allow model.frame to update them
211            fe$x <- fe$model <- fe$y <- TRUE            ## response vector
212            fe[[1]] <- as.name("glm")            Y <- model.response(mf, "numeric")
213            glm.fit <- eval(fe, parent.frame())            ## avoid problems with 1D arrays, but keep names
214            x <- glm.fit$x            if(length(dim(Y)) == 1) {
215            y <- as.double(glm.fit$y)                nm <- rownames(Y)
216            family <- glm.fit$family                dim(Y) <- NULL
217                  if(!is.null(nm)) names(Y) <- nm
218              }
219              ## null model support
220              X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts) else matrix(,NROW(Y),0)
221    
222              weights <- model.weights(mf)
223              offset <- model.offset(mf)
224              ## check weights and offset
225              if( !is.null(weights) && any(weights < 0) )
226                  stop("negative weights not allowed")
227              if(!is.null(offset) && length(offset) != NROW(Y))
228                  stop(gettextf("number of offsets is %d should equal %d (number of observations)", length(offset), NROW(Y)), domain = NA)
229              if (is.null(weights)) weights <- rep.int(1, length(Y))
230              if (is.null(offset)) offset <- numeric(length(Y))
231    
232              ## fit a glm model to the fixed formula
233    ##           fe$formula <- fixed.form
234    ##           fe$subset <- NULL             # subset has already been created in call to data.frame
235    ##           fe$data <- frm
236    ##           fe$x <- fe$model <- fe$y <- TRUE
237    ##           fe[[1]] <- as.name("glm")
238    ##           glmFit <- eval(fe, parent.frame())
239    ##           x <- glmFit$x
240    ##           y <- as.double(glmFit$y)
241    
242              if(is.character(family))
243                  family <- get(family, mode = "function", envir = parent.frame())
244              if(is.function(family)) family <- family()
245              if(is.null(family$family)) {
246                  print(family)
247                  stop("'family' not recognized")
248              }
249            ## check for a linear mixed model            ## check for a linear mixed model
250            lmm <- family$family == "gaussian" && family$link == "identity"            lmm <- family$family == "gaussian" && family$link == "identity"
251            if (lmm) { # linear mixed model            if (lmm) { # linear mixed model
# Line 236  Line 267 
267                                '\nUsing method = "PQL".\n')                                '\nUsing method = "PQL".\n')
268                }                }
269            }            }
270              if (method == "AGQ")
271                  stop('method = "AGQ" not yet implemented for supernodal representation')
272            ## create factor list for the random effects            ## create factor list for the random effects
273            bars <- findbars(formula[[3]])            bars <- findbars(formula[[3]])
274            names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))            names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
275            fl <- lapply(bars,            fl <- lapply(bars,
276                         function(x)                         function(x)
277                         eval(substitute(as.factor(fac)[,drop = TRUE],                         eval(substitute(as.factor(fac)[,drop = TRUE],
278                                         list(fac = x[[3]])), frm))                                         list(fac = x[[3]])), mf))
279            ## order factor list by decreasing number of levels            ## order factor list by decreasing number of levels
280            nlev <- sapply(fl, function(x) length(levels(x)))            nlev <- sapply(fl, function(x) length(levels(x)))
281            if (any(diff(nlev) > 0)) {            if (any(diff(nlev) > 0)) {
# Line 254  Line 287 
287            Ztl <- lapply(bars, function(x)            Ztl <- lapply(bars, function(x)
288                          t(model.matrix(eval(substitute(~ expr,                          t(model.matrix(eval(substitute(~ expr,
289                                                         list(expr = x[[2]]))),                                                         list(expr = x[[2]]))),
290                                         frm)))                                         mf)))
291            ## Create the mixed-effects representation (mer) object            ## Create the mixed-effects representation (mer) object
292            mer <- .Call("mer_create", fl,            mer <- .Call("mer_create", fl,
293                         .Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"),                         .Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"),
294                         x, y, method, sapply(Ztl, nrow),                         X, Y, method, sapply(Ztl, nrow),
295                         c(lapply(Ztl, rownames), list(.fixed = colnames(x))),                         c(lapply(Ztl, rownames), list(.fixed = colnames(X))),
296                         !(family$family %in% c("binomial", "poisson")),                         !(family$family %in% c("binomial", "poisson")),
297                         match.call(), family,                         match.call(), family,
298                         PACKAGE = "Matrix")                         PACKAGE = "Matrix")
299            if (lmm) {            if (lmm) {
300                .Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")                .Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
301                LMEoptimize(mer) <- cv                LMEoptimize(mer) <- cv
302                return(mer)                return(new("lmer", mer,
303                             frame = if (model) mf else data.frame(),
304                             terms = mt))
305            }            }
306    
307            ## The rest of the function applies to generalized linear mixed models            ## The rest of the function applies to generalized linear mixed models
308            gVerb <- getOption("verbose")            gVerb <- getOption("verbose")
309            eta <- glm.fit$linear.predictors            eta <- glm.fit(X, Y, weights = weights, offset = offset, family = family,
310            wts <- glm.fit$prior.weights                           intercept = attr(mt, "intercept") > 0))$fitted.values
311            wtssqr <- wts * wts            wtssqr <- weights * weights
           offset <- glm.fit$offset  
312            if (is.null(offset)) offset <- numeric(length(eta))            if (is.null(offset)) offset <- numeric(length(eta))
313            linkinv <- quote(family$linkinv(eta))            linkinv <- quote(family$linkinv(eta))
314            mu.eta <- quote(family$mu.eta(eta))            mu.eta <- quote(family$mu.eta(eta))
# Line 286  Line 320 
320    
321            GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")            GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
322            .Call("glmer_PQL", GSpt, PACKAGE = "Matrix")  # obtain PQL estimates            .Call("glmer_PQL", GSpt, PACKAGE = "Matrix")  # obtain PQL estimates
           fixInd <- seq(ncol(x))  
           ## pars[fixInd] == beta, pars[-fixInd] == theta  
323            PQLpars <- c(fixef(mer),            PQLpars <- c(fixef(mer),
324                         .Call("mer_coef", mer, 2, PACKAGE = "Matrix"))                         .Call("mer_coef", mer, 2, PACKAGE = "Matrix"))
325              if (method == "PQL") {
326            .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")            .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")
327                  .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
328                  return(new("lmer", mer,
329                             frame = if (model) mf else data.frame(),
330                             terms = mt))
331              }
332    
333              fixInd <- seq(ncol(x))
334              ## pars[fixInd] == beta, pars[-fixInd] == theta
335            ## indicator of constrained parameters            ## indicator of constrained parameters
336            const <- c(rep(FALSE, length(fixInd)),            const <- c(rep(FALSE, length(fixInd)),
337                       unlist(lapply(mer@nc[seq(along = fl)],                       unlist(lapply(mer@nc[seq(along = fl)],
# Line 302  Line 343 
343            optimRes <-            optimRes <-
344                nlminb(PQLpars, devLaplace,                nlminb(PQLpars, devLaplace,
345                       lower = ifelse(const, 5e-10, -Inf),                       lower = ifelse(const, 5e-10, -Inf),
346                       control = list(trace = getOption("verbose"),                       control = list(trace = cv$msVerbose,
347                       iter.max = cv$msMaxIter))                       iter.max = cv$msMaxIter))
348            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
349            return(mer)            return(new("lmer", mer,
350            deviance <- devAGQ(PQLpars, 1)                       frame = if (model) mf else data.frame(),
351                         terms = mt))
 ### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs  
 ### AGQ for nc > 1 first.  
           fxd <- PQLpars[fixInd]  
           loglik <- logLik(mer)  
   
           if (method %in% c("Laplace", "AGQ")) {  
               nAGQ <- 1  
               if (method == "AGQ") {    # determine nAGQ at PQL estimates  
                   dev11 <- devAGQ(PQLpars, 11)  
                   ## FIXME: Should this be an absolute or a relative tolerance?  
                   devTol <- sqrt(.Machine$double.eps) * abs(dev11)  
                   for (nAGQ in c(9, 7, 5, 3, 1))  
                       if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break  
                   nAGQ <- nAGQ + 2  
                   if (gVerb)  
                       cat(paste("Using", nAGQ, "quadrature points per column\n"))  
               }  
               obj <- function(pars)  
                   .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")  
               optimRes <-  
                   nlminb(PQLpars, obj,  
                          lower = ifelse(const, 5e-10, -Inf),  
                          control = list(trace = getOption("verbose"),  
                          iter.max = cv$msMaxIter))  
               optpars <- optimRes$par  
               if (optimRes$convergence != 0)  
                   warning("nlminb failed to converge")  
               deviance <- optimRes$objective  
               if (gVerb)  
                   cat(paste("convergence message", optimRes$message, "\n"))  
               fxd[] <- optpars[fixInd]  ## preserve the names  
               .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")  
           }  
352    
           .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")  
           loglik[] <- -deviance/2  
           new("lmer", mer,  
               frame = if (model) frm else data.frame(),  
               terms = glm.fit$terms,  
               assign = attr(glm.fit$x, "assign"),  
               call = match.call(), family = family,  
               logLik = loglik, fixed = fxd)  
353        })        })
354    
   
 ## Extract the permutation  
 setAs("mer", "pMatrix", function(from)  
       .Call("mer_pMatrix", from, PACKAGE = "Matrix"))  
   
355  ## Extract the L matrix  ## Extract the L matrix
356  setAs("mer", "dtCMatrix", function(from)  setAs("mer", "dtCMatrix", function(from)
357        .Call("mer_dtCMatrix", from, PACKAGE = "Matrix"))        .Call("mer_dtCMatrix", from, PACKAGE = "Matrix"))
# Line 369  Line 364 
364  ## Extract the random effects  ## Extract the random effects
365  setMethod("ranef", signature(object = "mer"),  setMethod("ranef", signature(object = "mer"),
366            function(object, ...)            function(object, ...)
367                .Call("mer_ranef", object, PACKAGE = "Matrix")                new("lmer.ranef", .Call("mer_ranef", object, PACKAGE = "Matrix"))
368            )            )
369    
370  ## Optimization for mer objects  ## Optimization for mer objects
# Line 381  Line 376 
376                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
377                   fn <- function(pars)                   fn <- function(pars)
378                       deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))                       deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
379                   gr <- if (value$analyticGradient)                   gr <- if (value$gradient)
380                       function(pars) {                       function(pars) {
381                           if (!isTRUE(all.equal(pars,                           if (!isTRUE(all.equal(pars,
382                                                 .Call("mer_coef", x,                                                 .Call("mer_coef", x,
# Line 413  Line 408 
408            function(object, n = 1, verbose = FALSE, saveb = FALSE,            function(object, n = 1, verbose = FALSE, saveb = FALSE,
409                     trans = TRUE, ...)                     trans = TRUE, ...)
410        {        {
411              family <- object@family
412              lmm <- family$family == "gaussian" && family$link == "identity"
413              if (!lmm)
414                  stop("mcmcsamp for GLMMs not yet implemented in supernodal representation")
415            ans <- t(.Call("mer_MCMCsamp", object, saveb, n,            ans <- t(.Call("mer_MCMCsamp", object, saveb, n,
416                           trans, PACKAGE = "Matrix"))                           trans, PACKAGE = "Matrix"))
417            attr(ans, "mcpar") <- as.integer(c(1, n, 1))            attr(ans, "mcpar") <- as.integer(c(1, n, 1))
# Line 516  Line 515 
515                fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")                fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")
516                useScale <- object@useScale                useScale <- object@useScale
517                corF <- vcov(object)@factors$correlation                corF <- vcov(object)@factors$correlation
518                DF <- getFixDF(object)                #DF <- getFixDF(object)
519                coefs <- cbind(fcoef, corF@sd, DF)                coefs <- cbind(fcoef, corF@sd)#, DF)
520                dimnames(coefs) <-                dimnames(coefs) <-
521                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))                    list(names(fcoef), c("Estimate", "Std. Error"))#, "DF"))
522                digits <- max(3, getOption("digits") - 2)                digits <- max(3, getOption("digits") - 2)
523                REML <- object@method == "REML"                REML <- object@method == "REML"
524                llik <- logLik(object, REML)                llik <- logLik(object, REML)
# Line 571  Line 570 
570                if (nrow(coefs) > 0) {                if (nrow(coefs) > 0) {
571                    if (useScale) {                    if (useScale) {
572                        stat <- coefs[,1]/coefs[,2]                        stat <- coefs[,1]/coefs[,2]
573                        pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)                        #pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
574                        nms <- colnames(coefs)                        nms <- colnames(coefs)
575                        coefs <- cbind(coefs, stat, pval)                        coefs <- cbind(coefs, stat) #, pval)
576                        colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")                        colnames(coefs) <- c(nms, "t value")#, "Pr(>|t|)")
577                    } else {                    } else {
578                        coefs <- coefs[, 1:2, drop = FALSE]                        coefs <- coefs[, 1:2, drop = FALSE]
579                        stat <- coefs[,1]/coefs[,2]                        stat <- coefs[,1]/coefs[,2]
# Line 584  Line 583 
583                        colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")                        colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
584                    }                    }
585                    cat("\nFixed effects:\n")                    cat("\nFixed effects:\n")
586                    printCoefmat(coefs, tst.ind = 4, zap.ind = 3)                    printCoefmat(coefs, zap.ind = 3)#, tst.ind = 4)
587                    rn <- rownames(coefs)                    rn <- rownames(coefs)
588                    if (!is.null(corF)) {                    if (!is.null(corF)) {
589                        p <- ncol(corF)                        p <- ncol(corF)
# Line 644  Line 643 
643            if (useScale)            if (useScale)
644                sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix")                sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix")
645            sc2 <- sc * sc            sc2 <- sc * sc
646            ans <- lapply(x@Omega, function(el) {            cnames <- x@cnames
647                el <- as(sc2 * solve(el), "dpoMatrix")            ans <- x@Omega
648              for (i in seq(a=ans)) {
649                  el <- as(sc2 * solve(ans[[i]]), "dpoMatrix")
650                  el@Dimnames <- list(cnames[[i]], cnames[[i]])
651                el@factors$correlation <- as(el, "correlation")                el@factors$correlation <- as(el, "correlation")
652                el                ans[[i]] <- el
653            })            }
654            attr(ans, "sc") <- sc            attr(ans, "sc") <- sc
655            ans            ans
656        })        })
# Line 660  Line 662 
662            dots <- list(...)            dots <- list(...)
663            modp <- logical(0)            modp <- logical(0)
664            if (length(dots))            if (length(dots))
665                modp <- sapply(dots, inherits, "mer") | sapply(dots, inherits, "lm")                modp <- sapply(dots, is, "mer") | sapply(dots, is, "lm")
666            if (any(modp)) {              # multiple models - form table            if (any(modp)) {              # multiple models - form table
667                opts <- dots[!modp]                opts <- dots[!modp]
668                mods <- c(list(object), dots[modp])                mods <- c(list(object), dots[modp])
# Line 697  Line 699 
699                return(val)                return(val)
700            } else {            } else {
701                foo <- object                foo <- object
702                foo@status["factored"] <- FALSE                #foo@status["factored"] <- FALSE
703                .Call("mer_factor", foo, PACKAGE="Matrix")                #.Call("mer_factor", foo, PACKAGE="Matrix")
704                dfr <- getFixDF(foo)                #dfr <- getFixDF(foo)
705                ss <- foo@rXy^2                ss <- foo@rXy^2
706                ssr <- exp(foo@devComp["logryy2"])                ssr <- exp(foo@devComp["logryy2"])
707                ss <- ss[seq(along = dfr)]                names(ss) <- object@cnames[[".fixed"]]
708                names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]                asgn <- attr(foo@X, "assign")
               asgn <- foo@assign  
709                terms <- foo@terms                terms <- foo@terms
710                nmeffects <- attr(terms, "term.labels")                nmeffects <- attr(terms, "term.labels")
711                if ("(Intercept)" %in% names(ss))                if ("(Intercept)" %in% names(ss))
# Line 721  Line 722 
722                    list(nmeffects,                    list(nmeffects,
723  #                       c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))  #                       c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
724                         c("Df", "Sum Sq", "Mean Sq"))                         c("Df", "Sum Sq", "Mean Sq"))
725                if ("(Intercept)" %in% nmeffects) table <- table[-1,]                if ("(Intercept)" %in% nmeffects) table <- table[-match("(Intercept)", nmeffects), ]
726                attr(table, "heading") <- "Analysis of Variance Table"                attr(table, "heading") <- "Analysis of Variance Table"
727                class(table) <- c("anova", "data.frame")                class(table) <- c("anova", "data.frame")
728                table                table
# Line 754  Line 755 
755            )            )
756    
757  setMethod("summary", signature(object = "mer"),  setMethod("summary", signature(object = "mer"),
758            function(object, ...)            function(object, ...) object
           .NotYetImplemented()  
759            )            )
760    
761  setMethod("update", signature(object = "mer"),  setMethod("update", signature(object = "mer"),
# Line 767  Line 767 
767  simss <- function(fm0, fma, nsim)  simss <- function(fm0, fma, nsim)
768  {  {
769      ysim <- simulate(fm0, nsim)      ysim <- simulate(fm0, nsim)
770      cv <- list(analyticGradient = FALSE, msMaxIter = 200:200,      cv <- list(gradient = FALSE, msMaxIter = 200:200,
771                 msVerbose = 0:0)                 msVerbose = 0:0)
772      sapply(ysim, function(yy) {      sapply(ysim, function(yy) {
773          .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")          .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")
# Line 778  Line 778 
778                Ha = fma@devComp[["logryy2"]]))                Ha = fma@devComp[["logryy2"]]))
779      })      })
780  }  }
781    
782    ## Some leftover code from the old AGQ method in lmer.
783    if (FALSE) {
784    ### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs
785    ### AGQ for nc > 1 first.
786        fxd <- PQLpars[fixInd]
787        loglik <- logLik(mer)
788    
789        if (method %in% c("Laplace", "AGQ")) {
790            nAGQ <- 1
791            if (method == "AGQ") {    # determine nAGQ at PQL estimates
792                dev11 <- devAGQ(PQLpars, 11)
793                ## FIXME: Should this be an absolute or a relative tolerance?
794                devTol <- sqrt(.Machine$double.eps) * abs(dev11)
795                for (nAGQ in c(9, 7, 5, 3, 1))
796                    if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break
797                nAGQ <- nAGQ + 2
798                if (gVerb)
799                    cat(paste("Using", nAGQ, "quadrature points per column\n"))
800            }
801            obj <- function(pars)
802                .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
803            optimRes <-
804                nlminb(PQLpars, obj,
805                       lower = ifelse(const, 5e-10, -Inf),
806                       control = list(trace = getOption("verbose"),
807                       iter.max = cv$msMaxIter))
808            optpars <- optimRes$par
809            if (optimRes$convergence != 0)
810                warning("nlminb failed to converge")
811            deviance <- optimRes$objective
812            if (gVerb)
813                cat(paste("convergence message", optimRes$message, "\n"))
814            fxd[] <- optpars[fixInd]  ## preserve the names
815            .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
816        }
817    
818        .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
819        loglik[] <- -deviance/2
820    }

Legend:
Removed from v.1165  
changed lines
  Added in v.1175

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