R/plot_results.R

In James-Thorson/FishStatsUtils: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox

#' 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[FishStatsUtils]{make_settings}} for generic settings, \code{\link[FishStatsUtils]{fit_model}} for model fitting, and \code{\link[FishStatsUtils]{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) )
}