Re: [R] mcmcsamp() in lmer

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Sun 22 Jan 2006 - 09:45:03 EST

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

    DougB> On 1/21/06, Spencer Graves <spencer.graves@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@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@flist) cnms <- object@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@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@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

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 10:13:03 2006

This archive was generated by hypermail 2.1.8 : Sun 22 Jan 2006 - 14:10:45 EST