Fix printing/tidying with random intercepts + time-varying #426

Since new ranef implementation

Author: seananderson

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: sdmTMB
 Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
-Version: 0.6.0.9033
+Version: 0.6.0.9034
 Authors@R: c(
     person(c("Sean", "C."), "Anderson", , "sean@seananderson.ca",
            role = c("aut", "cre"),

R/methods.R

@@ -186,12 +186,12 @@ ranef.sdmTMB <- function(object, ...) {
     for (j in 1:length(groups)) {
       sub <- .t[which(.t$group_name == groups[j]), ]
       level_ids <- unique(sub$level_ids)
-      sub <- sub[, c("group_name", "par_name", "estimate")]
+      sub <- sub[, c("group_name", "term", "estimate")]
       if (nrow(sub) > 0) {
         # convert long to wide, storing just estimates
-        split_data <- split(sub$estimate, sub$par_name)
+        split_data <- split(sub$estimate, sub$term)
         wide_df <- as.data.frame(split_data) # Convert to wide format
-        names(wide_df) <- unique(sub$par_name) # rename, fix .X issue
+        names(wide_df) <- unique(sub$term) # rename, fix .X issue
         rownames(wide_df) <- level_ids # add rownames, like lmer does
         # Create a list with the dataframe as an element named 'Dog'
         group_list[[j]] <- wide_df

R/tidy.R

@@ -305,12 +305,13 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals", "ran_vco
   }
 
   # optional time-varying random components
-  if(!is.null(x$time_varying)) {
+  if (!is.null(x$time_varying)) {
     tv_names <- colnames(model.matrix(x$time_varying, x$data))
     time_slices <- x$time_lu$time_from_data
     yrs <- rep(time_slices, times = length(tv_names))
 
     out_ranef_tv <- data.frame(
+      model = model,
       term = paste0(rep(tv_names, each = length(time_slices)), ":", yrs),
       estimate = c(est$b_rw_t),
       std.error = c(se$b_rw_t),
@@ -322,6 +323,8 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals", "ran_vco
     if(is.null(out_ranef)) {
       out_ranef <- out_ranef_tv
     } else {
+      out_ranef_tv$group_name <- NA
+      out_ranef_tv$level_ids <- NA
       out_ranef <- rbind(out_ranef, out_ranef_tv)
     }
   }
@@ -395,23 +398,23 @@ get_re_tidy_list <- function(x, crit) {
   re_b_df$estimate <- x$sd_report$value[re_indx][non_nas]
   re_b_df$std.error <- x$sd_report$sd[re_indx][non_nas]
   re_b_df$conf.low <- re_b_df$estimate - crit * re_b_df$std.error
-  re_b_df$conf.hi <- re_b_df$estimate + crit * re_b_df$std.error
+  re_b_df$conf.high <- re_b_df$estimate + crit * re_b_df$std.error
   re_b_df$index <- NULL
   re_b_df$group_name <- NA
-  re_b_df$par_name <- NA
+  re_b_df$term <- NA
   for (i in seq_len(length(x$split_formula))) {
     groupnames <- names(x$split_formula[[i]]$re_cov_terms$cnms)
     for (j in seq_len(length(x$split_formula[[i]]$barnames))) {
       model_grp <- which(re_b_df$model == i & re_b_df$group_id == j)
       re_b_df$group_name[model_grp] <- groupnames[j]
-      re_b_df$par_name[model_grp] <- rep(x$split_formula[[i]]$re_cov_terms$cnms[[j]], length.out = length(model_grp))
+      re_b_df$term[model_grp] <- rep(x$split_formula[[i]]$re_cov_terms$cnms[[j]], length.out = length(model_grp))
     }
   }
   group_key <- stats::aggregate(group_name ~ model + group_id, data = re_b_df, FUN = function(x) x[1])
 
   # more sensible re-ordering
   re_b_df$group_id <- NULL
-  re_b_df <- re_b_df[, c("model", "group_name", "par_name", "level_ids", "estimate", "std.error", "conf.low", "conf.hi")]
+  re_b_df <- re_b_df[, c("model", "group_name", "term", "level_ids", "estimate", "std.error", "conf.low", "conf.high")]
   # remove ":" in the level_ids
   re_b_df$level_ids <- sapply(strsplit(re_b_df$level_ids, ":"), function(x) x[2])
   out_ranef <- re_b_df
@@ -439,21 +442,21 @@ get_re_tidy_list <- function(x, crit) {
   re_cov_df$estimate <- x$sd_report$value[re_indx][non_nas]
   re_cov_df$std.error <- x$sd_report$sd[re_indx][non_nas]
   re_cov_df$conf.low <- re_cov_df$estimate - crit * re_cov_df$std.error
-  re_cov_df$conf.hi <- re_cov_df$estimate + crit * re_cov_df$std.error
+  re_cov_df$conf.high <- re_cov_df$estimate + crit * re_cov_df$std.error
   # the SD parameters are returned in log space -- use delta method to generate CIs
   est <- exp(re_cov_df$estimate) # estimate in normal space
   sd_est <- est * re_cov_df$std.error # SE in normal space
   sds <- which(re_cov_df$is_sd == 1) # index which elements are SDs
   re_cov_df$estimate[sds] <- est[sds]
   re_cov_df$conf.low[sds] <- est[sds] - crit * sd_est[sds]
-  re_cov_df$conf.hi[sds] <- est[sds] + crit * sd_est[sds]
+  re_cov_df$conf.high[sds] <- est[sds] + crit * sd_est[sds]
 
-  re_cov_df <- re_cov_df[, c("rows", "cols", "model", "group", "estimate", "std.error", "conf.low", "conf.hi")]
+  re_cov_df <- re_cov_df[, c("rows", "cols", "model", "group", "estimate", "std.error", "conf.low", "conf.high")]
   cov_matrices_lo = create_cov_matrices(re_cov_df, col_name = "conf.low")
-  cov_matrices_hi = create_cov_matrices(re_cov_df, col_name = "conf.hi")
+  cov_matrices_hi = create_cov_matrices(re_cov_df, col_name = "conf.high")
   list(out_ranef = out_ranef, cov_matrices = create_cov_matrices(re_cov_df),
        cov_matrices_lo = create_cov_matrices(re_cov_df, col_name = "conf.low"),
-       cov_matrices_hi = create_cov_matrices(re_cov_df, col_name = "conf.hi"))
+       cov_matrices_hi = create_cov_matrices(re_cov_df, col_name = "conf.high"))
 }
 
 create_cov_matrices <- function(df, col_name = "estimate") {

tests/testthat/test-tidy.R

@@ -86,3 +86,33 @@ test_that("tidy works", {
   expect_true("cv_split" %in% names(allmodels))
   expect_equal(allmodels$cv_split, sort(rep(1:2,2)))
 })
+
+test_that("printing/tidying works with a delta model that has random intercepts + an AR1 time series #426", {
+  skip_on_cran()
+  d <- pcod
+  d$fake <- rep(c("a", "b", "c"), 9999)[1:nrow(d)]
+  fit <- sdmTMB(
+    density ~ breakpt(depth_scaled) + (1|fake),
+    data = d,
+    time = "year",
+    extra_time = c(2006, 2008, 2010, 2012, 2014, 2016),
+    time_varying = ~1,
+    time_varying_type = "ar1",
+    mesh = mesh,
+    spatial = "off",
+    spatiotemporal = "off",
+    family = delta_gamma(type = "poisson-link")
+  )
+  b <- tidy(fit, effects = "ran_pars")
+  b <- tidy(fit, effects = "ran_vals")
+  expect_identical(b$term, c("(Intercept)", "(Intercept)", "(Intercept)", "(Intercept)",
+    "(Intercept)", "(Intercept)", "(Intercept):2003", "(Intercept):2004",
+    "(Intercept):2005", "(Intercept):2006", "(Intercept):2007", "(Intercept):2008",
+    "(Intercept):2009", "(Intercept):2010", "(Intercept):2011", "(Intercept):2012",
+    "(Intercept):2013", "(Intercept):2014", "(Intercept):2015", "(Intercept):2016",
+    "(Intercept):2017", "(Intercept):2003", "(Intercept):2004", "(Intercept):2005",
+    "(Intercept):2006", "(Intercept):2007", "(Intercept):2008", "(Intercept):2009",
+    "(Intercept):2010", "(Intercept):2011", "(Intercept):2012", "(Intercept):2013",
+    "(Intercept):2014", "(Intercept):2015", "(Intercept):2016", "(Intercept):2017"
+  ))
+})