-
Notifications
You must be signed in to change notification settings - Fork 16
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* Improve rotate_factors... ... to allow working with projected values of factors * small fix to plot_factors update... ... h/t Kristen Omori * fix bug with many extrapolation-grid cells with Area=0... ... this issue arises e.g. with Region = "Other" and h/t Ellen Yasumiishi for flagging the issue! the bug manifested in an uniformative error, so no change of previous users having untrustworthy results arising from it. * fix bug in last commit * update docs * improve docs * Enhance rotate_factors ... ... to work even when estimating a factor model where some factors have variance approaching zero. * fix a plotting bug in `plot_factors` ... ... that arose when plotting factor-model Omegas * fix to `plot_factors` when doing dual ordination * Fix summary.fit_model( ..., what="residuals") ... ... in response to #72, h/t Nicholas Ducharme-Barthe for reporting the bug (which threw an uninformative error) * Adding plots for spatially-varying catchability * improve panel-labels in density plots * add loc_s and latlon_s ... ... for potential future use in new feature for habitat-specific variance-covariance among variables * fix typo in last commit * fix typo in last commit * Change default colors for plots ... ... and add option to calculate subarea totals using Longitude for Region="eastern_bering_sea" * adding units package * update docs * adding amend_output ... ... which labels Report for use in subsequent plotting * Fix bug in fit_model ... ... caused during earlier push * Improve `amend_output` * Update plot_biomass_index ... ... to use metadata_ctz, and deprecate compatibility with MIST and SpatialDeltaGLMM * fix new bug in plot_biomass_index ... ... introduced a couple commits earlier in dev branch; also improve stability for DHARMa using PIT residuals * Continue to improve stability for DHARMa residuals ... ... see discussion here: florianhartig/DHARMa#286 * update plots * Change make_covariates ... ... to avoid passing unnecessarily large X_itp object, which exceeded R memory limits in some extreme cases * small update ... ... to deal with amend_output change where Density=0 such that log(Density) = -Inf, which was causing plotting problems * Update deprecated.R functions for old diagnostic plots ... ... to deal with case of not specifying explicit units * various updates to doxygen reference docs * Update deprecated.R * Update deprecated.R * fix DHARMa residuals when b_i=NA, and improve plot_data(.) to work with projargs * fix bug in dev-changes for project_coordinates * fix bugs in dev-branch for plot_data * improve DHARMa plots ... ... by removing useless p-values * fix dev-branch bug in plot_biomass_index ... ... which now requires `extrapolation_list` to pull appropriate units for plot * adding West Coast Groundfish combo survey ... ... to VAST shapefiles * small style-guide improvements * Improve default labelling in plots ... ... by calling `amend_output` within `plot_maps` * Avoid crash in `plot_maps` ... ... when asking to plot variable that is excluded from model * Fix annoying warnings about PROJ.4 * Eliminate (most) annoying warnings from PROJ.4 * small fixes .. ... to eliminate remaining PROJ.4 warnings, and fix plotting crash regarding contours in missing variables * small fix ... ... to allow for passing a color function `col` to `plot.fit_model` * remove dependency `plotKML` with `raster` h/t @Jason-Conner-NOAA for pointing out that `plotKML` was removed from CRAN * fix bug in last commit * Simplify plot_biomass_index ... ... by incorporating `amend_output` into function * fit issues in last committ ... ... such that amend_output is now responsible for handing units * adding units for mean_Z_ctm ... ... extracted using `sf` which is now a dependency * Fix use of units for range-shift metrics * fix bug in recent commit ... ... which caused `plot_maps` to not plot SEs properly by default * update docs * fix bug in amend_outputs ... ... regarding optional Report slot `mean_D_ctl` * update docs
- Loading branch information
1 parent
80d5cbe
commit 132ae46
Showing
58 changed files
with
1,386 additions
and
693 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,161 @@ | ||
|
|
||
|
|
||
| #' Amend output from VAST for user convenience | ||
| #' | ||
| #' \code{amend_output} add labels, units, and performs logical operations, | ||
| #' e.g., adds zeros as needed, to simplify user and downstream interpretation. | ||
| #' | ||
| #' @export | ||
| amend_output <- | ||
| function( fit = NULL, | ||
| TmbData = fit$data_list, | ||
| Report = fit$Report, | ||
| extrapolation_list = fit$extrapolation_list, | ||
| Map = fit$tmb_list$Map, | ||
| Sdreport = fit$parameter_estimates$SD, | ||
| year_labels = fit$year_labels, | ||
| category_names = fit$category_names, | ||
| strata_names = fit$strata_names ){ | ||
|
|
||
| # Local functions | ||
| add_dimnames = function( Report, report_names, dimnames ){ | ||
| for( i in seq_along(report_names) ){ | ||
| if( report_names[i] %in% names(Report) ){ | ||
| dimnames(Report[[report_names[i]]]) = dimnames | ||
| } | ||
| } | ||
| return(Report) | ||
| } | ||
|
|
||
| # Defaults | ||
| if( "treat_nonencounter_as_zero" %in% names(TmbData$Options_list$Options) ){ | ||
| treat_missing_as_zero = TmbData$Options_list$Options["treat_nonencounter_as_zero"] | ||
| }else{ | ||
| treat_missing_as_zero = FALSE | ||
| } | ||
|
|
||
| # Local function | ||
| process_labels = function( labels, prefix, length ){ | ||
| if(is.null(labels)){ | ||
| labels = paste0( prefix, "_", seq_len(length) ) | ||
| }else{ | ||
| if(length(labels)!=length) stop("Check labels") | ||
| } | ||
| return(labels) | ||
| } | ||
|
|
||
| # Fill in missing inputs | ||
| if( "D_xt" %in% names(Report)){ | ||
| # SpatialDeltaGLMM | ||
| year_labels = process_labels( year_labels, "Time", ncol(Report$D_xt) ) | ||
| category_names = "singlespecies" | ||
| } | ||
| if( "D_xct" %in% names(Report)){ | ||
| # VAST Version < 2.0.0 | ||
| year_labels = process_labels( year_labels, "Time", dim(Report$D_xct)[3] ) | ||
| category_names = process_labels( category_names, "Category", dim(Report$D_xct)[2] ) | ||
| } | ||
| if( "D_xcy" %in% names(Report)){ | ||
| # VAST Version >= 2.0.0 | ||
| year_labels = process_labels( year_labels, "Time", dim(Report$D_xcy)[3] ) | ||
| category_names = process_labels( category_names, "Category_", dim(Report$D_xcy)[2] ) | ||
| } | ||
| if( "D_gcy" %in% names(Report)){ | ||
| # VAST Version 8.0.0 through 9.3.0 | ||
| year_labels = process_labels( year_labels, "Time", dim(Report$D_gcy)[3] ) | ||
| category_names = process_labels( category_names, "Category", dim(Report$D_gcy)[2] ) | ||
| } | ||
| if( "D_gct" %in% names(Report)){ | ||
| # VAST Version >= 9.4.0 | ||
| year_labels = process_labels( year_labels, "Time", dim(Report$D_gct)[3] ) | ||
| category_names = process_labels( category_names, "Category", dim(Report$D_gct)[2] ) | ||
| } | ||
| if("dhat_ktp" %in% names(Report)){ | ||
| # MIST Version <= 14 | ||
| year_labels = process_labels( year_labels, "Time", dim(Report$dhat_ktp)[2] ) | ||
| category_names = process_labels( category_names, "Category", dim(Report$dhat_ktp)[3] ) | ||
| } | ||
| if("dpred_ktp" %in% names(Report)){ | ||
| # MIST Version >= 15 | ||
| year_labels = process_labels( year_labels, "Time", dim(Report$dpred_ktp)[2] ) | ||
| category_names = process_labels( category_names, "Category", dim(Report$dpred_ktp)[3] ) | ||
| } | ||
| strata_names = process_labels( strata_names, "Stratum", dim(Report$Index_ctl)[3] ) | ||
|
|
||
| # 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) | ||
| } | ||
| } | ||
|
|
||
| # No need to map off spatial statistics mean_Z_ctm or effective_area_ctl in any years | ||
| # These are valid when betas are mapped off, or even when epsilons are mapped off given the potential covariates | ||
| #if( any(is.na(Map$beta2_ft)) | any(is.na(Map$beta2_ft)) ){ | ||
| # if("mean_Z_ctm" %in% names(Report)) Report$mean_Z_ctm = ifelse(Num_ctm==0, NA, Report$mean_Z_ctm) | ||
| # if("effective_area_ctl" %in% names(Report)) Report$effective_area_ctl = ifelse(Num_ctl==0, NA, Report$effective_area_ctl) | ||
| #} | ||
|
|
||
| # 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) ) | ||
| Report = add_dimnames( Report = Report, | ||
| report_names = c("Omega1_gc","Omega2_gc"), | ||
| dimnames = list(NULL, "Category"=category_names) ) | ||
| Report = add_dimnames( Report = Report, | ||
| report_names = "Xi1_gcp", | ||
| dimnames = list(NULL, "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)) ) | ||
| Report = add_dimnames( Report = Report, | ||
| report_names = "Phi1_gk", | ||
| dimnames = list(NULL, "Covariate"=colnames(TmbData$Q1_ik)) ) | ||
| Report = add_dimnames( Report = Report, | ||
| report_names = "Phi2_gk", | ||
| dimnames = list(NULL, "Covariate"=colnames(TmbData$Q2_ik)) ) | ||
|
|
||
| # Add labels for other useful variables | ||
| Report = add_dimnames( Report = Report, | ||
| report_names = c("Index_ctl","effective_area_ctl","mean_D_ctl"), | ||
| dimnames = list("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, colnames(TmbData$Z_gm)) ) | ||
|
|
||
| # Modify Sdreport | ||
| if( !is.null(Sdreport) ){ | ||
| # See plot_biomass_index for how to efficiently use and combine SEs | ||
| } | ||
|
|
||
| # Add units | ||
| units(Report$Index_ctl) = units(TmbData$b_i / TmbData$a_i * extrapolation_list$Area_km2[1]) | ||
| units(Report$D_gct) = units(TmbData$b_i) | ||
|
|
||
| # Add units for COG, see: https://github.com/r-quantities/units/issues/291 | ||
| # In case of re-running amend_output on objects with existing units | ||
| if( "mean_Z_ctm" %in% names(Report) ){ | ||
| units(Report$mean_Z_ctm) = as_units("km") | ||
| } | ||
| if( "mean_D_ctl" %in% names(Report) ){ | ||
| units(Report$mean_D_ctl) = units(TmbData$b_i) | ||
| } | ||
| if( "effective_area_ctl" %in% names(Report) ){ | ||
| #units(Report$effective_area_ctl) = paste0( sf::st_crs( extrapolation_list$projargs )$units, "^2" ) | ||
| units(Report$effective_area_ctl) = as_units("km^2") | ||
| } | ||
|
|
||
| # Check for bad entries | ||
| return( Report ) | ||
| } | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.