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 1166, Sun Jan 15 19:18:20 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,
# Line 236  Line 236 
236                                '\nUsing method = "PQL".\n')                                '\nUsing method = "PQL".\n')
237                }                }
238            }            }
239              if (method == "AGQ")
240                  stop('method = "AGQ" not yet implemented for supernodal representation')
241            ## create factor list for the random effects            ## create factor list for the random effects
242            bars <- findbars(formula[[3]])            bars <- findbars(formula[[3]])
243            names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))            names(bars) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
# Line 266  Line 268 
268            if (lmm) {            if (lmm) {
269                .Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")                .Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
270                LMEoptimize(mer) <- cv                LMEoptimize(mer) <- cv
271                return(mer)                return(new("lmer", mer,
272                             frame = if (model) frm else data.frame(),
273                             terms = glm.fit$terms))
274            }            }
275    
276            ## The rest of the function applies to generalized linear mixed models            ## The rest of the function applies to generalized linear mixed models
# Line 286  Line 290 
290    
291            GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")            GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
292            .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  
293            PQLpars <- c(fixef(mer),            PQLpars <- c(fixef(mer),
294                         .Call("mer_coef", mer, 2, PACKAGE = "Matrix"))                         .Call("mer_coef", mer, 2, PACKAGE = "Matrix"))
295              if (method == "PQL") {
296            .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")            .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")
297                  .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
298                  return(new("lmer", mer,
299                             frame = if (model) frm else data.frame(),
300                             terms = glm.fit$terms))
301              }
302    
303              fixInd <- seq(ncol(x))
304              ## pars[fixInd] == beta, pars[-fixInd] == theta
305            ## indicator of constrained parameters            ## indicator of constrained parameters
306            const <- c(rep(FALSE, length(fixInd)),            const <- c(rep(FALSE, length(fixInd)),
307                       unlist(lapply(mer@nc[seq(along = fl)],                       unlist(lapply(mer@nc[seq(along = fl)],
# Line 302  Line 313 
313            optimRes <-            optimRes <-
314                nlminb(PQLpars, devLaplace,                nlminb(PQLpars, devLaplace,
315                       lower = ifelse(const, 5e-10, -Inf),                       lower = ifelse(const, 5e-10, -Inf),
316                       control = list(trace = getOption("verbose"),                       control = list(trace = cv$msVerbose,
317                       iter.max = cv$msMaxIter))                       iter.max = cv$msMaxIter))
318            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
319            return(mer)            return(new("lmer", mer,
           deviance <- devAGQ(PQLpars, 1)  
   
 ### 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")  
           }  
   
           .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")  
           loglik[] <- -deviance/2  
           new("lmer", mer,  
320                frame = if (model) frm else data.frame(),                frame = if (model) frm else data.frame(),
321                terms = glm.fit$terms,                       terms = glm.fit$terms))
               assign = attr(glm.fit$x, "assign"),  
               call = match.call(), family = family,  
               logLik = loglik, fixed = fxd)  
       })  
   
322    
323  ## Extract the permutation        })
 setAs("mer", "pMatrix", function(from)  
       .Call("mer_pMatrix", from, PACKAGE = "Matrix"))  
324    
325  ## Extract the L matrix  ## Extract the L matrix
326  setAs("mer", "dtCMatrix", function(from)  setAs("mer", "dtCMatrix", function(from)
# Line 369  Line 334 
334  ## Extract the random effects  ## Extract the random effects
335  setMethod("ranef", signature(object = "mer"),  setMethod("ranef", signature(object = "mer"),
336            function(object, ...)            function(object, ...)
337                .Call("mer_ranef", object, PACKAGE = "Matrix")                new("lmer.ranef", .Call("mer_ranef", object, PACKAGE = "Matrix"))
338            )            )
339    
340  ## Optimization for mer objects  ## Optimization for mer objects
# Line 381  Line 346 
346                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
347                   fn <- function(pars)                   fn <- function(pars)
348                       deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))                       deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
349                   gr <- if (value$analyticGradient)                   gr <- if (value$gradient)
350                       function(pars) {                       function(pars) {
351                           if (!isTRUE(all.equal(pars,                           if (!isTRUE(all.equal(pars,
352                                                 .Call("mer_coef", x,                                                 .Call("mer_coef", x,
# Line 516  Line 481 
481                fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")                fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")
482                useScale <- object@useScale                useScale <- object@useScale
483                corF <- vcov(object)@factors$correlation                corF <- vcov(object)@factors$correlation
484                DF <- getFixDF(object)                #DF <- getFixDF(object)
485                coefs <- cbind(fcoef, corF@sd, DF)                coefs <- cbind(fcoef, corF@sd)#, DF)
486                dimnames(coefs) <-                dimnames(coefs) <-
487                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))                    list(names(fcoef), c("Estimate", "Std. Error"))#, "DF"))
488                digits <- max(3, getOption("digits") - 2)                digits <- max(3, getOption("digits") - 2)
489                REML <- object@method == "REML"                REML <- object@method == "REML"
490                llik <- logLik(object, REML)                llik <- logLik(object, REML)
# Line 571  Line 536 
536                if (nrow(coefs) > 0) {                if (nrow(coefs) > 0) {
537                    if (useScale) {                    if (useScale) {
538                        stat <- coefs[,1]/coefs[,2]                        stat <- coefs[,1]/coefs[,2]
539                        pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)                        #pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
540                        nms <- colnames(coefs)                        nms <- colnames(coefs)
541                        coefs <- cbind(coefs, stat, pval)                        coefs <- cbind(coefs, stat) #, pval)
542                        colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")                        colnames(coefs) <- c(nms, "t value")#, "Pr(>|t|)")
543                    } else {                    } else {
544                        coefs <- coefs[, 1:2, drop = FALSE]                        coefs <- coefs[, 1:2, drop = FALSE]
545                        stat <- coefs[,1]/coefs[,2]                        stat <- coefs[,1]/coefs[,2]
# Line 584  Line 549 
549                        colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")                        colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
550                    }                    }
551                    cat("\nFixed effects:\n")                    cat("\nFixed effects:\n")
552                    printCoefmat(coefs, tst.ind = 4, zap.ind = 3)                    printCoefmat(coefs, zap.ind = 3)#, tst.ind = 4)
553                    rn <- rownames(coefs)                    rn <- rownames(coefs)
554                    if (!is.null(corF)) {                    if (!is.null(corF)) {
555                        p <- ncol(corF)                        p <- ncol(corF)
# Line 660  Line 625 
625            dots <- list(...)            dots <- list(...)
626            modp <- logical(0)            modp <- logical(0)
627            if (length(dots))            if (length(dots))
628                modp <- sapply(dots, inherits, "mer") | sapply(dots, inherits, "lm")                modp <- sapply(dots, is, "mer") | sapply(dots, is, "lm")
629            if (any(modp)) {              # multiple models - form table            if (any(modp)) {              # multiple models - form table
630                opts <- dots[!modp]                opts <- dots[!modp]
631                mods <- c(list(object), dots[modp])                mods <- c(list(object), dots[modp])
# Line 697  Line 662 
662                return(val)                return(val)
663            } else {            } else {
664                foo <- object                foo <- object
665                foo@status["factored"] <- FALSE                #foo@status["factored"] <- FALSE
666                .Call("mer_factor", foo, PACKAGE="Matrix")                #.Call("mer_factor", foo, PACKAGE="Matrix")
667                dfr <- getFixDF(foo)                #dfr <- getFixDF(foo)
668                ss <- foo@rXy^2                ss <- foo@rXy^2
669                ssr <- exp(foo@devComp["logryy2"])                ssr <- exp(foo@devComp["logryy2"])
670                ss <- ss[seq(along = dfr)]                names(ss) <- object@cnames[[".fixed"]]
671                names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]                asgn <- attr(foo@X, "assign")
               asgn <- foo@assign  
672                terms <- foo@terms                terms <- foo@terms
673                nmeffects <- attr(terms, "term.labels")                nmeffects <- attr(terms, "term.labels")
674                if ("(Intercept)" %in% names(ss))                if ("(Intercept)" %in% names(ss))
# Line 721  Line 685 
685                    list(nmeffects,                    list(nmeffects,
686  #                       c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))  #                       c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
687                         c("Df", "Sum Sq", "Mean Sq"))                         c("Df", "Sum Sq", "Mean Sq"))
688                if ("(Intercept)" %in% nmeffects) table <- table[-1,]                if ("(Intercept)" %in% nmeffects) table <- table[-match("(Intercept)", nmeffects), ]
689                attr(table, "heading") <- "Analysis of Variance Table"                attr(table, "heading") <- "Analysis of Variance Table"
690                class(table) <- c("anova", "data.frame")                class(table) <- c("anova", "data.frame")
691                table                table
# Line 754  Line 718 
718            )            )
719    
720  setMethod("summary", signature(object = "mer"),  setMethod("summary", signature(object = "mer"),
721            function(object, ...)            function(object, ...) object
           .NotYetImplemented()  
722            )            )
723    
724  setMethod("update", signature(object = "mer"),  setMethod("update", signature(object = "mer"),
# Line 767  Line 730 
730  simss <- function(fm0, fma, nsim)  simss <- function(fm0, fma, nsim)
731  {  {
732      ysim <- simulate(fm0, nsim)      ysim <- simulate(fm0, nsim)
733      cv <- list(analyticGradient = FALSE, msMaxIter = 200:200,      cv <- list(gradient = FALSE, msMaxIter = 200:200,
734                 msVerbose = 0:0)                 msVerbose = 0:0)
735      sapply(ysim, function(yy) {      sapply(ysim, function(yy) {
736          .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")          .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")
# Line 778  Line 741 
741                Ha = fma@devComp[["logryy2"]]))                Ha = fma@devComp[["logryy2"]]))
742      })      })
743  }  }
744    
745    ## Some leftover code from the old AGQ method in lmer.
746    if (FALSE) {
747    ### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs
748    ### AGQ for nc > 1 first.
749        fxd <- PQLpars[fixInd]
750        loglik <- logLik(mer)
751    
752        if (method %in% c("Laplace", "AGQ")) {
753            nAGQ <- 1
754            if (method == "AGQ") {    # determine nAGQ at PQL estimates
755                dev11 <- devAGQ(PQLpars, 11)
756                ## FIXME: Should this be an absolute or a relative tolerance?
757                devTol <- sqrt(.Machine$double.eps) * abs(dev11)
758                for (nAGQ in c(9, 7, 5, 3, 1))
759                    if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break
760                nAGQ <- nAGQ + 2
761                if (gVerb)
762                    cat(paste("Using", nAGQ, "quadrature points per column\n"))
763            }
764            obj <- function(pars)
765                .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
766            optimRes <-
767                nlminb(PQLpars, obj,
768                       lower = ifelse(const, 5e-10, -Inf),
769                       control = list(trace = getOption("verbose"),
770                       iter.max = cv$msMaxIter))
771            optpars <- optimRes$par
772            if (optimRes$convergence != 0)
773                warning("nlminb failed to converge")
774            deviance <- optimRes$objective
775            if (gVerb)
776                cat(paste("convergence message", optimRes$message, "\n"))
777            fxd[] <- optpars[fixInd]  ## preserve the names
778            .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
779        }
780    
781        .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
782        loglik[] <- -deviance/2
783    }

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

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