[R] mcmcsamp() in lmer

Spencer Graves spencer.graves at pdf.com
Sat Jan 21 23:35:43 CET 2006


Hi, Doug:

	  Thanks for all your hard work and creativity with mixed-effects 
models -- and for your quick reply to my question about how to use 
'showMethods' to view the source.

	  Best Wishes,
	  spencer graves

Douglas Bates wrote:

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




More information about the R-help mailing list