Merge pull request #26 from m-clark/extract-vc-zi-25

fix extract_vc for zi with no vc #25

R/extract_vc.R

@@ -67,11 +67,11 @@
 #' @export
 extract_vc <- function(
   model,
-  ci_level = .95,
-  ci_args = NULL,
-  ci_scale = 'sd',
-  show_cor = FALSE,
-  digits = 3,
+  ci_level  = .95,
+  ci_args   = NULL,
+  ci_scale  = 'sd',
+  show_cor  = FALSE,
+  digits    = 3,
   component = 'cond',
   ...
 ) {
@@ -92,11 +92,10 @@ extract_vc <- function(
 extract_vc.merMod <- function(
   model,
   ci_level = .95,
-  ci_args = NULL,
+  ci_args  = NULL,
   ci_scale = 'sd',
   show_cor = FALSE,
-  digits = 3,
-  # component,
+  digits   = 3,
   ...
 ) {
 
@@ -183,21 +182,26 @@ extract_vc.merMod <- function(
 #' @rdname extract_vc
 extract_vc.glmmTMB <- function(
   model,
-  ci_level = .95,
-  ci_args = NULL,
-  ci_scale = 'sd',
-  show_cor = FALSE,
-  digits = 3,
+  ci_level  = .95,
+  ci_args   = NULL,
+  ci_scale  = 'sd',
+  show_cor  = FALSE,
+  digits    = 3,
   component = 'cond',
   ...
 ) {
-  vc_mat <- glmmTMB::VarCorr(model)[[component]]
 
-  # no re allowed for dispersion formula
+  # note: disp formula doesn't allow re
   if (!component %in% c('cond', 'zi')) {
     stop('component must be one of "cond" or "zi".')
   }
 
+  vc_mat <- glmmTMB::VarCorr(model)[[component]]
+
+  if(is.null(vc_mat))
+    return(message(paste('No VarCorr for', component, 'component.')))
+
+
   # make dataframe and add names
   variance <- purrr::map(vc_mat, diag)
 
@@ -332,11 +336,10 @@ extract_vc.glmmTMB <- function(
 extract_vc.lme <- function(
   model,
   ci_level = .95,
-  ci_args = NULL,
+  ci_args  = NULL,
   ci_scale = 'sd',
   show_cor = FALSE,
-  digits = 3,
-  # component,
+  digits   = 3,
   ...
 ) {
   re_struct <- model$modelStruct$reStruct
@@ -469,11 +472,11 @@ extract_vc.lme <- function(
 #' @export
 extract_vc.brmsfit <- function(
   model,
-  ci_level = .95,
-  ci_args = NULL,
-  ci_scale = 'sd',
-  show_cor = FALSE,
-  digits = 3,
+  ci_level  = .95,
+  ci_args   = NULL,
+  ci_scale  = 'sd',
+  show_cor  = FALSE,
+  digits    = 3,
   component = NULL,
   ...
 ) {
@@ -559,11 +562,11 @@ extract_vc.brmsfit <- function(
 #' @rdname extract_vc
 extract_vc.stanreg <- function(
   model,
-  ci_level = .95,
-  ci_args = NULL,
-  ci_scale = 'sd',
-  show_cor = FALSE,
-  digits = 3,
+  ci_level  = .95,
+  ci_args   = NULL,
+  ci_scale  = 'sd',
+  show_cor  = FALSE,
+  digits    = 3,
   component = NULL,
   ...
 ) {
@@ -645,11 +648,10 @@ extract_vc.stanreg <- function(
 extract_vc.gam <- function(
   model,
   ci_level = .95,
-  ci_args = NULL,
+  ci_args  = NULL,
   ci_scale = 'sd',
   show_cor = FALSE,
-  digits = 3,
-  # component = 'cond',
+  digits   = 3,
   ...
   ) {
 

R/summarize_model.R

@@ -57,27 +57,67 @@ summarize_model <- function(
   vc <-
     extract_vc(
       model,
-      ci_level = ifelse(ci | inherits(model, 'gam'), .95, 0),
-      digits = digits,
-      show_cor = cor_re,
+      ci_level  = ifelse(ci | inherits(model, 'gam'), .95, 0),
+      digits    = digits,
+      show_cor  = cor_re,
       component = component
     )
 
-  if (cor_re == TRUE) {
-    if (inherits(model, 'gam')) {
-      cors <- 'Not estimated for gam.'
-    } else {
-      cors <- vc$Cor
-      vc   <- vc$`Variance Components`
+  if (!is.null(vc)) {
+
+    if (cor_re == TRUE) {
+      if (inherits(model, 'gam')) {
+        cors <- 'Not estimated for gam.'
+      } else {
+        cors <- vc$Cor
+        vc   <- vc$`Variance Components`
+      }
+    }
+
+    vc <- vc %>%
+      dplyr::rename_at(
+        dplyr::vars(dplyr::matches('group|effect|^var')),
+        totitle
+      ) %>%
+      dplyr::rename_at(dplyr::vars(dplyr::matches('^sd')), toupper)
+
+
+
+    ### Print re part ----
+
+    message("\nVariance Components:\n")
+
+    print(format(data.frame(vc), nsmall = digits), row.names = FALSE)
+
+    if (cor_re == TRUE) {
+      # correlations
+      message("\nCorrelation of Random Effects:\n")
+
+      if (inherits(model, 'gam'))
+        print(cors)
+      else
+        if (length(cors) == 1){
+
+          print(format(data.frame(cors[[1]]), nsmall = digits))
+
+        } else {
+          # pretty printing of multiple matrices
+          nams = names(cors)
+          nams[1] = paste0(nams[1], '\n')
+          nams[-1] = paste0('\n', nams[-1], '\n')
+
+          purrr::map2(cors, nams, function(mat, name) {
+            # cat('\n\n')
+            message(name)
+            print(format(data.frame(mat), nsmall = digits))
+          })
+        }
+
     }
   }
 
-  vc <- vc %>%
-    dplyr::rename_at(
-      dplyr::vars(dplyr::matches('group|effect|^var')),
-      totitle
-    ) %>%
-    dplyr::rename_at(dplyr::vars(dplyr::matches('^sd')), toupper)
+
+  ### Print fe part ----
 
   fe <-
     extract_fixed_effects(
@@ -92,40 +132,6 @@ summarize_model <- function(
     ) %>%
     dplyr::rename_at(dplyr::vars(dplyr::matches('se')), toupper)
 
-  ### Print re part ----
-
-  message("\nVariance Components:\n")
-
-  print(format(data.frame(vc), nsmall = digits), row.names = FALSE)
-
-  if (cor_re == TRUE) {
-    # correlations
-    message("\nCorrelation of Random Effects:\n")
-
-    if (inherits(model, 'gam'))
-      print(cors)
-    else
-      if (length(cors) == 1){
-
-        print(format(data.frame(cors[[1]]), nsmall = digits))
-
-      } else {
-        # pretty printing of multiple matrices
-        nams = names(cors)
-        nams[1] = paste0(nams[1], '\n')
-        nams[-1] = paste0('\n', nams[-1], '\n')
-
-        purrr::map2(cors, nams, function(mat, name) {
-          # cat('\n\n')
-          message(name)
-          print(format(data.frame(mat), nsmall = digits))
-        })
-      }
-
-  }
-
-  ### Print fe part ----
-
   message("\nFixed Effects:\n")
 
   print(format(data.frame(fe), nsmall = digits), row.names = FALSE)

tests/testthat/helper-load_data.R → tests/testthat/setup-load_data.R

@@ -93,6 +93,13 @@ load('tmb_cor_struct_results.RData')
 #   data = Salamanders
 # )
 #
+# tmb_zip_no_zi_re <- glmmTMB(
+#   count ~ spp + mined + (1 | site),
+#   zi =  ~ spp + mined,
+#   family = truncated_poisson,
+#   data = Salamanders
+# )
+#
 # tmb_disp <- update(tmb_2, . ~ ., dispformula =  ~Days)
 #
 # save(
@@ -102,6 +109,7 @@ load('tmb_cor_struct_results.RData')
 #   tmb_3,
 #   tmb_4,
 #   tmb_zip,
+#   tmb_zip_no_zi_re,
 #   tmb_disp,
 #   file = 'tests/testthat/tmb_results.RData'
 # )

tests/testthat/test-extract_vc.R

@@ -126,10 +126,16 @@ test_that('extract_vc.glmmTMB basic functionality: ints/slopes with multiple gro
   expect_s3_class(extract_vc(tmb_4, ci_level = 0), 'data.frame')
 })
 
+## zi
+
 test_that('extract_vc.glmmTMB basic functionality: zero-inflated', {
   expect_s3_class(extract_vc(tmb_zip, component = 'zi'), 'data.frame')
 })
 
+test_that('extract_vc.glmmTMB basic functionality: zero-inflated', {
+  expect_message(extract_vc(tmb_zip_no_zi_re, component = 'zi'))
+})
+
 ### ar and related
 
 test_that('extract_vc.glmmTMB basic functionality: ar', {