#' Plot results
#'
#' \code{plot_results} plots diagnostics, results, and indices for a given fitted model
#'
#' This function takes a fitted VAST model and generates a standard set of diagnostic and visualization plots.
#' It does this by calling a series of mid-level plotting functions; see list of functions in Value section of documentation.
#'
#' In particular, for making customized maps of output please see \code{\link{plot_variable}}
#'
#' @inheritParams fit_model
#' @inheritParams plot_maps
#' @inheritParams plot_variable
#' @inheritParams plot_quantile_residuals
#' @inheritParams plot_range_edge
#' @inheritParams simulate_data
#' @inheritParams plot_factors
#' @inheritParams plot_similarity
#' @inheritParams plot_clusters
#'
#' @param fit Output from \code{fit_model}
#' @param check_residuals Boolean indicating whether to run or skip residual diagnostic plots (which can be slow as currently implemented)
#' @param cluster_results Boolean whether to run \code{\link{plot_clusters}}, which can be slow for analyses involving a large number
#' of extrapolation-grid cells
#' @param ... additional settings to pass to \code{\link{plot_maps}}
#'
#' @return Invisibly returns a tagged list of outputs generated by standard plots. See linked functions for details
#' \describe{
#' \item{\code{dharmaRes}}{Output from \code{\link{summary.fit_model}}, representing quantile residuals calculated using package DHARMa}
#' \item{\code{Enc_prob}}{Output from \code{\link{plot_encounter_diagnostic}}}
#' \item{\code{Index}}{Output from \code{\link{plot_biomass_index}}}
#' \item{\code{Proportions}}{Output from \code{\link{calculate_proportion}}}
#' \item{\code{Range}}{Output from \code{\link{plot_range_index}}}
#' \item{\code{Dens_xt}}{Output from \code{\link{plot_maps}}}
#' \item{\code{Edge}}{Output from \code{\link{plot_range_edge}}}
#' \item{\code{Factors}}{Output from \code{\link{plot_factors}}}
#' \item{\code{Clusters}}{Output from \code{\link{plot_clusters}}}
#' }
#'
#' @family wrapper functions
#' @seealso \code{\link[VAST]{VAST}} for general documentation, \code{\link{make_settings}} for generic settings, \code{\link{fit_model}} for model fitting, and \code{\link{plot_results}} for generic plots
#'
#' @export
plot_results <-
function( fit,
settings = fit$settings,
plot_set = 3,
working_dir = getwd(),
year_labels = fit$year_labels,
years_to_plot = fit$years_to_plot,
category_names = fit$category_names,
strata_names = fit$strata_names,
use_biascorr = TRUE,
map_list = NULL,
check_residuals = TRUE,
cluster_results = TRUE,
#projargs = fit$extrapolation_list$projargs,
projargs = '+proj=longlat',
zrange,
n_samples = 100,
calculate_relative_to_average = FALSE,
type = 1,
n_cells = NULL,
n_cells_residuals = NULL,
RotationMethod = "PCA",
quantiles = c(0.05,0.5,0.95),
similarity_metric = c("hclust", "Correlation", "Dissimilarity", "Covariance")[1],
... ){
# Check for known issues
if( is.null(fit$Report)) stop("`fit$Report` is missing, please check inputs")
if( is.null(category_names)) category_names = paste0( "Category_", 1:fit$data_list$n_c )
# Make directory
dir.create(working_dir, showWarnings=FALSE, recursive=TRUE)
message("\n### Creating plots in directory ", working_dir )
# plot data
message("\n### Making plots of data availability and knots")
plot_data_args = list(...)
plot_data_args = combine_lists( "input"=plot_data_args, "args_to_use"=formalArgs(plot_data),
"default" = list( Extrapolation_List = fit$extrapolation_list,
Spatial_List = fit$spatial_list,
Lat_i = fit$data_frame[,'Lat_i'],
Lon_i = fit$data_frame[,'Lon_i'],
Year_i = fit$data_frame[,'t_i'],
PlotDir = working_dir,
year_labels = year_labels,
projargs = projargs) )
do.call( what=plot_data, args=plot_data_args )
# PLot settings
if( is.null(map_list) ){
message("\n### Obtaining default settings for plotting maps")
map_list = make_map_info( Region = settings$Region,
spatial_list = fit$spatial_list,
Extrapolation_List = fit$extrapolation_list )
}
# Plot anisotropy
message("\n### Making plot of anisotropy")
plot_anisotropy( FileName = file.path(working_dir,"Aniso.png"),
Obj = fit$tmb_list$Obj )
# Plot index
plot_biomass_index_args = list(...)
if( !is.null(fit$parameter_estimates$SD) ){
message("\n### Making plot of abundance index")
#if( !all(is.numeric(year_labels)) ) stop("`plot_biomass_index` isn't built to handle non-numeric `year_labels`")
plot_biomass_index_args = combine_lists( "input"=plot_biomass_index_args, "args_to_use"=formalArgs(plot_biomass_index),
"default" = list(DirName = working_dir,
fit = fit,
year_labels = year_labels,
years_to_plot = years_to_plot,
use_biascorr = use_biascorr,
category_names = category_names,
strata_names = strata_names) )
Index = do.call( what=plot_biomass_index, args=plot_biomass_index_args )
}else{
Index = "Not run"
message("\n### Skipping plot of abundance index; must re-run with standard errors to plot")
}
# Plot covariance/dissimilarity matrices
plot_similarity_args = list(...)
message("\n### Making plot of covariance/dissimilarity matrices")
#if( !all(is.numeric(year_labels)) ) stop("`plot_biomass_index` isn't built to handle non-numeric `year_labels`")
plot_similarity_args = combine_lists( "input"=plot_similarity_args, "args_to_use"=formalArgs(plot_similarity),
"default" = list(fit = fit,
year_labels = year_labels,
category_names = category_names,
similarity_metric = similarity_metric,
working_dir = working_dir) )
do.call( what=plot_similarity, args=plot_similarity_args )
# Plot comps
if( !is.null(fit$parameter_estimates$SD) & fit$data_list$n_c>1 ){
message("\n### Making plot of composition data")
#if( !all(is.numeric(year_labels)) ) stop("`plot_biomass_index` isn't built to handle non-numeric `year_labels`")
calculate_proportion_args = list(...)
calculate_proportion_args = combine_lists( "input"=calculate_proportion_args, "args_to_use"=formalArgs(calculate_proportion),
"default" = list( fit = fit,
TmbData = fit$data_list,
Index = Index,
year_labels = year_labels,
years_to_plot = years_to_plot,
use_biascorr = use_biascorr,
category_names = category_names,
DirName = working_dir,
n_samples = n_samples ) )
Proportions = do.call( what=calculate_proportion, args=calculate_proportion_args )
#Compositions = plot_biomass_index( DirName=working_dir, TmbData=fit$data_list, Sdreport=fit$parameter_estimates$SD, year_labels=year_labels,
# years_to_plot=years_to_plot, use_biascorr=use_biascorr, category_names=category_names )
}else{
Proportions = "Not run"
message("\n### Skipping plot of composition data; must re-run with standard errors and multiple categories to plot")
}
# Plot range indices
if( !is.null(fit$parameter_estimates$SD) ){
message("\n### Making plot of spatial indices")
#if( !all(is.numeric(year_labels)) ) stop("`plot_range_index` isn't built to handle non-numeric `year_labels`")
Range = plot_range_index( Report = fit$Report,
TmbData = fit$data_list,
Sdreport = fit$parameter_estimates$SD,
Znames = colnames(fit$data_list$Z_xm),
PlotDir = working_dir,
year_labels = year_labels,
years_to_plot = years_to_plot,
use_biascorr = use_biascorr,
category_names = category_names )
}else{
Range = "Not run"
message("\n### Skipping plot of spatial indices; must re-run with standard errors to plot")
}
# Plot range edges
if( "jointPrecision"%in%names(fit$parameter_estimates$SD) & n_samples>0 & fit$data_list$Options_list$Options['Calculate_Range']==TRUE ){
message("\n### Making plot of range edges")
Edge = plot_range_edge( Obj = fit$tmb_list$Obj,
Sdreport = fit$parameter_estimates$SD,
working_dir = working_dir,
year_labels = year_labels,
years_to_plot = years_to_plot,
category_names = category_names,
n_samples = n_samples,
quantiles = quantiles,
calculate_relative_to_average = calculate_relative_to_average )
}else{
Edge = "Not run"
message("\n### Skipping plot of range edge; only possible if `getJointPrecision=TRUE`, `Options['Calculate_Range']=TRUE`, and `n_samples`>0")
}
# Plot densities
message("\n### Making plots of spatial predictions")
# Including Report as input to plot_maps_args so that it doesn't need to re-run it when loading results without TMB linked
# Not restricting to named arguments using args_to_use so that ... passes settings to plot_variable even that aren't formal arguments to plot_maps
plot_maps_args = list(...)
plot_maps_args = combine_lists( "input"=plot_maps_args, # "args_to_use"=formalArgs(plot_maps)
"default" = list( fit = fit,
plot_set = plot_set,
category_names = category_names,
PlotDF = map_list[["PlotDF"]],
MapSizeRatio = map_list[["MapSizeRatio"]],
working_dir = working_dir,
year_labels = year_labels,
years_to_plot = years_to_plot,
legend_x = map_list[["Legend"]]$x/100,
legend_y = map_list[["Legend"]]$y/100,
projargs = projargs,
n_cells = n_cells) )
Dens_xt = do.call( what=plot_maps, args=plot_maps_args )
# Plot factors
message("\n### Making plots for factors (if present)")
# Not restricting to named arguments using args_to_use so that ... and n_cells pass to plot_variable
plot_factors_args = list(...)
plot_factors_args = combine_lists( "input"=plot_factors_args, # "args_to_use"=formalArgs(plot_factors),
"default" = list( fit = fit,
mapdetails_list = map_list,
projargs = projargs,
n_cells = n_cells,
RotationMethod = RotationMethod,
plotdir = working_dir,
category_names = category_names ) )
Factors = do.call( what=plot_factors, args=plot_factors_args )
# Plot factors
if( cluster_results == TRUE ){
message("\n### Making plots for spatial cluster analysis")
plot_clusters_args = list(...)
plot_clusters_args = combine_lists( "input"=plot_clusters_args, "args_to_use"=c(formalArgs(plot_clusters),formalArgs(plot_variable)),
"default" = list( fit = fit,
year_labels = year_labels,
category_names = category_names,
map_list = map_list,
working_dir = working_dir,
n_cells = n_cells,
projargs = projargs ) )
Clusters = do.call( what=plot_clusters, args=plot_clusters_args )
}else{
Clusters = NULL
}
# Plot quantile-quantile plot
if( check_residuals == TRUE ){
# Plot diagnostic for encounter probability
#message("\n### Making plot of encounter probability")
#Enc_prob = plot_encounter_diagnostic( Report=fit$Report, Data_Geostat=cbind("Catch_KG"=fit$data_frame[,'b_i']), DirName=working_dir)
#
#message("\n### Making Q-Q plot")
#Q = plot_quantile_diagnostic( TmbData=fit$data_list, Report=fit$Report, FileName_PP="Posterior_Predictive",
# FileName_Phist="Posterior_Predictive-Histogram", FileName_QQ="Q-Q_plot", FileName_Qhist="Q-Q_hist", save_dir=working_dir )
#
## Pearson residuals
#message("\n### Making plot of Pearson residuals")
#plot_residuals(Lat_i=fit$data_frame[,'Lat_i'], Lon_i=fit$data_frame[,'Lon_i'], TmbData=fit$data_list, Report=fit$Report,
# Q=Q, working_dir=working_dir, spatial_list=fit$spatial_list, extrapolation_list=fit$extrapolation_list,
# year_labels=year_labels, years_to_plot=years_to_plot, zrange=zrange,
# legend_x=map_list[["Legend"]]$x/100, legend_y=map_list[["Legend"]]$y/100 )
# Plotting quantile residuals
message("\n### Making quantile residuals using conditional simulation and package DHARMa")
dharmaRes = summary( fit,
what = "residuals",
working_dir = working_dir,
type = type,
... )
# Mapping quantile residuals
message("\n### Plotting quantile residuals ")
dharma_raster = plot_quantile_residuals( dharmaRes = dharmaRes,
fit = fit,
working_dir = working_dir,
year_labels = year_labels,
years_to_plot = years_to_plot,
n_cells_residuals = n_cells_residuals,
projargs = projargs,
... )
# Semivariance for quantile residuals
# Disabled due to problems with raster plots in V >= 4.2.1
#if( fit$data_list$n_t > 1 ){
#message("\n### Plotting semivariance for normal-transformed quantile residuals ")
#residual_semivariance = plot_residual_semivariance( fit = fit,
# dharma_raster = dharma_raster,
# dharmaRes = dharmaRes,
# working_dir = working_dir )
#}else{
message("\n### Skipping plot of semivariance for normal-transformed quantile residuals")
residual_semivariance = NULL
#}
}else{
#Q = "Not run"
#message("\n### Skipping Q-Q plot")
#message("\n### Skipping plot of Pearson residuals")
message("\n### Skipping quantile residuals using conditional simulation and package DHARMa")
message("\n### Skipping plot of quantile residuals ")
message("\n### Skipping plot of semivariance for normal-transformed quantile residuals")
dharmaRes = NULL
dharma_raster = NULL
residual_semivariance = NULL
}
# return
Return = list( "dharmaRes" = dharmaRes,
"dharma_raster" = dharma_raster,
"residual_semivariance" = residual_semivariance,
"Index" = Index,
"Proportions" = Proportions,
"Range" = Range,
"Dens_xt" = Dens_xt,
"Edge" = Edge,
"map_list" = map_list,
"plot_maps_args" = plot_maps_args,
"plot_biomass_index_args" = plot_biomass_index_args,
"Factors" = Factors,
"Clusters" = Clusters )
return( invisible(Return) )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.