Release 0.17.2 on CRAN

Author: strengejacke

DESCRIPTION

@@ -2,8 +2,8 @@ Package: sjstats
 Type: Package
 Encoding: UTF-8
 Title: Collection of Convenient Functions for Common Statistical Computations
-Version: 0.17.1.9000
-Date: 2018-10-05
+Version: 0.17.2
+Date: 2018-11-15
 Authors@R: person("Daniel", "Lüdecke", role = c("aut", "cre"), email = "[email protected]", comment = c(ORCID = "0000-0002-8895-3206"))
 Maintainer: Daniel Lüdecke <[email protected]>
 Description: Collection of convenient functions for common statistical computations,

NAMESPACE

@@ -20,6 +20,8 @@ S3method(icc,brmsfit)
 S3method(icc,glmmTMB)
 S3method(icc,merMod)
 S3method(icc,stanreg)
+S3method(is_singular,glmmTMB)
+S3method(is_singular,merMod)
 S3method(mcse,brmsfit)
 S3method(mcse,stanmvreg)
 S3method(mcse,stanreg)
@@ -80,6 +82,7 @@ S3method(print,sj_pval)
 S3method(print,sj_r2)
 S3method(print,sj_resample)
 S3method(print,sj_revar)
+S3method(print,sj_revar_adjust)
 S3method(print,sj_rope)
 S3method(print,sj_se_icc)
 S3method(print,sj_splithalf)

R/S3-methods.R

@@ -186,7 +186,7 @@ print.sj_r2 <- function(x, digits = 3, ...) {
 #' @importFrom purrr map_chr map2_chr
 #' @export
 print.sj_icc <- function(x, digits = 4, ...) {
-  cat("\nIntra-Class Correlation Coefficient for Generalized Linear Mixed Model\n\n")
+  cat("\nIntraclass Correlation Coefficient for Generalized Linear Mixed Model\n\n")
   print_icc_and_r2(x, digits, ...)
 }
 
@@ -216,7 +216,7 @@ print_icc_and_r2 <- function(x, digits, ...) {
 #' @export
 print.sj_icc_merMod <- function(x, comp, ...) {
   # print model information
-  cat(sprintf("\n%s\n\n", attr(x, "model", exact = T)))
+  cat(sprintf("\nIntraclass Correlation Coefficient for %s\n\n", attr(x, "model", exact = T)))
 
   cat(crayon::blue(
     sprintf("Family : %s (%s)\nFormula: %s\n\n",
@@ -287,6 +287,8 @@ print.sj_icc_merMod <- function(x, comp, ...) {
                   as.vector(x[i])))
     }
   }
+
+  cat("\n")
 }
 
 
@@ -1424,6 +1426,30 @@ print.sj_grpmeans <- function(x, ...) {
 }
 
 
+#' @importFrom crayon blue cyan
+#' @export
+print.sj_revar_adjust <- function(x, ...) {
+  cat("\nVariance Components of Mixed Models\n\n")
+  cat(crayon::blue(sprintf("Family : %s (%s)\nFormula: %s\n\n", x$family, x$link, deparse(x$formula))))
+
+  vals <- c(
+    sprintf("%.3f", x$var.fixef),
+    sprintf("%.3f", x$var.ranef),
+    sprintf("%.3f", x$var.disp),
+    sprintf("%.3f", x$var.dist),
+    sprintf("%.3f", x$var.resid)
+  )
+
+  vals <- format(vals, justify = "right")
+
+  cat(sprintf("         fixed: %s\n", vals[1]))
+  cat(sprintf("        random: %s\n", vals[2]))
+  cat(sprintf("      residual: %s\n", vals[5]))
+  cat(crayon::cyan(sprintf("    dispersion: %s\n", vals[3])))
+  cat(crayon::cyan(sprintf("  distribution: %s\n\n", vals[4])))
+}
+
+
 #' @export
 print.sj_revar <- function(x, ...) {
   # get parameters

R/converge_ok.R

@@ -12,23 +12,40 @@
 #' @param tolerance Indicates up to which value the convergence result is
 #'          accepted. The smaller \code{tolerance} is, the stricter the test
 #'          will be.
+#' @param ... Currently not used.
 #'
 #' @return For \code{converge_ok()}, a logical vector, which is \code{TRUE} if
 #'           convergence is fine and \code{FALSE} if convergence is suspicious.
 #'           Additionally, the convergence value is returned as return value's name.
 #'           \code{is_singluar()} returns \code{TRUE} if the model fit is singular.
 #'
 #' @details \code{converge_ok()} provides an alternative convergence test for
-#'                \code{\link[lme4]{merMod}}-objects, as discussed
-#'                \href{https://github.com/lme4/lme4/issues/120}{here}
-#'                and suggested by Ben Bolker in
-#'                \href{https://github.com/lme4/lme4/issues/120#issuecomment-39920269}{this comment}.
-#'                \cr \cr
-#'                \code{is_singular()} checks if a model fit is singular, and can
-#'                be used in case of post-fitting convergence warnings, such as
-#'                warnings about negative eigenvalues of the Hessian. If the fit
-#'                is singular (i.e. \code{is_singular()} returns \code{TRUE}), these
-#'                warnings can most likely be ignored.
+#'   \code{\link[lme4]{merMod}}-objects, as discussed
+#'   \href{https://github.com/lme4/lme4/issues/120}{here}
+#'   and suggested by Ben Bolker in
+#'   \href{https://github.com/lme4/lme4/issues/120#issuecomment-39920269}{this comment}.
+#'   \cr \cr
+#'   If a model is "singular", this means that some dimensions of the variance-covariance
+#'   matrix have been estimated as exactly zero. \code{is_singular()} checks if
+#'   a model fit is singular, and can be used in case of post-fitting convergence
+#'   warnings, such as warnings about negative eigenvalues of the Hessian. If the fit
+#'   is singular (i.e. \code{is_singular()} returns \code{TRUE}), these warnings
+#'   can most likely be ignored.
+#'   \cr \cr
+#'   There is no gold-standard about how to deal with singularity and which
+#'   random-effects specification to choose. Beside using fully Bayesian methods
+#'   (with informative priors), proposals in a frequentist framework are:
+#'   \itemize{
+#'   \item avoid fitting overly complex models, such that the variance-covariance matrices can be estimated precisely enough (\cite{Matuschek et al. 2017})
+#'   \item use some form of model selection to choose a model that balances predictive accuracy and overfitting/type I error (\cite{Bates et al. 2015}, \cite{Matuschek et al. 2017})
+#'   \item \dQuote{keep it maximal}, i.e. fit the most complex model consistent with the experimental design, removing only terms required to allow a non-singular fit (\cite{Barr et al. 2013})
+#'   }
+#'
+#' @references \itemize{
+#'   \item Bates D, Kliegl R, Vasishth S, Baayen H. Parsimonious Mixed Models. arXiv:1506.04967, June 2015.
+#'   \item Barr DJ, Levy R, Scheepers C, Tily HJ. Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3):255–278, April 2013.
+#'   \item Matuschek H, Kliegl R, Vasishth S, Baayen H, Bates D. Balancing type I error and power in linear mixed models. Journal of Memory and Language, 94:305–315, 2017.
+#'   }
 #'
 #' @examples
 #' library(sjmisc)
@@ -74,23 +91,25 @@ converge_ok <- function(x, tolerance = 0.001) {
 }
 
 
-#' @importFrom lme4 getME
-#' @importFrom glmmTMB getME
 #' @rdname converge_ok
 #' @export
-is_singular <- function(x, tolerance = 1e-5) {
-  if (is_merMod(x)) {
-    theta <- lme4::getME(x, "theta")
-    # diagonal elements are identifiable because they are fitted
-    #  with a lower bound of zero ...
-    diag.element <- lme4::getME(x, "lower") == 0
-    any(abs(theta[diag.element]) < tolerance)
-  } else if (inherits(x, "glmmTMB")) {
-    theta <- glmmTMB::getME(x, "theta")
-    # diagonal elements are identifiable because they are fitted
-    #  with a lower bound of zero ...
-    diag.element <- glmmTMB::getME(x, "lower") == 0
-    any(abs(theta[diag.element]) < tolerance)
-  } else
-    warning("`x` must be a merMod- or glmmTMB-object.", call. = F)
+is_singular <- function(x, tolerance = 1e-5, ...) {
+  UseMethod("is_singular")
+}
+
+#' @importFrom lme4 getME
+#' @export
+is_singular.merMod <- function(x, tolerance = 1e-5, ...) {
+  theta <- lme4::getME(x, "theta")
+  # diagonal elements are identifiable because they are fitted
+  #  with a lower bound of zero ...
+  diag.element <- lme4::getME(x, "lower") == 0
+  any(abs(theta[diag.element]) < tolerance)
+}
+
+#' @importFrom lme4 VarCorr
+#' @export
+is_singular.glmmTMB <- function(x, tolerance = 1e-5, ...) {
+  vc <- collapse_cond(lme4::VarCorr(x))
+  any(sapply(vc, function(.x) any(abs(diag(.x)) < tolerance)))
 }

R/icc.R

@@ -19,11 +19,9 @@
 #' @param ppd Logical, if \code{TRUE}, variance decomposition is based on the
 #'   posterior predictive distribution, which is the correct way for Bayesian
 #'   non-Gaussian models.
-#' @param adjusted Logical, if \code{TRUE}, the adjusted (and conditional) ICC
-#'   is calculated, which reflects the uncertainty of all random effects (see
-#'   'Details'). \strong{Note} that if \code{adjusted = TRUE}, \strong{no}
-#'   additional information on the variance components is returned.
-#'
+#' @param adjusted Logical, if \code{TRUE}, the adjusted (and
+#'   conditional) ICC is calculated, which reflects the uncertainty of all
+#'   random effects (see 'Details').
 #'
 #' @inheritParams hdi
 #'
@@ -103,7 +101,7 @@
 #'    To get a meaningful ICC also for models with random slopes, use \code{adjusted = TRUE}.
 #'    The adjusted ICC used the mean random effect variance, which is based
 #'    on the random effect variances for each value of the random slope
-#'    (see \cite{Johnson 2014}).
+#'    (see \cite{Johnson et al. 2014}).
 #'    \cr \cr
 #'    \strong{ICC for models with multiple or nested random effects}
 #'    \cr \cr
@@ -759,15 +757,23 @@ icc.brmsfit <- function(x, re.form = NULL, typical = "mean", prob = .89, ppd = F
 #'                objects are supported.
 #'
 #' @param x Fitted mixed effects model (of class \code{merMod}, \code{glmmTMB},
-#'          \code{stanreg} or \code{brmsfit}). \code{get_re_var()} also accepts
-#'           an object of class \code{icc.lme4}, as returned by the
-#'           \code{\link{icc}} function.
+#'   \code{stanreg} or \code{brmsfit}). \code{get_re_var()} also accepts
+#'   an object of class \code{icc.lme4}, as returned by the
+#'   \code{\link{icc}} function.
 #' @param comp Name of the variance component to be returned. See 'Details'.
+#' @param adjusted Logical, if \code{TRUE}, returns the variance of the fixed
+#'   and random effects as well as of the additive dispersion and
+#'   distribution-specific variance, which are used to calculate the
+#'   adjusted and conditional \code{\link{r2}} and \code{\link{icc}}.
 #'
 #' @return \code{get_re_var()} returns the value of the requested variance component,
 #'           \code{re_var()} returns all random effects variances.
 #'
-#' @references Aguinis H, Gottfredson RK, Culpepper SA. 2013. Best-Practice Recommendations for Estimating Cross-Level Interaction Effects Using Multilevel Modeling. Journal of Management 39(6): 1490–1528 (\doi{10.1177/0149206313478188})
+#' @references \itemize{
+#'    \item Aguinis H, Gottfredson RK, Culpepper SA. 2013. Best-Practice Recommendations for Estimating Cross-Level Interaction Effects Using Multilevel Modeling. Journal of Management 39(6): 1490–1528 (\doi{10.1177/0149206313478188})
+#'    \item Johnson PC, O'Hara RB. 2014. Extension of Nakagawa & Schielzeth's R2GLMM to random slopes models. Methods Ecol Evol, 5: 944-946. (\doi{10.1111/2041-210X.12225})
+#'    \item Nakagawa S, Johnson P, Schielzeth H (2017) The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisted and expanded. J. R. Soc. Interface 14. \doi{10.1098/rsif.2017.0213}
+#'  }
 #'
 #' @details The random effect variances indicate the between- and within-group
 #'         variances as well as random-slope variance and random-slope-intercept
@@ -785,6 +791,18 @@ icc.brmsfit <- function(x, re.form = NULL, typical = "mean", prob = .89, ppd = F
 #'         direct effects) affect the between-group-variance. Cross-level
 #'         interaction effects are group-level factors that explain the
 #'         variance in random slopes (Aguinis et al. 2013).
+#'         \cr \cr
+#'         If \code{adjusted = TRUE}, the variance of the fixed and random
+#'         effects as well as of the additive dispersion and
+#'         distribution-specific variance are returned (see \cite{Johnson et al. 2014}
+#'         and \cite{Nakagawa et al. 2017}):
+#'         \describe{
+#'          \item{\code{"fixed"}}{variance attributable to the fixed effects}
+#'          \item{\code{"random"}}{variance of random effects}
+#'          \item{\code{"dispersion"}}{variance due to additive dispersion}
+#'          \item{\code{"distribution"}}{distribution-specific variance}
+#'          \item{\code{"residual"}}{sum of dispersion and distribution}
+#'         }
 #'
 #' @seealso \code{\link{icc}}
 #'
@@ -794,6 +812,7 @@ icc.brmsfit <- function(x, re.form = NULL, typical = "mean", prob = .89, ppd = F
 #'
 #' # all random effect variance components
 #' re_var(fit1)
+#' re_var(fit1, adjusted = TRUE)
 #'
 #' # just the rand. slope-intercept covariance
 #' get_re_var(fit1, "tau.01")
@@ -806,20 +825,40 @@ icc.brmsfit <- function(x, re.form = NULL, typical = "mean", prob = .89, ppd = F
 #' @importFrom purrr map map2 flatten_dbl flatten_chr
 #' @importFrom sjmisc trim
 #' @export
-re_var <- function(x) {
-  # iterate all attributes and return them as vector
-  rv <- c("sigma_2", "tau.00", "tau.11", "tau.01", "rho.01")
+re_var <- function(x, adjusted = FALSE) {
 
-  # compute icc
-  icc_ <- suppressMessages(icc(x))
+  if (adjusted) {
 
-  rv_ <- purrr::map(rv, ~ attr(icc_, .x, exact = TRUE))
-  rn <- purrr::map2(1:length(rv_), rv, ~ sjmisc::trim(paste(names(rv_[[.x]]), .y, sep = "_")))
-  rv_ <- purrr::flatten_dbl(rv_)
+    rv <- r2(x)
 
-  names(rv_) <- purrr::flatten_chr(rn)[1:length(rv_)]
+    rv_ <- list(
+      var.fixef = attr(rv, "var.fixef", exact = TRUE),
+      var.ranef = attr(rv, "var.ranef", exact = TRUE),
+      var.disp = attr(rv, "var.disp", exact = TRUE),
+      var.dist = attr(rv, "var.dist", exact = TRUE),
+      var.resid = attr(rv, "var.resid", exact = TRUE),
+      formula = attr(rv, "formula", exact = TRUE),
+      family = attr(rv, "family", exact = TRUE),
+      link = attr(rv, "link", exact = TRUE)
+    )
+
+    class(rv_) <- c("sj_revar_adjust", class(rv_))
+
+  } else {
+    # iterate all attributes and return them as vector
+    rv <- c("sigma_2", "tau.00", "tau.11", "tau.01", "rho.01")
 
-  class(rv_) <- c("sj_revar", class(rv_))
+    # compute icc
+    icc_ <- suppressMessages(icc(x))
+
+    rv_ <- purrr::map(rv, ~ attr(icc_, .x, exact = TRUE))
+    rn <- purrr::map2(1:length(rv_), rv, ~ sjmisc::trim(paste(names(rv_[[.x]]), .y, sep = "_")))
+    rv_ <- purrr::flatten_dbl(rv_)
+
+    names(rv_) <- purrr::flatten_chr(rn)[1:length(rv_)]
+
+    class(rv_) <- c("sj_revar", class(rv_))
+  }
 
   rv_
 }

R/pred_vars.R

@@ -3,7 +3,7 @@
 #'
 #' @description Several functions to retrieve information from model objects,
 #'    like variable names, link-inverse function, model frame,
-#'    model_family etc., in a tidy and consistent way.
+#'    model family etc., in a tidy and consistent way.
 #'
 #' @param x A fitted model; for \code{var_names()}, \code{x} may also be a
 #'    character vector.