small updates and fixed some bugs particularly in stdranef

JWiley · Apr 18, 2013 · 95ad495 · 95ad495
1 parent b05a292
commit 95ad495
Show file tree

Hide file tree

Showing 5 changed files with 26 additions and 11 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: postMCMCglmm
 Type: Package
-Version: 0.1-1
-Date: 2013-1-30
+Version: 0.1-2
+Date: 2013-4-17
 Title: Average marginal predicted probabilities from Bayesian ordered
     probit models
 Author: Joshua Wiley <jwiley.psych@gmail.com>

diff --git a/R/extraction.R b/R/extraction.R
@@ -200,7 +200,7 @@ ranef.MCMCglmm <- function(object, use = c("all", "mean"), ...) {
 #'   deviation on the linear predictor metric, \sQuote{lp} or
 #'   response, \sQuote{response}.
 #' @param \dots Not currently used.
-#' @return A list of class postMCMCglmmRE with individual estimates and means
+#' @return A list of class postMCMCglmmRE with means (\code{M}) and individual estimates (\code{Data})
 #' @export
 #' @seealso \code{\link{print.postMCMCglmmRE}}, \code{\link{predict2.MCMCglmm}}, \code{\link{ranef.MCMCglmm}}
 #' @examples
@@ -238,14 +238,19 @@ ranef.MCMCglmm <- function(object, use = c("all", "mean"), ...) {
 #'
 #'   stdranef(m, which = list("FSfamily", "plate"), type = "lp")
 #'
-#'   # this does not work, check zero setting
-#'   #stdranef(m, type = "response")
+#'   # mean standard deviations on the probability metric
+#'   # also the full distributions, if desired in the Data slot.
+#'   res <- stdranef(m, type = "response")
+#'   res$M # means
+#'   hist(res$Data$FSfamily[, 1]) # histogram
 #' }
 stdranef <- function(object, which, type = c("lp", "response"), ...) {
   type <- match.arg(type)
 
   if (is.null(object$Z)) stop("Z matrix must be saved")
-  z <- object$Z
+  ## z <- object$Z
+  z <- diag(ncol(object$Z))
+  colnames(z) <- colnames(object$Z)
 
   re <- paramNamesMCMCglmm(object)$random
 
@@ -277,7 +282,9 @@ stdranef <- function(object, which, type = c("lp", "response"), ...) {
     }, response = {
       yhat <- lapply(index, function(i) {
         tmp <- z
-        tmp[, -i] <- 0L # bug here??
+        if (length(i) < ncol(tmp)) {
+          tmp[, -i] <- 0L # zero out all random effects we are not interested in
+        }
         predict2(object, X = NULL, Z = tmp, use = "all", type = type)
       })
       lapply(yhat, function(m) sapply(m, function(n) apply(n, 1, sd)))

diff --git a/demo/postMCMCglmm.R b/demo/postMCMCglmm.R
@@ -93,3 +93,8 @@ finaltable
 
 ## write to clipboard on Windows (e.g., for collaborators who want in Excel)
 #write.table(finaltable, file = "clipboard", na = "", sep = "\t", row.names=FALSE, col.names=TRUE)
+
+## now look at the average standard deviation of the distribution of predicted
+## probabilities including the random effects, given a fixed effect
+## probability of .5 see ?stdranef for more details
+stdranef(m, type = "response")
diff --git a/inst/doc/postMCMCglmm-manual.pdf b/inst/doc/postMCMCglmm-manual.pdf
diff --git a/man/stdranef.Rd b/man/stdranef.Rd
@@ -19,8 +19,8 @@
   \item{\dots}{Not currently used.}
 }
 \value{
-  A list of class postMCMCglmmRE with individual estimates
-  and means
+  A list of class postMCMCglmmRE with means (\code{M}) and
+  individual estimates (\code{Data})
 }
 \description{
   Function designed to extract the standard deviation of
@@ -70,8 +70,11 @@
 
   stdranef(m, which = list("FSfamily", "plate"), type = "lp")
 
-  # this does not work, check zero setting
-  #stdranef(m, type = "response")
+  # mean standard deviations on the probability metric
+  # also the full distributions, if desired in the Data slot.
+  res <- stdranef(m, type = "response")
+  res$M # means
+  hist(res$Data$FSfamily[, 1]) # histogram
 }
 }
 \seealso{