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 878, Sun Aug 28 13:57:39 2005 UTC revision 879, Tue Aug 30 14:35:05 2005 UTC
# Line 833  Line 833 
833            print(reMat, quote = FALSE)            print(reMat, quote = FALSE)
834        })        })
835    
836  setMethod("mcmcsamp", signature(obj = "lmer"),  setMethod("mcmcsamp", signature(object = "lmer"),
837            function(obj, nsamp = 1, verbose = FALSE, saveb = FALSE,            function(object, n = 1, verbose = FALSE, saveb = FALSE,
838                     trans = TRUE, ...)                     trans = TRUE, ...)
839        {        {
840            if (obj@family$family == "gaussian" &&            if (object@family$family == "gaussian" &&
841                obj@family$link == "identity") {                object@family$link == "identity") {
842                glmer <- FALSE                glmer <- FALSE
843                ans <- .Call("lmer_MCMCsamp", obj, saveb, nsamp, trans,                ans <- .Call("lmer_MCMCsamp", object, saveb, n, trans,
844                              PACKAGE = "Matrix")                              PACKAGE = "Matrix")
845            } else {            } else {
846                glmer <- TRUE                glmer <- TRUE
# Line 848  Line 848 
848                    warning("trans option not currently allowed for generalized models")                    warning("trans option not currently allowed for generalized models")
849                trans <- FALSE                trans <- FALSE
850                ## Check arguments                ## Check arguments
851                if (length(obj@Omega) > 1 || obj@nc[1] > 1)                if (length(object@Omega) > 1 || object@nc[1] > 1)
852                    stop("mcmcsamp currently defined for glmm models with only one variance component")                    stop("mcmcsamp currently defined for glmm models with only one variance component")
853                cv <- Matrix:::lmerControl()                cv <- Matrix:::lmerControl()
854                if (verbose) cv$msVerbose <- 1                if (verbose) cv$msVerbose <- 1
855                family <- obj@family                family <- object@family
856                frm <- obj@frame                frm <- object@frame
857    
858                ## recreate model matrices                ## recreate model matrices
859                fixed.form <- Matrix:::nobars(obj@call$formula)                fixed.form <- Matrix:::nobars(object@call$formula)
860                if (inherits(fixed.form, "name")) # RHS is empty - use a constant                if (inherits(fixed.form, "name")) # RHS is empty - use a constant
861                    fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))                    fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
862                glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE,                glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE,
863                               y = TRUE)                               y = TRUE)
864                x <- glm.fit$x                x <- glm.fit$x
865                y <- as.double(glm.fit$y)                y <- as.double(glm.fit$y)
866                bars <- Matrix:::findbars(obj@call$formula[[3]])                bars <- Matrix:::findbars(object@call$formula[[3]])
867                random <-                random <-
868                    lapply(bars,                    lapply(bars,
869                           function(x) list(model.matrix(eval(substitute(~term,                           function(x) list(model.matrix(eval(substitute(~term,
# Line 872  Line 872 
872                                            eval(substitute(as.factor(fac)[,drop = TRUE],                                            eval(substitute(as.factor(fac)[,drop = TRUE],
873                                                            list(fac = x[[3]])), frm)))                                                            list(fac = x[[3]])), frm)))
874                names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))                names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
875                if (any(names(random) != names(obj@flist)))                if (any(names(random) != names(object@flist)))
876                    random <- random[names(obj@flist)]                    random <- random[names(object@flist)]
877                mmats <- c(lapply(random, "[[", 1),                mmats <- c(lapply(random, "[[", 1),
878                           .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))                           .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
879                mer <- as(obj, "mer")                mer <- as(object, "mer")
880    
881                ## establish the GS object and the ans matrix                ## establish the GS object and the ans matrix
882                eta <- glm.fit$linear.predictors # perhaps later change this to obj@fitted?                eta <- glm.fit$linear.predictors # perhaps later change this to object@fitted?
883                wts <- glm.fit$prior.weights                wts <- glm.fit$prior.weights
884                wtssqr <- wts * wts                wtssqr <- wts * wts
885                offset <- glm.fit$offset                offset <- glm.fit$offset
# Line 893  Line 893 
893                LMEopt <- getAnywhere("LMEoptimize<-")                LMEopt <- getAnywhere("LMEoptimize<-")
894                doLMEopt <- quote(LMEopt(x = mer, value = cv))                doLMEopt <- quote(LMEopt(x = mer, value = cv))
895                GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")                GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")
896                fixed <- obj@fixed                fixed <- object@fixed
897                varc <- .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")                varc <- .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")
898                b <- .Call("lmer_ranef", mer, PACKAGE = "Matrix")                b <- .Call("lmer_ranef", mer, PACKAGE = "Matrix")
899                ans <- .Call("glmer_MCMCsamp", GSpt, b, fixed, varc, saveb, nsamp,                ans <- .Call("glmer_MCMCsamp", GSpt, b, fixed, varc, saveb, n,
900                             PACKAGE = "Matrix")                             PACKAGE = "Matrix")
901                .Call("glmer_finalize", GSpt, PACKAGE = "Matrix");                .Call("glmer_finalize", GSpt, PACKAGE = "Matrix");
902            }            }
903            gnms <- names(obj@flist)            gnms <- names(object@flist)
904            cnms <- obj@cnames            cnms <- object@cnames
905            ff <- fixef(obj)            ff <- fixef(object)
906            colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",            colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
907                      unlist(lapply(seq(along = gnms),                      unlist(lapply(seq(along = gnms),
908                                    function(i)                                    function(i)
# Line 929  Line 929 
929  rWishart <- function(n, df, invScal)  rWishart <- function(n, df, invScal)
930    .Call("Matrix_rWishart", n, df, invScal)    .Call("Matrix_rWishart", n, df, invScal)
931    
932  setMethod("simulate", signature(obj = "lmer"),  setMethod("simulate", signature(object = "lmer"),
933            function(obj, nsamp = 1, verbose = FALSE, ...)            function(object, n = 1, ...)
934        {        {
935            family <- obj@family            family <- object@family
936            if (family$family != "gaussian" ||            if (family$family != "gaussian" ||
937                family$link != "identity")                family$link != "identity")
938                stop("simulation of generalized linear mixed models not implemented yet")                stop("simulation of generalized linear mixed models not implemented yet")
939            ## create the mean from the fixed effects            ## create the mean from the fixed effects
940            frm <- obj@frame            frm <- object@frame
941            n <- nrow(frm)            fixed.form <- Matrix:::nobars(object@call$formula)
           fixed.form <- Matrix:::nobars(obj@call$formula)  
942            if (inherits(fixed.form, "name")) # RHS is empty - use a constant            if (inherits(fixed.form, "name")) # RHS is empty - use a constant
943                fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))                fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
944            glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE, y = TRUE)            glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE, y = TRUE)
945            fxd <- drop(glm.fit$x %*% fixef(obj))            fxd <- unname(drop(glm.fit$x %*% fixef(object)))
946            bars <- Matrix:::findbars(obj@call$formula[[3]])            bars <- Matrix:::findbars(object@call$formula[[3]])
947            random <-            random <-
948                lapply(bars,                lapply(bars,
949                       function(x) list(model.matrix(eval(substitute(~term,                       function(x) list(model.matrix(eval(substitute(~term,
# Line 954  Line 953 
953                                                        list(fac = x[[3]])), frm)))                                                        list(fac = x[[3]])), frm)))
954            names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))            names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
955            ## re-order the random effects pairs if necessary            ## re-order the random effects pairs if necessary
956            if (any(names(random) != names(obj@flist)))            if (any(names(random) != names(object@flist)))
957                random <- random[names(obj@flist)]                random <- random[names(object@flist)]
958            rmmats <- c(lapply(random, "[[", 1))            mer <- as(object, "mer")
959            vc <- VarCorr(obj)            ## similate the linear predictors
960            ans <- fxd + matrix(rnorm(n * nsamp, sd = vc@scale), nr = n)            lpred <- .Call("lmer_simulate", mer, n, fxd,
961            ans                           lapply(random, "[[", 1), TRUE, PACKAGE = "Matrix")
962              ## add per-observation noise term
963              lpred + rnorm(prod(dim(lpred)),
964                            sd = .Call("lmer_sigma", mer,
965                            object@method == "REML", PACKAGE = "Matrix"))
966          })
967    
968    simulate2 <- function(object, n = 1, ...)
969    {
970        family <- object@family
971        if (family$family != "gaussian" ||
972            family$link != "identity")
973            stop("simulation of generalized linear mixed models not implemented yet")
974    
975        ## create the mean from the fixed effects
976        frm <- object@frame
977        fixed.form <- Matrix:::nobars(object@call$formula)
978        if (inherits(fixed.form, "name")) # RHS is empty - use a constant
979            fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
980        glm.fit <- glm(eval(fixed.form), family, frm, x = TRUE, y = TRUE)
981        lpred <- matrix(glm.fit$x %*% fixef(object), nr = nrow(frm), nc = n)
982    
983        ## Create the random effects model matrices
984        bars <- Matrix:::findbars(object@call$formula[[3]])
985        random <-
986            lapply(bars,
987                   function(x) list(model.matrix(eval(substitute(~term,
988                                                                 list(term=x[[2]]))),
989                                                 frm),
990                                    eval(substitute(as.factor(fac)[,drop = TRUE],
991                                                    list(fac = x[[3]])), frm)))
992        names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
993        ## re-order the random effects pairs if necessary
994        flist <- object@flist
995        if (any(names(random) != names(flist)))
996            random <- random[names(flist)]
997        mmats <- lapply(random, "[[", 1)
998    
999        ## simulate the random effects
1000        scale <- .Call("lmer_sigma", object, object@method == "REML",
1001                       PACKAGE = "Matrix")
1002        Omega <- object@Omega
1003        re <- lapply(seq(along = Omega),
1004                     function(i) {
1005                         om <- Omega[[i]]
1006                         nr <- nrow(om)
1007                         nlev <- length(levels(flist[[i]]))
1008                         scale * array(solve(chol(new("dpoMatrix", Dim = dim(om),
1009                                                      uplo = "U", x = c(om))),
1010                                             matrix(rnorm(nr * n * nlev),
1011                                                    nr = nr))@x, c(nr, n, nlev))
1012        })        })
1013        ## apply the random effects
1014        for (j in seq(along = Omega)) {
1015            for (i in 1:nrow(lpred))
1016            lpred[i,] <- lpred[i,] + mmats[[j]][i,] %*% re[[j]][, , as.integer(flist[[j]])[i]]
1017        }
1018        ## add per-observation noise term
1019        lpred <- lpred + rnorm(prod(dim(lpred)), sd = scale)
1020        attr(lpred, "re") <- re
1021        lpred
1022    }
1023    
1024    refdist <- function(fm1, fm2, n, ...)
1025    {
1026        frm <- fm1@frame
1027        newy <- simulate(fm2, n)
1028        ref <- numeric(n)
1029        f1 <- eval(fm1@call$formula)
1030        f2 <- eval(fm2@call$formula)
1031        nm <- as.character(as.name(f1[[2]]))
1032        for (j in 1:n) {
1033            frm[[nm]] <- newy[, j]
1034            ref[j] <- deviance(lmer(f2, frm)) - deviance(lmer(f1, frm))
1035        }
1036        attr(ref, "observed") <- deviance(fm2) - deviance(fm1)
1037        ref
1038    }

Legend:
Removed from v.878  
changed lines
  Added in v.879

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