Add support for group-level estimates from coxme (#1087)

* Add support for group-level estimates from coxme

* fix

* refactor

* use insight

* fix

* fix

* as.data.frame method

* ...

* fixes

* ...

* fix

* add test

* seed

* add test

* fix

* fix

* update news

* wordlist

* fix test

* fix

* Update DESCRIPTION

* fix

* comment

* fix

* test

* fix warning

* fix

Author: strengejacke

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: parameters
 Title: Processing of Model Parameters
-Version: 0.27.0.3
+Version: 0.27.0.5
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(as.data.frame,VarCorr.coxme)
 S3method(as.data.frame,VarCorr.lme)
 S3method(as.double,n_clusters)
 S3method(as.double,n_factors)
@@ -225,6 +226,7 @@ S3method(model_parameters,clmm2)
 S3method(model_parameters,coeftest)
 S3method(model_parameters,compare.loo)
 S3method(model_parameters,comparisons)
+S3method(model_parameters,coxme)
 S3method(model_parameters,cpglmm)
 S3method(model_parameters,data.frame)
 S3method(model_parameters,dbscan)

NEWS.md

@@ -16,6 +16,9 @@
   generic `display()` method can be used, too, which has a `format` argument that
   also supports `"tt"` for `tinytable`.
 
+* Improved support for *coxme* models in `model_parameters()`. Random effects
+  and group level estimates are now returned as well.
+
 ## Bug fixes
 
 * Fixed issue with models of class `selection` with multiple outcomes.

R/extract_random_parameters.R

@@ -122,3 +122,57 @@
 .extract_random_parameters.MixMod <- function(model, ...) {
   NULL
 }
+
+
+.extract_random_parameters.coxme <- function(model, ci = NULL, effects = "random", ...) {
+  insight::check_if_installed(c("lme4", "coxme"))
+
+  # extract random effects
+  re_grp <- insight::find_random(model, split_nested = TRUE, flatten = TRUE)
+
+  do.call(rbind, lapply(re_grp, function(grp) {
+    # coerce to data frame
+    out <- as.data.frame(lme4::ranef(model)[[grp]], stringsAsFactors = FALSE)
+    colnames(out)[1] <- "(Intercept)"
+
+    # reshape, to have a long format
+    out$.grp <- unique(insight::get_data(model)[[grp]])
+    out <- datawizard::reshape_longer(out, select = -".grp")
+    out$grpvar <- grp
+
+    # rename columns
+    colnames(out) <- c("Level", "Parameter", "Coefficient", "Group")
+    out$SE <- NA
+
+    # re-order columns
+    out <- out[c("Group", "Parameter", "Level", "Coefficient", "SE")]
+    out$Parameter[out$Parameter == "Intercept"] <- "(Intercept)"
+
+    # sort
+    out <- datawizard::data_arrange(out, c("Group", "Parameter", "Level"))
+
+    # coerce to character
+    out$Parameter <- as.character(out$Parameter)
+    out$Level <- as.character(out$Level)
+    out$Group <- as.character(out$Group)
+    out$Effects <- "random"
+
+    out$CI_low <- out$CI_high <- NA
+    ci_cols <- c("CI_low", "CI_high")
+
+    stat_column <- gsub("-statistic", "", insight::find_statistic(model), fixed = TRUE)
+
+    # to match rbind
+    out[[stat_column]] <- NA
+    out$df_error <- NA
+    out$p <- NA
+
+    out <- out[c("Parameter", "Level", "Coefficient", "SE", ci_cols, stat_column, "df_error", "p", "Effects", "Group")]
+
+    if (effects == "random") {
+      out[c(stat_column, "df_error", "p")] <- NULL
+    }
+
+    out
+  }))
+}

R/extract_random_variances.R

@@ -159,8 +159,18 @@
                                              ci_random = NULL,
                                              verbose = FALSE,
                                              ...) {
-  varcorr <- .get_variance_information(model, component)
-  if (!inherits(model, "lme")) {
+  # special handling for lme objects
+  if (inherits(model, "lme")) {
+    insight::check_if_installed("lme4")
+    varcorr <- lme4::VarCorr(model)
+    class(varcorr) <- "VarCorr.lme"
+  } else {
+    varcorr <- insight::get_mixed_info(model, component = component, verbose = FALSE)$vc
+  }
+
+  # we have own data frame methods for VarCorr objects from lme and coxme, so
+  # only change class attribute for other models
+  if (!inherits(model, c("lme", "coxme"))) {
     class(varcorr) <- "VarCorr.merMod"
   }
 
@@ -332,6 +342,47 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ...
 }
 
 
+#' @export
+as.data.frame.VarCorr.coxme <- function(x, row.names = NULL, optional = FALSE, ...) {
+  # extract variances from VarCorr object
+  variances <- lapply(x, diag)
+  # create data frame, similar to as.data.frame.VarCorr.merMod
+  out <- do.call(rbind, lapply(names(variances), function(i) {
+    # information on variances
+    d <- data.frame(
+      grp = i,
+      var1 = names(variances[[i]]),
+      var2 = NA_character_,
+      vcov = as.numeric(variances[[i]]),
+      sdcor = sqrt(as.numeric(variances[[i]])),
+      stringsAsFactors = FALSE
+    )
+    # add correlations, if any
+    if (nrow(x[[i]]) > 1) {
+      d <- rbind(d, data.frame(
+        grp = i,
+        var1 = "(Intercept)",
+        var2 = rownames(x[[i]])[2],
+        vcov = NA_real_,
+        sdcor = as.numeric(x[[i]][2, 1]),
+        stringsAsFactors = FALSE
+      ))
+    }
+    d
+  }))
+
+  # bind residual variance
+  rbind(out, data.frame(
+    grp = "Residual",
+    var1 = NA_character_,
+    var2 = NA_character_,
+    vcov = NA_real_,
+    sdcor = NA_real_,
+    stringsAsFactors = FALSE
+  ))
+}
+
+
 # extract CI for random SD ------------------------
 
 .random_sd_ci <- function(model,
@@ -684,188 +735,6 @@ as.data.frame.VarCorr.lme <- function(x, row.names = NULL, optional = FALSE, ...
 }
 
 
-# Extract Variance and Correlation Components ----
-
-# store essential information about variance components...
-# basically, this function should return lme4::VarCorr(x)
-.get_variance_information <- function(model, model_component = "conditional") {
-  # check if packages are available
-  .check_mixedmodels_namespace(model)
-
-  # stanreg
-  # ---------------------------
-  if (inherits(model, "stanreg")) {
-    varcorr <- lme4::VarCorr(model)
-
-    # GLMMapdative
-    # ---------------------------
-  } else if (inherits(model, "MixMod")) {
-    vc1 <- vc2 <- NULL
-    re_names <- insight::find_random(model)
-
-    vc_cond <- !startsWith(colnames(model$D), "zi_")
-    if (any(vc_cond)) {
-      vc1 <- model$D[vc_cond, vc_cond, drop = FALSE]
-      attr(vc1, "stddev") <- sqrt(diag(vc1))
-      attr(vc1, "correlation") <- stats::cov2cor(model$D[vc_cond, vc_cond, drop = FALSE])
-    }
-
-    vc_zi <- startsWith(colnames(model$D), "zi_")
-    if (any(vc_zi)) {
-      colnames(model$D) <- gsub("^zi_(.*)", "\\1", colnames(model$D))
-      rownames(model$D) <- colnames(model$D)
-      vc2 <- model$D[vc_zi, vc_zi, drop = FALSE]
-      attr(vc2, "stddev") <- sqrt(diag(vc2))
-      attr(vc2, "correlation") <- stats::cov2cor(model$D[vc_zi, vc_zi, drop = FALSE])
-    }
-
-    vc1 <- list(vc1)
-    names(vc1) <- re_names[[1]]
-    attr(vc1, "sc") <- sqrt(insight::get_deviance(model, verbose = FALSE) / insight::get_df(model, type = "residual", verbose = FALSE)) # nolint
-    attr(vc1, "useSc") <- TRUE
-
-    if (!is.null(vc2)) {
-      vc2 <- list(vc2)
-      names(vc2) <- re_names[[2]]
-      attr(vc2, "sc") <- sqrt(insight::get_deviance(model, verbose = FALSE) / insight::get_df(model, type = "residual", verbose = FALSE)) # nolint
-      attr(vc2, "useSc") <- FALSE
-    }
-
-    varcorr <- insight::compact_list(list(vc1, vc2))
-    names(varcorr) <- c("cond", "zi")[seq_along(varcorr)]
-
-    # joineRML
-    # ---------------------------
-  } else if (inherits(model, "mjoint")) {
-    re_names <- insight::find_random(model, flatten = TRUE)
-    varcorr <- summary(model)$D
-    attr(varcorr, "stddev") <- sqrt(diag(varcorr))
-    attr(varcorr, "correlation") <- stats::cov2cor(varcorr)
-    varcorr <- list(varcorr)
-    names(varcorr) <- re_names[1]
-    attr(varcorr, "sc") <- model$coef$sigma2[[1]]
-    attr(varcorr, "useSc") <- TRUE
-
-    # nlme
-    # ---------------------------
-  } else if (inherits(model, "lme")) {
-    varcorr <- lme4::VarCorr(model)
-
-    # ordinal
-    # ---------------------------
-  } else if (inherits(model, "clmm")) {
-    varcorr <- ordinal::VarCorr(model)
-    attr(varcorr, "useSc") <- FALSE
-
-    # glmmadmb
-    # ---------------------------
-  } else if (inherits(model, "glmmadmb")) {
-    varcorr <- lme4::VarCorr(model)
-
-    # brms
-    # ---------------------------
-  } else if (inherits(model, "brmsfit")) {
-    varcorr <- lapply(names(lme4::VarCorr(model)), function(i) {
-      element <- lme4::VarCorr(model)[[i]]
-      if (i != "residual__") {
-        if (is.null(element$cov)) {
-          out <- as.matrix(drop(element$sd[, 1])^2)
-          colnames(out) <- rownames(out) <- gsub("Intercept", "(Intercept)", rownames(element$sd), fixed = TRUE)
-        } else {
-          out <- as.matrix(drop(element$cov[, 1, ]))
-          colnames(out) <- rownames(out) <- gsub("Intercept", "(Intercept)", rownames(element$cov), fixed = TRUE)
-        }
-        attr(out, "sttdev") <- element$sd[, 1]
-      } else {
-        out <- NULL
-      }
-      out
-    })
-    varcorr <- insight::compact_list(varcorr)
-    names(varcorr) <- setdiff(names(lme4::VarCorr(model)), "residual__")
-    attr(varcorr, "sc") <- lme4::VarCorr(model)$residual__$sd[1, 1]
-
-    # cpglmm
-    # ---------------------------
-  } else if (inherits(model, "cpglmm")) {
-    varcorr <- cplm::VarCorr(model)
-    attr(varcorr, "useSc") <- FALSE
-
-    # lme4 / glmmTMB
-    # ---------------------------
-  } else {
-    varcorr <- lme4::VarCorr(model)
-  }
-
-
-  # for glmmTMB, tell user that dispersion model is ignored
-
-  if (inherits(model, c("glmmTMB", "MixMod"))) {
-    if (is.null(model_component) || model_component == "conditional") {
-      varcorr <- .collapse_cond(varcorr)
-    } else {
-      varcorr <- .collapse_zi(varcorr)
-    }
-  }
-
-  varcorr
-}
-
-
-.check_mixedmodels_namespace <- function(model) {
-  # reason to be installed
-  reason <- "to compute random effect variances for mixed models"
-
-  # installed?
-  insight::check_if_installed("lme4", reason = reason)
-
-  if (inherits(model, "lme")) {
-    insight::check_if_installed("nlme", reason = reason)
-  }
-
-  if (inherits(model, "glmmTMB")) {
-    insight::check_if_installed("glmmTMB", reason = reason)
-  }
-
-  if (inherits(model, "clmm")) {
-    insight::check_if_installed("ordinal", reason = reason)
-  }
-
-  if (inherits(model, "brmsfit")) {
-    insight::check_if_installed("brms", reason = reason)
-  }
-
-  if (inherits(model, "cpglmm")) {
-    insight::check_if_installed("cplm", reason = reason)
-  }
-
-  if (inherits(model, "rstanarm")) {
-    insight::check_if_installed("rstanarm", reason = reason)
-  }
-}
-
-
-# glmmTMB returns a list of model information, one for conditional
-# and one for zero-inflation part, so here we "unlist" it, returning
-# only the conditional part.
-.collapse_cond <- function(x) {
-  if (is.list(x) && "cond" %in% names(x)) {
-    x[["cond"]]
-  } else {
-    x
-  }
-}
-
-
-.collapse_zi <- function(x) {
-  if (is.list(x) && "zi" %in% names(x)) {
-    x[["zi"]]
-  } else {
-    x
-  }
-}
-
-
 # this is used to only temporarily load merDeriv and to point registered
 # methods from merDeriv to lme4-methods. if merDeriv was loaded before,
 # nothing will be changed. If merDeriv was not loaded, vcov-methods registered