making ran_vals work for MCMCglmm; check cleanup

R/MCMCglmm_tidiers.R

@@ -22,7 +22,7 @@
 #'   # a simple MCMCglmm model
 #'   if (require(MCMCglmm)) {
 #'      data(PlodiaPO)
-#'      m <- MCMCglmm(PO ~ 1, random = ~ FSfamily, data = PlodiaPO, verbose=FALSE)
+#'      m <- MCMCglmm(PO ~ 1, random = ~ FSfamily, data = PlodiaPO, verbose=FALSE, pr=TRUE)
 #'   }
 #'   # extract the parameter names
 #'   paramNamesMCMCglmm(m)
@@ -313,13 +313,14 @@ stdranef <- function(object, which, type = c("lp", "response"), ...) {
 ##' tidy(mm0)
 ##' tidy(mm1)
 ##' tidy(mm2)
+##' tail(tidy(mm0,effects="ran_vals"))
 ##' }
 ##' @export
 tidy.MCMCglmm <- function(x,effects=c("fixed","ran_pars"),
                           scales = NULL, ## c("sdcor","vcov",NA),
                           ...) {
     ## FIXME: allow scales= parameter to get varcov on sd/corr scale?
-    clist <- c(fixed="Sol",ran_pars="VCV",ran_vals="Liab")
+    clist <- c(fixed="Sol",ran_pars="VCV",ran_vals="Sol")
     comp <- clist[effects]
     if (!is.null(scales)) {
         if (length(scales) != length(effects)) {
@@ -332,8 +333,25 @@ tidy.MCMCglmm <- function(x,effects=c("fixed","ran_pars"),
     ## FIXME:: have to work harder to retrieve group/term information
     ##  about random parameters
     ## individual components are mcmc objects: call tidy on them
-    retList <- (purrr::map(x[comp],tidy)
+    retList <- (purrr::map(x[comp],tidy,...)
         %>% setNames(effects))
+    fnames <- paramNamesMCMCglmm(x)$fixed
+
+    if ("fixed" %in% effects) {
+        retList$fixed <- filter(retList$fixed, term %in% fnames)
+    }
+    if ("ran_vals" %in% effects) {
+        retList$ran_vals <- filter(retList$ran_vals, !(term %in% fnames))
+        if (nrow(retList$ran_vals)==0) {
+            stop("for tidying random effects values, must run MCMglmm with pr=TRUE")
+        }
+        ss <- strsplit(retList$ran_vals$term,"\\.")
+        retList$ran_vals$level <- sapply(ss,utils::tail,1)
+        retList$ran_vals$group <- sapply(ss,function(x) x[length(x)-1])
+        retList$ran_vals$term <- sapply(ss,
+                function(x) if (length(x)==3) x[1] else "(Intercept)")
+
+    }
 
     if ("ran_pars" %in% effects) {
         ss <- strsplit(retList$ran_pars$term,"(:|\\.)")

R/lme4_tidiers.R

@@ -362,8 +362,9 @@ augment.merMod <- function(x, data = stats::model.frame(x), newdata, ...) {
 #' @importFrom broom glance
 #' @export
 glance.merMod <- function(x, ...) {
+    deviance <- NULL  ## false-positive code checks
     ff <- finish_glance(x=x)
-    if (isREML(x)) ff <- rename(ff,REMLcrit=deviance)
+    if (lme4::isREML(x)) ff <- rename(ff,REMLcrit=deviance)
     return(ff)
 }
 

R/mcmc_tidiers.R

@@ -15,7 +15,8 @@
 #' @param rhat,ess (logical) include Rhat and/or effective sample size estimates?
 #' @param index Add index column, remove index from term. For example, 
 #' \code{term a[13]} becomes \code{term a} and \code{index 13}.
-#' @param ... unused
+#' @param ... mostly unused; for \code{tidy.MCMCglmm}, these represent options
+#' passed through to \code{tidy.mcmc} (e.g. \code{robust}, \code{conf.int}, \code{conf.method}, ...)
 #' 
 #' @name mcmc_tidiers
 #' @importFrom broom tidy

broom.mixed.Rproj

@@ -1,19 +1 @@
-Version: 1.0
-
-RestoreWorkspace: Default
-SaveWorkspace: Default
-AlwaysSaveHistory: Default
-
-EnableCodeIndexing: Yes
-UseSpacesForTab: Yes
-NumSpacesForTab: 2
-Encoding: UTF-8
-
-RnwWeave: Sweave
-LaTeX: pdfLaTeX
-
-BuildType: Package
-PackageUseDevtools: Yes
-PackageInstallArgs: --no-multiarch --with-keep.source
-PackageCheckArgs: --no-manual
-PackageRoxygenize: rd,collate,namespace
+Fatal error: cannot open file 'run_examples.R': No such file or directory

inst/extdata/run_examples.R

@@ -84,16 +84,17 @@ run_pkg("glmmADMB",
 run_pkg("MCMCglmm",
 {
     data("sleepstudy",package="lme4")
-    mm0 <- MCMCglmm(Reaction ~ Days, random = ~ Subject, data=sleepstudy)
+    mm0 <- MCMCglmm(Reaction ~ Days, random = ~ Subject, data=sleepstudy,
+                    pr=TRUE)
     mm1 <- MCMCglmm(Reaction ~ Days, random = ~us(1+Days):Subject,
                     ## parameter-expanded priors
                     ## V is 2x2 identity wlog
                     ## t(2) with standard dev of 2
                     prior=list(G=list(list(nu=2,V=diag(2),
                                            alpha.mu=rep(0,2),
                                            alpha.V=diag(rep(4,2))))),
-                    data=sleepstudy)
-    mm2 <- MCMCglmm(Reaction ~ Days, random = ~idh(1+Days):Subject, data=sleepstudy)
+                    data=sleepstudy, pr=TRUE)
+    mm2 <- MCMCglmm(Reaction ~ Days, random = ~idh(1+Days):Subject, data=sleepstudy, pr=TRUE)
     save_file(mm0, mm1, mm2, pkg="MCMCglmm", type = "rda")
 })