Merge pull request #407 from pbs-assess/issue_374

issue 374

Author: seananderson

DESCRIPTION

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

NEWS.md

@@ -67,6 +67,8 @@
 * Add option for `area` to be passed in as the name of a column in the 
   dataframe to be used for area weighting. Used in `get_index()`, 
   `get_cog()`, `get_eao()`, etc.
+  
+* Add time varying AR(1) correlation to tidy() and print() #374
 
 # sdmTMB 0.6.0
 

R/print.R

@@ -181,6 +181,15 @@ print_time_varying <- function(x, m = 1) {
     colnames(mm_tv) <- c("coef.est", "coef.se")
     time_slices <- x$time_lu$time_from_data
     row.names(mm_tv) <- paste(rep(tv_names, each = length(time_slices)), time_slices, sep = "-")
+
+    # append AR(1) parameters if they are estimated
+    p <- tidy(x, effects = "ran_pars", model = m, silent = TRUE)
+    if(any(p$term == "rho_time")) {
+      rho_tv <- as.matrix(p[p$term=="rho_time",c("estimate","std.error")])
+      colnames(rho_tv) <- colnames(mm_tv)
+      row.names(rho_tv) <- paste("rho", tv_names, sep = "-")
+      mm_tv <- rbind(mm_tv, rho_tv)
+    }
   } else {
     mm_tv <- NULL
   }

R/tidy.R

@@ -107,8 +107,12 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model =
     p$sigma_O <- as.numeric(p$sigma_O[1,model])
     p$sigma_Z <- as.numeric(p$sigma_Z[,model])
     p$sigma_G <- as.numeric(p$sigma_G[,model])
+
+    # if delta, a single AR1 -> rho_time_unscaled is a 1x2 matrix
+    p$rho_time <- 2 * plogis(p$rho_time_unscaled[,model]) - 1
     p
   }
+
   est <- subset_pars(est, model)
   se <- subset_pars(se, model)
 
@@ -186,6 +190,10 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model =
     log_name <- c(log_name, "ln_tau_G")
     name <- c(name, "sigma_G")
   }
+  if (!all(est$rho_time==0)) {
+    log_name <- c(log_name, "rho_time_unscaled")
+    name <- c(name, "rho_time")
+  }
 
   j <- 0
   if (!"log_range" %in% names(est)) {
@@ -204,6 +212,8 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model =
       this <- non_log_name[j]
       if (this == "tau_G") this <- "sigma_G"
       if (this == "tau_V") this <- "sigma_V"
+      if (this == "rho_time_unscaled") this <- "rho_time"
+
       this_se <- as.numeric(se[[this]])
       this_est <- as.numeric(est[[this]])
       if (length(this_est) && !(all(this_se == 0) && all(this_est == 0))) {
@@ -213,6 +223,19 @@ tidy.sdmTMB <- function(x, effects = c("fixed", "ran_pars", "ran_vals"), model =
           conf.high = exp(.e + crit * .se),
           stringsAsFactors = FALSE
         )
+        if(this == "rho_time") {
+          out_re[[i]] <- data.frame(
+            term = i,
+            estimate = this_est,
+            # use delta method to get SE in normal space
+            std.error = 2 * plogis (.e[,model]) * (1 - plogis (.e[,model])) * .se[,model],
+            # don't use delta-method for CIs, because they can be outside (-1,1)
+            conf.low = 2 * plogis(.e[,model] - crit * .se[,model]) - 1,
+            conf.high = 2 * plogis(.e[,model] + crit * .se[,model]) - 1,
+            stringsAsFactors = FALSE
+          )
+
+        }
       }
       ii <- ii + 1
     }

tests/testthat/test-time-varying.R

@@ -68,6 +68,79 @@ test_that("AR1 time-varying works", {
   expect_equal(dim(ss), c(4000L, 1L))
   expect_gt(exp(s$ln_tau_V)[1,1], 0.75)
   expect_gt(m121(s$rho_time_unscaled)[1,1], 0.65)
+
+  # test that tidy works
+  fit <- sdmTMB(density ~ 1, time = "year",
+                time_varying = ~ 1,
+                time_varying_type = "ar1",
+                data = pcod_2011, spatial="off",
+                spatiotemporal = "off",
+                family = tweedie())
+  s <- tidy(fit, "ran_pars")
+  expect_equal(dim(s), c(3L, 5L))
+  expect_equal(s$term, c("phi","rho_time","tweedie_p"))
+
+  # test that tidy works -- RW
+  fit <- sdmTMB(density ~ as.factor(year), time = "year",
+                time_varying = ~ -1 + depth_scaled,
+                data = pcod_2011, spatial="off",
+                spatiotemporal = "off",
+                family = tweedie())
+  s <- tidy(fit, "ran_pars")
+  expect_equal(dim(s), c(2L, 5L))
+  expect_equal(s$term, c("phi","tweedie_p"))
+
+  # test that tidy works -- AR1
+  fit <- sdmTMB(density ~ as.factor(year), time = "year",
+                time_varying = ~ -1 + depth_scaled,
+                data = pcod_2011, spatial="off",
+                time_varying_type = "ar1",
+                spatiotemporal = "off",
+                family = tweedie())
+  s <- tidy(fit, "ran_pars")
+  expect_equal(dim(s), c(3L, 5L))
+  expect_equal(s$term, c("phi", "rho_time", "tweedie_p"))
+
+  # test that tidy works -- no time varying
+  fit <- sdmTMB(density ~ as.factor(year), time = "year",
+                data = pcod_2011, spatial="off",
+                spatiotemporal = "off",
+                family = tweedie())
+  s <- tidy(fit, "ran_pars")
+  expect_equal(dim(s), c(2L, 5L))
+  expect_equal(s$term, c("phi","tweedie_p"))
+
+  # test that tidy works with 2 parameters
+  fit <- sdmTMB(density ~ 1, time = "year",
+                time_varying = ~ -1+depth_scaled + I(depth_scaled^2),
+                time_varying_type = "ar1",
+                data = pcod_2011, spatial="off",
+                spatiotemporal = "off",
+                family = tweedie())
+  s <- tidy(fit, "ran_pars")
+  expect_equal(dim(s), c(4L, 5L))
+  expect_equal(s$term, c("phi","rho_time","rho_time","tweedie_p"))
+
+  # test that tidy works with delta
+  fit <- sdmTMB(density ~ 1, time = "year",
+                time_varying = ~ -1+depth_scaled + I(depth_scaled^2),
+                time_varying_type = "ar1",
+                data = pcod_2011, spatial="off",
+                spatiotemporal = "off",
+                family = delta_gamma())
+  s <- tidy(fit, "ran_pars")
+  expect_equal(dim(s), c(2L, 5L))
+  expect_equal(s$term, c("rho_time","rho_time"))
+
+  s <- tidy(fit, "ran_pars", model = 1)
+  expect_equal(dim(s), c(2L, 5L))
+  expect_equal(s$term, c("rho_time","rho_time"))
+  expect_equal(s$estimate, c(1.0, 0.9582931))
+
+  s <- tidy(fit, "ran_pars", model = 2)
+  expect_equal(s$term, c("phi", "rho_time","rho_time"))
+  expect_equal(s$estimate, c(0.652, 1.0, 0.868), tolerance = 0.001)
+
 })