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 1171, Mon Jan 16 16:59:46 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 413  Line 378 
378            function(object, n = 1, verbose = FALSE, saveb = FALSE,            function(object, n = 1, verbose = FALSE, saveb = FALSE,
379                     trans = TRUE, ...)                     trans = TRUE, ...)
380        {        {
381              family <- object@family
382              lmm <- family$family == "gaussian" && family$link == "identity"
383              if (!lmm)
384                  stop("mcmcsamp for GLMMs not yet implemented in supernodal representation")
385            ans <- t(.Call("mer_MCMCsamp", object, saveb, n,            ans <- t(.Call("mer_MCMCsamp", object, saveb, n,
386                           trans, PACKAGE = "Matrix"))                           trans, PACKAGE = "Matrix"))
387            attr(ans, "mcpar") <- as.integer(c(1, n, 1))            attr(ans, "mcpar") <- as.integer(c(1, n, 1))
# Line 516  Line 485 
485                fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")                fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")
486                useScale <- object@useScale                useScale <- object@useScale
487                corF <- vcov(object)@factors$correlation                corF <- vcov(object)@factors$correlation
488                DF <- getFixDF(object)                #DF <- getFixDF(object)
489                coefs <- cbind(fcoef, corF@sd, DF)                coefs <- cbind(fcoef, corF@sd)#, DF)
490                dimnames(coefs) <-                dimnames(coefs) <-
491                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))                    list(names(fcoef), c("Estimate", "Std. Error"))#, "DF"))
492                digits <- max(3, getOption("digits") - 2)                digits <- max(3, getOption("digits") - 2)
493                REML <- object@method == "REML"                REML <- object@method == "REML"
494                llik <- logLik(object, REML)                llik <- logLik(object, REML)
# Line 571  Line 540 
540                if (nrow(coefs) > 0) {                if (nrow(coefs) > 0) {
541                    if (useScale) {                    if (useScale) {
542                        stat <- coefs[,1]/coefs[,2]                        stat <- coefs[,1]/coefs[,2]
543                        pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)                        #pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
544                        nms <- colnames(coefs)                        nms <- colnames(coefs)
545                        coefs <- cbind(coefs, stat, pval)                        coefs <- cbind(coefs, stat) #, pval)
546                        colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")                        colnames(coefs) <- c(nms, "t value")#, "Pr(>|t|)")
547                    } else {                    } else {
548                        coefs <- coefs[, 1:2, drop = FALSE]                        coefs <- coefs[, 1:2, drop = FALSE]
549                        stat <- coefs[,1]/coefs[,2]                        stat <- coefs[,1]/coefs[,2]
# Line 584  Line 553 
553                        colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")                        colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
554                    }                    }
555                    cat("\nFixed effects:\n")                    cat("\nFixed effects:\n")
556                    printCoefmat(coefs, tst.ind = 4, zap.ind = 3)                    printCoefmat(coefs, zap.ind = 3)#, tst.ind = 4)
557                    rn <- rownames(coefs)                    rn <- rownames(coefs)
558                    if (!is.null(corF)) {                    if (!is.null(corF)) {
559                        p <- ncol(corF)                        p <- ncol(corF)
# Line 644  Line 613 
613            if (useScale)            if (useScale)
614                sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix")                sc <- .Call("mer_sigma", x, REML, PACKAGE = "Matrix")
615            sc2 <- sc * sc            sc2 <- sc * sc
616            ans <- lapply(x@Omega, function(el) {            cnames <- x@cnames
617                el <- as(sc2 * solve(el), "dpoMatrix")            ans <- x@Omega
618              for (i in seq(a=ans)) {
619                  el <- as(sc2 * solve(ans[[i]]), "dpoMatrix")
620                  el@Dimnames <- list(cnames[[i]], cnames[[i]])
621                el@factors$correlation <- as(el, "correlation")                el@factors$correlation <- as(el, "correlation")
622                el                ans[[i]] <- el
623            })            }
624            attr(ans, "sc") <- sc            attr(ans, "sc") <- sc
625            ans            ans
626        })        })
# Line 660  Line 632 
632            dots <- list(...)            dots <- list(...)
633            modp <- logical(0)            modp <- logical(0)
634            if (length(dots))            if (length(dots))
635                modp <- sapply(dots, inherits, "mer") | sapply(dots, inherits, "lm")                modp <- sapply(dots, is, "mer") | sapply(dots, is, "lm")
636            if (any(modp)) {              # multiple models - form table            if (any(modp)) {              # multiple models - form table
637                opts <- dots[!modp]                opts <- dots[!modp]
638                mods <- c(list(object), dots[modp])                mods <- c(list(object), dots[modp])
# Line 697  Line 669 
669                return(val)                return(val)
670            } else {            } else {
671                foo <- object                foo <- object
672                foo@status["factored"] <- FALSE                #foo@status["factored"] <- FALSE
673                .Call("mer_factor", foo, PACKAGE="Matrix")                #.Call("mer_factor", foo, PACKAGE="Matrix")
674                dfr <- getFixDF(foo)                #dfr <- getFixDF(foo)
675                ss <- foo@rXy^2                ss <- foo@rXy^2
676                ssr <- exp(foo@devComp["logryy2"])                ssr <- exp(foo@devComp["logryy2"])
677                ss <- ss[seq(along = dfr)]                names(ss) <- object@cnames[[".fixed"]]
678                names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]                asgn <- attr(foo@X, "assign")
               asgn <- foo@assign  
679                terms <- foo@terms                terms <- foo@terms
680                nmeffects <- attr(terms, "term.labels")                nmeffects <- attr(terms, "term.labels")
681                if ("(Intercept)" %in% names(ss))                if ("(Intercept)" %in% names(ss))
# Line 721  Line 692 
692                    list(nmeffects,                    list(nmeffects,
693  #                       c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))  #                       c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
694                         c("Df", "Sum Sq", "Mean Sq"))                         c("Df", "Sum Sq", "Mean Sq"))
695                if ("(Intercept)" %in% nmeffects) table <- table[-1,]                if ("(Intercept)" %in% nmeffects) table <- table[-match("(Intercept)", nmeffects), ]
696                attr(table, "heading") <- "Analysis of Variance Table"                attr(table, "heading") <- "Analysis of Variance Table"
697                class(table) <- c("anova", "data.frame")                class(table) <- c("anova", "data.frame")
698                table                table
# Line 754  Line 725 
725            )            )
726    
727  setMethod("summary", signature(object = "mer"),  setMethod("summary", signature(object = "mer"),
728            function(object, ...)            function(object, ...) object
           .NotYetImplemented()  
729            )            )
730    
731  setMethod("update", signature(object = "mer"),  setMethod("update", signature(object = "mer"),
# Line 767  Line 737 
737  simss <- function(fm0, fma, nsim)  simss <- function(fm0, fma, nsim)
738  {  {
739      ysim <- simulate(fm0, nsim)      ysim <- simulate(fm0, nsim)
740      cv <- list(analyticGradient = FALSE, msMaxIter = 200:200,      cv <- list(gradient = FALSE, msMaxIter = 200:200,
741                 msVerbose = 0:0)                 msVerbose = 0:0)
742      sapply(ysim, function(yy) {      sapply(ysim, function(yy) {
743          .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")          .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")
# Line 778  Line 748 
748                Ha = fma@devComp[["logryy2"]]))                Ha = fma@devComp[["logryy2"]]))
749      })      })
750  }  }
751    
752    ## Some leftover code from the old AGQ method in lmer.
753    if (FALSE) {
754    ### FIXME: For nf == 1 change this to an AGQ evaluation.  Needs
755    ### AGQ for nc > 1 first.
756        fxd <- PQLpars[fixInd]
757        loglik <- logLik(mer)
758    
759        if (method %in% c("Laplace", "AGQ")) {
760            nAGQ <- 1
761            if (method == "AGQ") {    # determine nAGQ at PQL estimates
762                dev11 <- devAGQ(PQLpars, 11)
763                ## FIXME: Should this be an absolute or a relative tolerance?
764                devTol <- sqrt(.Machine$double.eps) * abs(dev11)
765                for (nAGQ in c(9, 7, 5, 3, 1))
766                    if (abs(dev11 - devAGQ(PQLpars, nAGQ - 2)) > devTol) break
767                nAGQ <- nAGQ + 2
768                if (gVerb)
769                    cat(paste("Using", nAGQ, "quadrature points per column\n"))
770            }
771            obj <- function(pars)
772                .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")
773            optimRes <-
774                nlminb(PQLpars, obj,
775                       lower = ifelse(const, 5e-10, -Inf),
776                       control = list(trace = getOption("verbose"),
777                       iter.max = cv$msMaxIter))
778            optpars <- optimRes$par
779            if (optimRes$convergence != 0)
780                warning("nlminb failed to converge")
781            deviance <- optimRes$objective
782            if (gVerb)
783                cat(paste("convergence message", optimRes$message, "\n"))
784            fxd[] <- optpars[fixInd]  ## preserve the names
785            .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
786        }
787    
788        .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")
789        loglik[] <- -deviance/2
790    }

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

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