Merging 2.12.0

* adding `plot_clusters` * removing deprecated functions from NAMESPACE * small fixes * fix bug in `plot_similarity` ... ... where it would throw an error when using a univariate model * improve plot_similarity ... ... to use smarter color-limits than the current default in corrplot package * update docs * fix plotting bugs in `plot_clusters` and `plot_similarity` * update docs for `plot_maps` * fix bug in last update to `plot_range_index` * fix `plot_clusters` to work well with NAs * add informative error for `plot_range_edge` ... ... for cases when normal approximation used by `sample_variable` performs poorly * prune dependencies using `renv::dependencies` * Try adding INLA to IMPORTS using `Additional_repositories` * moving SpatialDeltaGLMM fits to FishStatsUtils ... ... so that I can then remove SpatialDeltaGLMM from dependencies of VAST * noting `summarize_covariance` is deprecated * housekeeping: remove .onAttach action, and deprecate map_hypervariance * eliminate dependency `gap` * small fixes * adding Lindegren et al. ref to `plot_clusters` * fix obj$report() when using intern=TRUE * improved tracking for run times * update messages * fix bug * Deprecate Calc_Anisotropic_Mesh * a * fix bug in plot_similarity ... ... where plot.hclust(.) throws an error with only 2 groups * Update version number * adding `form` to plot.fit_model ... .. so that residuals can be easily plotted by category * improve `plot_clusters` and `summary.fit_model` ... ... in support of h/t Joe Bizzarro and West Coast diet analysis * update docs * fix `plot_similarity` bug ... ... which threw uniformative error when doing EpsilonTime for a model with only one time-interval * better fix for plot_similarity .. ... when n_t = 1 * add warning for incorrect class in `plot_variable` * adding `plot_residual_semivariance` ... ... as additional diagnostic for model fit * small updates * fix typo in `plot_residual_semivariance` * Update fit_model ... ... to allow user to pass data_list, spatial_list, or extrapolation_list and thereby override default functions. See #83 * fix bug in `plot_residual_semivariance` and ... ... add backdoor for new maps for barrier feature, see #84 * fix bug in last commit * Improve `plot_results` ... ... by allowing `Expansion_cz` to be passed to `calculation_proportion` * adding sf option for plotting maps * adding new epsilon-correction by default * small plotting fix for format="sf" ... ... h/t Sarah Gaichas for flagging the issue * adding demo code for SVC paper * fix to work with suplus-production example to use Bratio and Fratio * increase eval.max and iter.max defaults * various doc updates
James-Thorson-NOAA · Nov 18, 2022 · 57e18bb · 57e18bb
1 parent c8bf4df
commit 57e18bb
Show file tree

Hide file tree

Showing 77 changed files with 1,612 additions and 942 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: FishStatsUtils
 Type: Package
 Title: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox
-Version: 2.11.0
-Date: 2022-03-17
+Version: 2.12.0
+Date: 2022-11-18
 Authors@R: 
     c(person(given = "James",
              family = "Thorson",
@@ -14,47 +14,46 @@ Description: FishStatsUtils contains utilities (shared code and data) used by mu
   packages (VAST, SpatialDeltaGLMM, MIST, Spatial_FA, SpatialDFA, surplus_production, EOFR) that are designed
   for spatio-temporal analysis of ecological data.  
 Imports:
-    graphics,
-    utils,
-    mapproj,
-    maptools,
-    deldir,
+    abind,
+    alphahull,
+    corpcor,
+    corrplot,
+    DHARMa,
+    ecodist,
+    fastcluster,
+    INLA,
+    plotrix,
     RANN,
-    stats,
-    colorspace,
-    shape,
-    devtools,
-    mixtools,
-    sp,
     sf,
     raster,
-    plotrix,
-    TMB,
-    MatrixModels,
+    reshape2,
     rgdal,
-    ThorsonUtilities,
-    TMBhelper,
-    abind,
-    corpcor,
-    pander,
     rnaturalearth,
     rnaturalearthdata,
-    formatR,
-    splancs,
-    DHARMa,
-    viridisLite
+    seriation,
+    shape,
+    sp,
+    stats,
+    TMB (>= 1.8.0),
+    TMBhelper (>= 1.4.0),
+    utils,
+    viridisLite,
+    spacetime,
+    gstat
 Depends:
     R (>= 3.5.0),
     units
+Enhances:
+    tidyr
 Suggests:
     testthat
 Remotes: 
-    james-thorson/utilities,
     kaskr/TMB_contrib_R/TMBhelper
+Additional_repositories: https://inla.r-inla-download.org/R/stable
 License: GPL-3
 LazyData: yes
 BuildVignettes: yes
 Encoding: UTF-8
-RoxygenNote: 7.1.2
+RoxygenNote: 7.2.2
 URL: http://github.com/james-thorson-NOAA/FishStatsUtils
 BugReports: http://github.com/james-thorson-NOAA/FishStatsUtils/issues
diff --git a/NAMESPACE b/NAMESPACE
@@ -9,11 +9,7 @@ S3method(print,make_extrapolation_info)
 S3method(print,make_spatial_info)
 S3method(summary,fit_model)
 export(Calc_Anisotropic_Mesh)
-export(Calc_Kmeans)
 export(Calc_Polygon_Areas_and_Polygons_Fn)
-export(Geostat_Sim)
-export(PlotLoadings)
-export(PlotMap_Fn)
 export(Prepare_AI_Extrapolation_Data_Fn)
 export(Prepare_BC_Coast_Extrapolation_Data_Fn)
 export(Prepare_BFISH_MHI_Extrapolation_Data_Fn)
@@ -34,16 +30,15 @@ export(Prepare_User_Extrapolation_Data_Fn)
 export(Prepare_WCGBTS_Extrapolation_Data_Fn)
 export(Prepare_WCGHL_Extrapolation_Data_Fn)
 export(Print_Message)
-export(Rotate_Fn)
 export(amend_output)
 export(calc_cov)
+export(calculate_knot_areas)
 export(calculate_proportion)
 export(combine_extrapolation_info)
 export(combine_lists)
 export(convert_shapefile)
 export(convert_version_name)
 export(fit_model)
-export(format_covariates)
 export(get_latest_version)
 export(load_example)
 export(make_covariates)
@@ -53,31 +48,27 @@ export(make_map_info)
 export(make_mesh)
 export(make_settings)
 export(make_spatial_info)
-export(map_hypervariance)
 export(match_strata_fn)
 export(plot_anisotropy)
 export(plot_biomass_index)
-export(plot_cov)
+export(plot_clusters)
 export(plot_data)
 export(plot_dharma)
-export(plot_encounter_diagnostic)
 export(plot_factors)
 export(plot_index)
-export(plot_lines)
 export(plot_loadings)
 export(plot_maps)
 export(plot_overdispersion)
-export(plot_quantile_diagnostic)
 export(plot_quantile_residuals)
 export(plot_range_edge)
 export(plot_range_index)
-export(plot_residuals)
+export(plot_residual_semivariance)
 export(plot_results)
+export(plot_similarity)
 export(plot_timeseries)
 export(plot_variable)
 export(project_coordinates)
 export(rotate_factors)
 export(sample_variable)
 export(simulate_data)
 export(strip_units)
-export(summarize_covariance)
diff --git a/R/Calc_Anisotropic_Mesh.R b/R/Calc_Anisotropic_Mesh.R
diff --git a/R/amend_output.R b/R/amend_output.R
@@ -85,18 +85,24 @@ function( fit = NULL,
   }
 
   # Determine year-category pairs with no data
-  Num_gct = rep(1,TmbData$n_g) %o% abind::adrop(TmbData$Options_list$metadata_ctz[,,'num_notna',drop=FALSE], drop=3)
-  Num_ctl = abind::adrop(TmbData$Options_list$metadata_ctz[,,'num_notna',drop=FALSE], drop=3) %o% rep(1,TmbData$n_l)
-  Num_ctm = abind::adrop(TmbData$Options_list$metadata_ctz[,,'num_notna',drop=FALSE], drop=3) %o% rep(1,TmbData$n_m)
-  if( treat_missing_as_zero==TRUE ){
-    # if treat_missing_as_zero==TRUE, then switch density from year-categories with no data to zero
-    Report$D_gct = ifelse(Num_gct==0, 0, Report$D_gct)
-    Report$Index_ctl = ifelse(Num_ctl==0, 0, Report$Index_ctl)
-  }else{
-    # If some intercepts are mapped off, then switch density from year-categories with no data to NA
-    if( any(is.na(Map$beta2_ft)) | any(is.na(Map$beta2_ft)) ){
-      Report$D_gct = ifelse(Num_gct==0, NA, Report$D_gct)
-      Report$Index_ctl = ifelse(Num_ctl==0, NA, Report$Index_ctl)
+  if( "metadata_ctz" %in% names(TmbData$Options_list) ){
+    Num_gct = rep(1,TmbData$n_g) %o% abind::adrop(TmbData$Options_list$metadata_ctz[,,'num_notna',drop=FALSE], drop=3)
+    Num_ctl = abind::adrop(TmbData$Options_list$metadata_ctz[,,'num_notna',drop=FALSE], drop=3) %o% rep(1,TmbData$n_l)
+    Num_ctm = abind::adrop(TmbData$Options_list$metadata_ctz[,,'num_notna',drop=FALSE], drop=3) %o% rep(1,TmbData$n_m)
+    if( treat_missing_as_zero==TRUE ){
+      # if treat_missing_as_zero==TRUE, then switch density from year-categories with no data to zero
+      if("D_gct"%in%names(Report)) Report$D_gct = ifelse(Num_gct==0, 0, Report$D_gct)
+      if("D_gcy"%in%names(Report)) Report$D_gcy = ifelse(Num_gct==0, 0, Report$D_gcy)
+      if("Index_ctl"%in%names(Report)) Report$Index_ctl = ifelse(Num_ctl==0, 0, Report$Index_ctl)
+      if("Index_cyl"%in%names(Report)) Report$Index_cyl = ifelse(Num_ctl==0, 0, Report$Index_cyl)
+    }else{
+      # If some intercepts are mapped off, then switch density from year-categories with no data to NA
+      if( any(is.na(Map$beta2_ft)) | any(is.na(Map$beta2_ft)) ){
+        if("D_gct"%in%names(Report)) Report$D_gct = ifelse(Num_gct==0, NA, Report$D_gct)
+        if("D_gcy"%in%names(Report)) Report$D_gcy = ifelse(Num_gct==0, NA, Report$D_gcy)
+        if("Index_ctl"%in%names(Report)) Report$Index_ctl = ifelse(Num_ctl==0, NA, Report$Index_ctl)
+        if("Index_cyl"%in%names(Report)) Report$Index_cyl = ifelse(Num_ctl==0, NA, Report$Index_cyl)
+      }
     }
   }
 
@@ -110,39 +116,42 @@ function( fit = NULL,
   # Add labels for all variables plotted using `plot_maps`
   Report = add_dimnames( Report = Report,
                          report_names = c("P1_gct","P2_gct","R1_gct","R2_gct","D_gct","Epsilon1_gct","Epsilon2_gct","eta1_gct","eta2_gct"),
-                         dimnames = list(NULL, "Category"=category_names, "Time"=year_labels) )
+                         dimnames = list("Site"=seq_len(TmbData$n_g), "Category"=category_names, "Time"=year_labels) )
   Report = add_dimnames( Report = Report,
                          report_names = c("Omega1_gc","Omega2_gc"),
-                         dimnames = list(NULL, "Category"=category_names) )
+                         dimnames = list("Site"=seq_len(TmbData$n_g), "Category"=category_names) )
   Report = add_dimnames( Report = Report,
                          report_names = "Xi1_gcp",
-                         dimnames = list(NULL, "Category"=category_names, "Covariate"=colnames(TmbData$X1_ip)) )
+                         dimnames = list("Site"=seq_len(TmbData$n_g), "Category"=category_names, "Covariate"=colnames(TmbData$X1_ip)) )
   Report = add_dimnames( Report = Report,
                          report_names = "Xi2_gcp",
-                         dimnames = list(NULL, "Category"=category_names, "Covariate"=colnames(TmbData$X2_ip)) )
+                         dimnames = list("Site"=seq_len(TmbData$n_g), "Category"=category_names, "Covariate"=colnames(TmbData$X2_ip)) )
   Report = add_dimnames( Report = Report,
                          report_names = "Phi1_gk",
-                         dimnames = list(NULL, "Covariate"=colnames(TmbData$Q1_ik)) )
+                         dimnames = list("Site"=seq_len(TmbData$n_g), "Covariate"=colnames(TmbData$Q1_ik)) )
   Report = add_dimnames( Report = Report,
                          report_names = "Phi2_gk",
-                         dimnames = list(NULL, "Covariate"=colnames(TmbData$Q2_ik)) )
+                         dimnames = list("Site"=seq_len(TmbData$n_g), "Covariate"=colnames(TmbData$Q2_ik)) )
   Report = add_dimnames( Report = Report,
                          report_names = c("L_omega1_cf","L_omega2_cf","L_beta1_cf","L_beta2_cf","L_epsilon1_cf","L_epsilon2_cf"),
                          dimnames = list("Category"=category_names, NULL) )
   Report = add_dimnames( Report = Report,
-                         report_names = c("Ltime_epsilon2_tf","Ltime_epsilon2_tf"),
+                         report_names = c("Ltime_epsilon1_tf","Ltime_epsilon2_tf"),
                          dimnames = list("Time"=year_labels, NULL) )
   Report = add_dimnames( Report = Report,
                          report_names = c("beta1_tc","beta2_tc"),
                          dimnames = list("Time"=year_labels, "Category"=category_names) )
 
   # Add labels for other useful variables
   Report = add_dimnames( Report = Report,
-                         report_names = c("Index_ctl","effective_area_ctl","mean_D_ctl"),
+                         report_names = c("Index_ctl","effective_area_ctl","mean_D_ctl","Bratio_ctl"),
                          dimnames = list("Category"=category_names, "Time"=year_labels, "Stratum"=strata_names) )
+  Report = add_dimnames( Report = Report,
+                         report_names = c("Fratio_ct"),
+                         dimnames = list("Category"=category_names, "Time"=year_labels) )
   Report = add_dimnames( Report = Report,
                          report_names = c("Index_gctl"),
-                         dimnames = list(NULL, "Category"=category_names, "Time"=year_labels, "Stratum"=strata_names) )
+                         dimnames = list("Site"=seq_len(TmbData$n_g), "Category"=category_names, "Time"=year_labels, "Stratum"=strata_names) )
   Report = add_dimnames( Report = Report,
                          report_names = "mean_Z_ctm",
                          dimnames = list("Category"=category_names, "Time"=year_labels, "Spatial_axis"=colnames(TmbData$Z_gm)) )
@@ -154,6 +163,8 @@ function( fit = NULL,
 
   # Add units
   if("Index_ctl" %in% names(Report)) units(Report$Index_ctl) = units(TmbData$b_i / TmbData$a_i * extrapolation_list$Area_km2[1])
+  if("Fratio_ct" %in% names(Report)) units(Report$Fratio_ct) = units(unitless)
+  if("Bratio_ctl" %in% names(Report)) units(Report$Bratio_ctl) = units(unitless)
   if("D_gct" %in% names(Report)) units(Report$D_gct) = units(TmbData$b_i / TmbData$a_i)
 
   # Add units for COG, see: https://github.com/r-quantities/units/issues/291

diff --git a/R/Calc_Polygon_Areas_and_Polygons_Fn.R → R/calculate_knot_areas.R b/R/Calc_Polygon_Areas_and_Polygons_Fn.R → R/calculate_knot_areas.R
@@ -1,9 +1,9 @@
 
 #' Calculate areas and covariates for a given triangulated mesh
 #'
-#' \code{Calc_Polygon_Areas_and_Polygons_Fn} builds outputs for a given triangulated mesh used for approximating spatial variation
+#' \code{calculate_knot_areas} builds outputs for a given triangulated mesh used for approximating spatial variation
 #'
-#' @inheritParams Calc_Anisotropic_Mesh
+#' @inheritParams make_mesh
 
 #' @param Data_Extrap the output from e.g., \code{Extrapolation_List <- Prepare_WCGBTS_Extrapolation_Data_Fn()} using slot \code{Extrapolation_List$Data_Extrap}
 #' @param Covariates character vector giving names of columns from Data_Extrap that should be used as covariates (default="none", which results in a design matrix with one columns of 1s)
@@ -15,8 +15,11 @@
 #' }
 
 #' @export
-Calc_Polygon_Areas_and_Polygons_Fn <-
-function( loc_x, Data_Extrap, Covariates="none", a_el=NULL ){
+calculate_knot_areas <-
+function( loc_x,
+          Data_Extrap,
+          Covariates = "none",
+          a_el = NULL ){
 
   # Calculate a_el if missing
   if( is.null(a_el) ){