Re: [R] mcmcsamp() in lmer

From: Douglas Bates <dmbates_at_gmail.com>
Date: Sun 22 Jan 2006 - 09:13:37 EST

On 1/21/06, Spencer Graves <spencer.graves@pdf.com> wrote:
> Dear Prof. Gelman:
>
> Thanks for providing such a simple, replicable example. I got the
> same results you describe under Windows XP Pro with:
> > sessionInfo()
> R version 2.2.1, 2005-12-20, i386-pc-mingw32
>
> attached base packages:
> [1] "methods" "stats" "graphics" "grDevices" "utils" "datasets"
> [7] "base"
>
> other attached packages:
> lme4 lattice Matrix
> "0.995-1" "0.12-11" "0.995-1"
>
> Unfortunately, I'm missing something in my attempts to move beyond
> this. First, I tried "traceback()", which gave me a standard
> cryptogram, which I couldn't decypher. Then I typed "mcmcsamp" and
> learned that it consists solely of a call to
> 'standardGeneric("mcmcsamp")'. The help file for 'standardGeneric' sent
> me to "?GenericFunctions", which sent me further to "showMethods", which
> produced the following:
>
> showMethods("mcmcsamp")
>
> Function "mcmcsamp":
> object = "mer"
> object = "lmer"
> (inherited from object = "mer")
>
> Then I confirmed that the object "M1" created by the "lmer" call
> indeed had class "lmer". From there, I tried several things without
> success. The help("GenericFunctions") file mentioned 'dumpMethod' and
> 'dumpMethods'.
>
> > dumpMethods("mcmcsamp") # produced nothing I could find.
> > dumpMethods("mcmcsamp", file="mcmcsamp.R")
> # produced a file named "mcmcsamp.R" of official length 0 KB
> # containing nothing, as far as I could tell.
> > dumpMethods("mcmcsamp", file="mcmcsamp.R", signature="lmer")
> # also generated an empty file.
>
> > dumpMethod("mcmcsamp", "lmer")
> [1] "mcmcsamp.lmer.R"
> # Produced a file 'mcmcsamp.lmer.R' in the working directory,
> #which contained only the following:
>
> setMethod("mcmcsamp", "lmer",
> NULL
> )
>
> I also tried 'trace(mcmcsamp)' and 'trace("mcmcsamp", browser, exit =
> browser)' before running the function giving the error message. Nothing
> I saw from that seemed useful.
>
> I'd be much obliged to anyone who could help understand how I could
> diagnose this issue.

First you try

> showMethods("mcmcsamp", classes = "mer", includeDefs = TRUE)

Function "mcmcsamp":
object = "mer":
structure(function (object, n = 1, verbose = FALSE, ...) {

    .local <- function (object, n = 1, verbose = FALSE, saveb = FALSE,

        trans = TRUE, ...)
    {

        family <- object@family
        lmm <- family$family == "gaussian" && family$link ==
            "identity"
        if (!lmm)
            stop("mcmcsamp for GLMMs not yet implemented in supernodal
representation")
        ans <- t(.Call("mer_MCMCsamp", object, saveb, n, trans,
            PACKAGE = "Matrix"))
        attr(ans, "mcpar") <- as.integer(c(1, n, 1))
        class(ans) <- "mcmc"
        glmer <- FALSE
        gnms <- names(object@flist)
        cnms <- object@cnames
        ff <- fixef(object)
        colnms <- c(names(ff), if (glmer) character(0) else "sigma^2",
            unlist(lapply(seq(along = gnms), function(i) abbrvNms(gnms[i],
                cnms[[i]]))))
        if (trans) {
            ptyp <- c(integer(length(ff)), if (glmer) integer(0) else 1:1,
                unlist(lapply(seq(along = gnms), function(i) {
                  k <- length(cnms[[i]])
                  rep(1:2, c(k, (k * (k - 1))/2))
                })))
            colnms[ptyp == 1] <- paste("log(", colnms[ptyp ==
                1], ")", sep = "")
            colnms[ptyp == 2] <- paste("atanh(", colnms[ptyp ==
                2], ")", sep = "")
        }
        colnames(ans) <- colnms
        ans

    }
    .local(object, n, verbose, ...)
}, class = structure("MethodDefinition", package = "methods"), target = structure("mer", .Names = "object", class = structure("signature", package = "methods")), defined = structure("mer", .Names = "object", class = structure("signature", package = "methods")))

which tells you that the real work is being done inside a C function called mer_mcmcsamp. The sources for that function are in Matrix/src/lmer.c from the source package.

I had code for the saveb = TRUE option in there but hadn't tested it out yet. I believe that Martin has enabled it in the sources for the 0.995-3 release. That's the good news. The bad news is that we are still running tests on that version trying to find the cause of the segfault on Windows.

>
> Best Wishes,
> Spencer Graves
>
> Andrew Gelman wrote:
>
> > I am working with lmer() in the latest release of Matrix, doing various
> > things including writing a function called mcsamp() that acts as a
> > wrapper for mcmcsamp() and automatically runs multiple chains, diagnoses
> > convergence, and stores the result as a bugs object so it can be
> > plotted. I recognize that at this point, mcmcsamp() is somewhat of a
> > placeholder (since it doesn't work on a lot of models) but I'm sure it
> > will continue to be improved so I'd like to be able to work with it, as
> > a starting point if necessary.
> >
> > Anyway, I couldn't get mcmcsamp() to work with the saveb=TRUE option.
> > Here's a simple example:
> >
> > y <- 1:10
> > group <- rep (c(1,2), c(5,5))
> > M1 <- lmer (y ~ 1 + (1 | group)) # works fine
> > mcmcsamp (M1) # works fine
> > mcmcsamp (M1, saveb=TRUE)
> >
> > This last gives an error message:
> >
> > Error in "colnames<-"(`*tmp*`, value = c("(Intercept)", "log(sigma^2)", :
> > length of 'dimnames' [2] not equal to array extent
> >
> > Thanks for your help.
> > Andrew
> >
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sun Jan 22 09:36:52 2006

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:42:09 EST