[R] mcmcsamp() in lmer

Martin Maechler maechler at stat.math.ethz.ch
Sat Jan 21 23:45:03 CET 2006 

>>>>> "DougB" == Douglas Bates <dmbates at gmail.com>
>>>>>     on Sat, 21 Jan 2006 16:13:37 -0600 writes:

    DougB> On 1/21/06, Spencer Graves <spencer.graves at pdf.com>
    DougB> 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.

    DougB> First you try

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

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

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

    DougB> I had code for the saveb = TRUE option in there but
    DougB> hadn't tested it out yet.  I believe that Martin has
    DougB> enabled it in the sources for the 0.995-3 release.

yes, indeed, that's been fixed (and tests/lmer.R now also has
a test case).

I hadn't got around to reply to Andrew's e-mail,..
Martin

    DougB> That's the good news.  The bad news is that we are
    DougB> still running tests on that version trying to find
    DougB> 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 at 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
    >> 

    DougB> ______________________________________________
    DougB> R-help at stat.math.ethz.ch mailing list
    DougB> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE
    DougB> do read the posting guide!
    DougB> http://www.R-project.org/posting-guide.html

More information about the R-help mailing list