From c4db9d75246d48e5af9161d981753e7285bedd71 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Fri, 7 Aug 2020 16:30:16 -0700 Subject: [PATCH 01/15] Expanding formula interface --- DESCRIPTION | 6 +- NAMESPACE | 3 + R/Calc_Polygon_Areas_and_Polygons_Fn.R | 2 +- R/calculate_proportion.R | 2 +- R/convert_shapefile.R | 2 +- R/fit_model.R | 126 +++++++++++++++++++--- R/make_covariates.R | 20 +++- R/make_extrapolation_info.R | 9 +- R/{Calc_Kmeans.R => make_kmeans.R} | 21 +++- R/make_mesh.R | 112 +++++++++++++++++++ R/make_settings.R | 2 +- R/make_spatial_info.R | 25 +++-- R/plot_results.R | 2 +- R/project_coordinates.R | 2 +- R/zzz.R | 27 ++++- man/Calc_Anisotropic_Mesh.Rd | 16 ++- man/Calc_Kmeans.Rd | 35 ++---- man/Calc_Polygon_Areas_and_Polygons_Fn.Rd | 8 +- man/Coherence.Rd | 8 +- man/Crossvalidate_Fn.Rd | 14 ++- man/FishStatsUtils.Rd | 1 - man/Geostat_Sim.Rd | 11 +- man/PlotMap_Fn.Rd | 37 +++++-- man/Rerun_Fn.Rd | 16 ++- man/calc_coherence.Rd | 8 +- man/calc_synchrony.Rd | 3 +- man/calculate_proportion.Rd | 25 +++-- man/combine_lists.Rd | 3 +- man/convert_shapefile.Rd | 17 ++- man/fit_model.Rd | 98 ++++++++++++----- man/format_covariates.Rd | 14 ++- man/make_covariates.Rd | 9 +- man/make_extrapolation_info.Rd | 36 +++++-- man/make_kmeans.Rd | 46 ++++++++ man/make_map_info.Rd | 9 +- man/make_mesh.Rd | 32 ++++++ man/make_settings.Rd | 40 +++++-- man/make_spatial_info.Rd | 37 +++++-- man/plot.make_extrapolation_info.Rd | 3 +- man/plot_biomass_index.Rd | 26 +++-- man/plot_cov.Rd | 11 +- man/plot_data.Rd | 27 +++-- man/plot_encounter_diagnostic.Rd | 12 ++- man/plot_factors.Rd | 22 +++- man/plot_index.Rd | 37 +++++-- man/plot_loadings.Rd | 16 ++- man/plot_maps.Rd | 29 +++-- man/plot_overdispersion.Rd | 14 ++- man/plot_quantile_diagnostic.Rd | 11 +- man/plot_quantile_residuals.Rd | 11 +- man/plot_range_edge.Rd | 22 +++- man/plot_range_index.Rd | 25 +++-- man/plot_residuals.Rd | 19 +++- man/plot_results.Rd | 31 ++++-- man/plot_variable.Rd | 37 +++++-- man/predict.fit_model.Rd | 47 ++++++++ man/project_coordinates.Rd | 12 ++- man/rotate_factors.Rd | 10 +- man/sample_variable.Rd | 10 +- man/smallPlot.Rd | 19 +++- man/summarize_covariance.Rd | 21 +++- man/summary.fit_model.Rd | 10 +- 62 files changed, 1089 insertions(+), 277 deletions(-) rename R/{Calc_Kmeans.R => make_kmeans.R} (71%) create mode 100644 R/make_mesh.R create mode 100644 man/make_kmeans.Rd create mode 100644 man/make_mesh.Rd create mode 100644 man/predict.fit_model.Rd diff --git a/DESCRIPTION b/DESCRIPTION index f7c21ba..15f82b8 100644 --- 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.7.0 -Date: 2020-04-23 +Version: 2.8.0 +Date: 2020-08-07 Authors@R: c(person(given = "James", family = "Thorson", @@ -53,6 +53,6 @@ Remotes: License: GPL-3 LazyData: yes BuildVignettes: yes -RoxygenNote: 7.0.2 +RoxygenNote: 7.1.1 URL: http://github.com/james-thorson/FishStatsUtils BugReports: http://github.com/james-thorson/FishStatsUtils/issues diff --git a/NAMESPACE b/NAMESPACE index 0ea1302..e83f2d8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ S3method(plot,fit_model) S3method(plot,make_extrapolation_info) S3method(plot,make_spatial_info) +S3method(predict,fit_model) S3method(print,fit_model) S3method(print,make_extrapolation_info) S3method(print,make_spatial_info) @@ -44,7 +45,9 @@ export(get_latest_version) export(load_example) export(make_covariates) export(make_extrapolation_info) +export(make_kmeans) export(make_map_info) +export(make_mesh) export(make_settings) export(make_spatial_info) export(map_hypervariance) diff --git a/R/Calc_Polygon_Areas_and_Polygons_Fn.R b/R/Calc_Polygon_Areas_and_Polygons_Fn.R index 587a066..6e44531 100644 --- a/R/Calc_Polygon_Areas_and_Polygons_Fn.R +++ b/R/Calc_Polygon_Areas_and_Polygons_Fn.R @@ -3,10 +3,10 @@ #' #' \code{Calc_Polygon_Areas_and_Polygons_Fn} builds outputs for a given triangulated mesh used for approximating spatial variation #' +#' @inheritParams Calc_Anisotropic_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) #' @param a_el a matrix with \code{nrow(Data_Extrap)} rows and l columns, giving the area within the l-th stratum for the e-th row of Data_Extrap -#' @inheritParams Calc_Anisotropic_Mesh #' @return Tagged list containing distance metrics #' \describe{ diff --git a/R/calculate_proportion.R b/R/calculate_proportion.R index b2df0db..5c447d6 100644 --- a/R/calculate_proportion.R +++ b/R/calculate_proportion.R @@ -4,9 +4,9 @@ #' #' \code{calculate_proportion} takes output from a VAST run and calculates the proportion of biomass in different categories #' -#' @param Index output from \code{FishStatsUtils::plot_biomass_index} #' @inheritParams plot_biomass_index #' @inheritParams VAST::make_data +#' @param Index output from \code{FishStatsUtils::plot_biomass_index} #' @param ... list of arguments to pass to \code{plot_index} #' #' @return Tagged list of output diff --git a/R/convert_shapefile.R b/R/convert_shapefile.R index 7d59247..6c29477 100644 --- a/R/convert_shapefile.R +++ b/R/convert_shapefile.R @@ -5,7 +5,7 @@ #' \code{convert_shapefile} reads in a shapefile and creates an extrapolation with a chosen resolution #' #' @inheritParams sp::CRS -#' @inheritParams make_extraploation_info +#' @inheritParams make_extrapolation_info #' #' @param file_path path for shapefile on harddrive #' @param make_plots Boolean indicating whether to visualize inputs and outputs as maps diff --git a/R/fit_model.R b/R/fit_model.R index 4d8fc47..7d54cd1 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -55,7 +55,7 @@ #' } #' #' @family wrapper functions -#' @seealso \code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots +#' @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 #' @seealso \code{\link{summary.fit_model}} for methods to summarize output, including obtain a dataframe of estimated densities #' #' @examples @@ -70,14 +70,20 @@ #' example = load_example( data_set="EBS_pollock" ) #' #' # Make settings -#' settings = make_settings( n_x=50, Region=example$Region, purpose="index", -#' strata.limits=example$strata.limits ) +#' settings = make_settings( n_x=50, +#' Region=example$Region, +#' purpose="index", +#' strata.limits=example$strata.limits ) #' #' # Run model -#' fit = fit_model( "settings"=settings, "Lat_i"=example$sampling_data[,'Lat'], -#' "Lon_i"=example$sampling_data[,'Lon'], "t_i"=example$sampling_data[,'Year'], -#' "c_i"=rep(0,nrow(example$sampling_data)), "b_i"=example$sampling_data[,'Catch_KG'], -#' "a_i"=example$sampling_data[,'AreaSwept_km2'], "v_i"=example$sampling_data[,'Vessel'] ) +#' fit = fit_model( "settings"=settings, +#' "Lat_i"=example$sampling_data[,'Lat'], +#' "Lon_i"=example$sampling_data[,'Lon'], +#' "t_i"=example$sampling_data[,'Year'], +#' "c_i"=rep(0,nrow(example$sampling_data)), +#' "b_i"=example$sampling_data[,'Catch_KG'], +#' "a_i"=example$sampling_data[,'AreaSwept_km2'], +#' "v_i"=example$sampling_data[,'Vessel'] ) #' #' # Plot results #' plot_results( settings=settings, fit=fit ) @@ -88,8 +94,9 @@ # Using https://cran.r-project.org/web/packages/roxygen2/vignettes/rd-formatting.html for guidance on markdown-enabled documentation fit_model = function( settings, Lat_i, Lon_i, t_i, b_i, a_i, c_iz=rep(0,length(b_i)), v_i=rep(0,length(b_i)), working_dir=paste0(getwd(),"/"), - Xconfig_zcp=NULL, covariate_data, formula=~0, Q_ik=NULL, newtonsteps=1, - silent=TRUE, build_model=TRUE, run_model=TRUE, test_fit=TRUE, ... ){ + X1config_cp=NULL, X2config_cp=NULL, covariate_data, X1_formula=~0, X2_formula=~0, + Q1config_k=NULL, Q2config_k=NULL, catchability_data, Q1_formula=~0, Q2_formula=~0, + newtonsteps=1, silent=TRUE, build_model=TRUE, run_model=TRUE, test_fit=TRUE, ... ){ # Capture extra arguments to function extra_args = list(...) @@ -125,12 +132,15 @@ fit_model = function( settings, Lat_i, Lon_i, t_i, b_i, a_i, c_iz=rep(0,length(b # Do *not* restrict inputs to formalArgs(make_data) because other potential inputs are still parsed by make_data for backwards compatibility message("\n### Making data object") # VAST:: if(missing(covariate_data)) covariate_data = NULL + if(missing(catchability_data)) catchability_data = NULL data_args_default = list("Version"=settings$Version, "FieldConfig"=settings$FieldConfig, "OverdispersionConfig"=settings$OverdispersionConfig, "RhoConfig"=settings$RhoConfig, "VamConfig"=settings$VamConfig, "ObsModel"=settings$ObsModel, "c_iz"=c_iz, "b_i"=b_i, "a_i"=a_i, "v_i"=v_i, "s_i"=spatial_list$knot_i-1, "t_i"=t_i, "spatial_list"=spatial_list, "Options"=settings$Options, "Aniso"=settings$use_anisotropy, - Xconfig_zcp=Xconfig_zcp, covariate_data=covariate_data, formula=formula, Q_ik=Q_ik) + "X1config_cp"=X1config_cp, "X2config_cp"=X2config_cp, "covariate_data"=covariate_data, "X1_formula"=X1_formula, "X2_formula"=X2_formula, + "Q1config_k"=Q1config_k, "Q2config_k"=Q2config_k, "catchability_data"=catchability_data, "Q1_formula"=Q1_formula, "Q2_formula"=Q2_formula ) data_args_input = combine_lists( input=extra_args, default=data_args_default ) # Do *not* use args_to_use data_list = do.call( what=make_data, args=data_args_input ) + #return(data_list) } # Build object message("\n### Making TMB object") @@ -216,8 +226,9 @@ fit_model = function( settings, Lat_i, Lon_i, t_i, b_i, a_i, c_iz=rep(0,length(b "data_args_input"=data_args_input ) Return = list("data_frame"=data_frame, "extrapolation_list"=extrapolation_list, "spatial_list"=spatial_list, "data_list"=data_list, "tmb_list"=tmb_list, "parameter_estimates"=parameter_estimates, "Report"=Report, - "ParHat"=ParHat, "year_labels"=year_labels, "years_to_plot"=years_to_plot, "settings"=settings, - "input_args"=input_args ) + "ParHat"=ParHat, "year_labels"=year_labels, "years_to_plot"=years_to_plot, "settings"=settings, "input_args"=input_args, + "X1config_cp"=X1config_cp, "X2config_cp"=X2config_cp, "covariate_data"=covariate_data, "X1_formula"=X1_formula, "X2_formula"=X2_formula, + "Q1config_k"=Q1config_k, "Q2config_k"=Q1config_k, "catchability_data"=catchability_data, "Q1_formula"=Q1_formula, "Q2_formula"=Q2_formula) class(Return) = "fit_model" return( Return ) } @@ -408,4 +419,95 @@ summary.fit_model <- function(x, what="density", n_samples=250, } +#' Predict density for new samples (\emph{Beta version; may change without notice}) +#' +#' \code{predict.fit_model} calculates predictions given new data +#' +#' When the user does not supply \code{covariate_data}, predictions are based upon covariate values +#' interpolated from covariates supplied when fitting the model. When the user does supply \code{covariate_data}, +#' predictions are created when interpolating based on those new supplied values. +#' +#' @inheritParams make_covariates +#' @inheritParams VAST::make_data +#' @param x Output from \code{\link{fit_model}} +#' @param what Which output from \code{fit$Report} should be extracted; default is predicted density +#' +#' @return NULL +#' @method predict fit_model +#' @export +predict.fit_model <- function(x, what="D_i", Lat_i, Lon_i, t_i, a_i, c_iz=rep(0,length(t_i)), + v_i=rep(0,length(t_i)), covariate_data=NULL, catchability_data=NULL, + working_dir=paste0(getwd(),"/")) +{ + warning("`predict.fit_model(.)` is still in development") + + # Process inputs + PredTF_i = c( x$data_list$PredTF_i, rep(1,length(t_i)) ) + b_i = c( x$data_frame[,"b_i"], rep(1,length(t_i)) ) + c_iz = rbind( matrix(x$data_frame[,grep("c_iz",names(x$data_frame))]), matrix(c_iz) ) + Lat_i = c( x$data_frame[,"Lat_i"], Lat_i ) + Lon_i = c( x$data_frame[,"Lon_i"], Lon_i ) + a_i = c( x$data_frame[,"a_i"], a_i ) + v_i = c( x$data_frame[,"v_i"], v_i ) + t_i = c( x$data_frame[,"t_i"], t_i ) + catchability_data = rbind( x$catchability_data, catchability_data ) + assign("b_i", b_i, envir=.GlobalEnv) + + # Populate covariate_data + # Default: use original covariate values + # When user provides new values, deal with missing years in covariate_data so that it doesn't throw an unnecessary error + if( is.null(covariate_data) ){ + message("Using `covariate_data` supplied during original fit to interpolate covariate values for predictions") + covariate_data = x$covariate_data + }else{ + if( !(all(t_i %in% covariate_data[,'Year']) | any(is.na(covariate_data[,'Year']))) ){ + stop("Some `t_i` values are supplied without covariates") + }else{ + message("Covariates not provided for all modeled years, so filling in zeros for covariates in other years") + zeros_data = covariate_data[match(covariate_data[,'Year'],unique(covariate_data[,'Year'])),] + zeros_data[,setdiff(names(zeros_data),c('Lat','Lon','Year'))] = 0 + zeros_data[,c('Lat','Lon')] = Inf + covariate_data = rbind( covariate_data, zeros_data) + } + } + + # Build information regarding spatial location and correlation + message("\n### Re-making spatial information") + spatial_args_new = list("anisotropic_mesh"=x$spatial_list$MeshList$anisotropic_mesh, "Kmeans"=x$spatial_list$Kmeans, + "Lon_i"=Lon_i, "Lat_i"=Lat_i ) + spatial_args_input = combine_lists( input=spatial_args_new, default=x$input_args$spatial_args_input ) + spatial_list = do.call( what=make_spatial_info, args=spatial_args_input ) + + # Build data + # Do *not* restrict inputs to formalArgs(make_data) because other potential inputs are still parsed by make_data for backwards compatibility + message("\n### Re-making data object") + data_args_new = list( "c_iz"=c_iz, "b_i"=b_i, "a_i"=a_i, "v_i"=v_i, "PredTF_i"=PredTF_i, + "t_i"=t_i, "spatial_list"=spatial_list, + "covariate_data"=covariate_data, "catchability_data"=catchability_data ) + data_args_input = combine_lists( input=data_args_new, default=x$input_args$data_args_input ) # Do *not* use args_to_use + data_list = do.call( what=make_data, args=data_args_input ) + data_list$n_g = 0 + + # Build object + message("\n### Re-making TMB object") + model_args_default = list("TmbData"=data_list, "RunDir"=working_dir, "Version"=x$settings$Version, + "RhoConfig"=x$settings$RhoConfig, "loc_x"=spatial_list$loc_x, "Method"=spatial_list$Method) + model_args_input = combine_lists( input=list("Parameters"=x$ParHat), + default=model_args_default, args_to_use=formalArgs(make_model) ) + tmb_list = do.call( what=make_model, args=model_args_input ) + + # Extract output + Report = tmb_list$Obj$report() + Y_i = Report[[what]][(1+nrow(x$data_frame)):length(Report$D_i)] + + # sanity check + if( all.equal(covariate_data,x$covariate_data) & Report$jnll!=x$Report$jnll){ + message("Problem detected in `predict.fit_model`; returning outputs for diagnostic purposes") + Return = list("Report"=Report, "data_list"=data_list) + return(Return) + } + + # return prediction + return(Y_i) +} diff --git a/R/make_covariates.R b/R/make_covariates.R index 1fd899b..2e03b9d 100644 --- a/R/make_covariates.R +++ b/R/make_covariates.R @@ -46,9 +46,10 @@ make_covariates = function( formula, covariate_data, Year_i, spatial_list, extra # Create data frame of necessary size DF_zp = NULL - DF_ip = cbind( sample_data, covariate_data[rep(1,nrow(sample_data)),covariate_names] ) + #DF_ip = cbind( sample_data, covariate_data[rep(1,nrow(sample_data)),covariate_names] ) + DF_ip = data.frame( sample_data, covariate_data[rep(1,nrow(sample_data)),covariate_names] ) colnames(DF_ip) = c( names(sample_data) ,covariate_names ) - DF_ip[,covariate_names] = NA + #DF_ip[,covariate_names] = NA # Loop through data and extrapolation-grid for( tI in seq_along(Year_Set) ){ @@ -78,7 +79,17 @@ make_covariates = function( formula, covariate_data, Year_i, spatial_list, extra # Convert to dimensions requested DF = rbind( DF_ip, DF_zp ) - X = model.matrix( update.formula(formula, ~.+0), data=DF )[,,drop=FALSE] + + # Make model.matrix + # To ensure identifiability given betas (intercepts), add intercept to formula + # and then remove that term from model.matrix. This will not fix identifiability + # issues arising when both conditions are met: + # factor(Year) has an interaction with another factor, and + # betas vary among years (are not constant) + Model_matrix = model.matrix( update.formula(formula, ~.+1), data=DF ) + Columns_to_keep = which( attr(Model_matrix,"assign") != 0 ) + coefficient_names = attr(Model_matrix,"dimnames")[[2]][Columns_to_keep] + X = Model_matrix[,Columns_to_keep,drop=FALSE] # Make X_ip X_ip = X[ 1:nrow(DF_ip), , drop=FALSE ] @@ -102,7 +113,8 @@ make_covariates = function( formula, covariate_data, Year_i, spatial_list, extra } # return stuff - Return = list( "X_gtp"=X_gtp, "X_itp"=X_itp, "covariate_names"=covariate_names ) + Return = list( "X_gtp"=X_gtp, "X_itp"=X_itp, "covariate_names"=covariate_names, + "coefficient_names"=coefficient_names ) return( Return ) } diff --git a/R/make_extrapolation_info.R b/R/make_extrapolation_info.R index f018248..d9315f8 100644 --- a/R/make_extrapolation_info.R +++ b/R/make_extrapolation_info.R @@ -19,7 +19,7 @@ #' projection errors regarding total area than rnaturalearth. #' #' @inheritParams sp::CRS -#' @inheritParams Calc_Kmeans +#' @inheritParams make_kmeans #' @inheritParams convert_shapefile #' #' @param Region a character vector, where each element that is matched against potential values to determine the region for the extrapolation grid. Current options are "california_current", "west_coast_hook_and_line", "british_columbia", "eastern_bering_sea", "northern_bering_sea", "bering_sea_slope", "st_matthews_island", "aleutian_islands", "gulf_of_alaska", "northwest_atlantic", "south_africa", "gulf_of_st_lawrence", "new_zealand", "habcam", "gulf_of_mexico", "ATL-IBTS-Q1", "ATL-IBTS-Q4", "BITS", "BTS", "BTS-VIIA", "EVHOE", "IE-IGFS", "NIGFS", "NS_IBTS", "PT-IBTS", "SP-ARSA", "SP-NORTH", "SP-PORC", "stream_network", "user", "other", or the absolute path and file name for a GIS shapefile @@ -57,7 +57,7 @@ make_extrapolation_info = function( Region, projargs=NA, zone=NA, strata.limits= region=c("south_coast","west_coast"), strata_to_use=c('SOG','WCVI','QCS','HS','WCHG'), epu_to_use=c('All','Georges_Bank','Mid_Atlantic_Bight','Scotian_Shelf','Gulf_of_Maine','Other')[1], survey="Chatham_rise", surveyname='propInWCGBTS', flip_around_dateline, nstart=100, - area_tolerance=0.05, ... ){ + area_tolerance=0.05, backwards_compatible_kmeans=FALSE, ... ){ # Note: flip_around_dateline must appear in arguments for argument-matching in fit_model # However, it requires a different default value for different regions; hence the input format being used. @@ -183,8 +183,9 @@ make_extrapolation_info = function( Region, projargs=NA, zone=NA, strata.limits= # Run K-means only on grid-cells with nonzero area loc_orig = Return$Data_Extrap[,c("E_km","N_km")] loc_orig = loc_orig[ which(Return$Area_km2_x>0), ] - Kmeans = Calc_Kmeans( n_x=max_cells, loc_orig=loc_orig, nstart=nstart, - randomseed=1, iter.max=1000, DirPath=paste0(getwd(),"/"), Save_Results=FALSE ) + Kmeans = make_kmeans( n_x=max_cells, loc_orig=loc_orig, nstart=nstart, + randomseed=1, iter.max=1000, DirPath=paste0(getwd(),"/"), Save_Results=FALSE, + backwards_compatible_kmeans=backwards_compatible_kmeans ) Kmeans[["cluster"]] = RANN::nn2( data=Kmeans[["centers"]], query=Return$Data_Extrap[,c("E_km","N_km")], k=1)$nn.idx[,1] # Transform Extrapolation_List aggregate_vector = function( values_x, index_x, max_index, FUN=sum ){ diff --git a/R/Calc_Kmeans.R b/R/make_kmeans.R similarity index 71% rename from R/Calc_Kmeans.R rename to R/make_kmeans.R index df01c28..75c4608 100644 --- a/R/Calc_Kmeans.R +++ b/R/make_kmeans.R @@ -1,7 +1,7 @@ #' Calculate location for knots approximating spatial variation #' -#' \code{Calc_Kmeans} determines the location for a set of knots for approximating spatial variation +#' \code{make_kmeans} determines the location for a set of knots for approximating spatial variation #' #' @param n_x the number of knots to select #' @param loc_orig a matrix with two columns where each row gives the 2-dimensional coordinates to be approximated @@ -10,6 +10,9 @@ #' @param iter.max the number of iterations used per k-means algorithm (default=1000) #' @param DirPath a directory where the algorithm looks for a previously-saved output (default is working directory) #' @param Save_Results a boolean stating whether to save the output (Default=TRUE) +#' @param backwards_compatible_kmeans a boolean stating how to deal with changes in the kmeans algorithm implemented in R version 3.6.0, +#' where \code{backwards_compatible_kmeans==TRUE} modifies the default algorithm to maintain backwards compatibility, and +#' where \code{backwards_compatible_kmeans==FALSE} breaks backwards compatibility between R versions prior to and after R 3.6.0. #' @return Tagged list containing outputs #' \describe{ @@ -18,8 +21,9 @@ #' } #' @export -Calc_Kmeans <- -function( n_x, loc_orig, nstart=100, randomseed=1, iter.max=1000, DirPath=paste0(getwd(),"/"), Save_Results=TRUE ){ +make_kmeans <- +function( n_x, loc_orig, nstart=100, randomseed=1, iter.max=1000, DirPath=paste0(getwd(),"/"), + Save_Results=TRUE, backwards_compatible_kmeans=FALSE ){ # get old seed oldseed = ceiling(runif(1,min=1,max=1e6)) @@ -29,6 +33,17 @@ function( n_x, loc_orig, nstart=100, randomseed=1, iter.max=1000, DirPath=paste0 options( "warn" = -1 ) on.exit( options(old.options) ) + # Backwards compatibility + if( backwards_compatible_kmeans==TRUE ){ + if( identical(formalArgs(RNGkind), c("kind","normal.kind","sample.kind")) ){ + RNGkind_orig = RNGkind() + on.exit( RNGkind(kind=RNGkind_orig[1], normal.kind=RNGkind_orig[2], sample.kind=RNGkind_orig[3]), add=TRUE ) + RNGkind( sample.kind="Rounding" ) + }else if( !identical(formalArgs(RNGkind), c("kind","normal.kind")) ){ + stop("Assumptions about `RNGkind` are not met within `make_kmeans`; please report problem to package developers") + } + } + # Calculate knots for SPDE mesh if( length(unique(paste(loc_orig[,1],loc_orig[,2],sep="_")))<=n_x ){ # If number of knots is less than number of sample locations diff --git a/R/make_mesh.R b/R/make_mesh.R new file mode 100644 index 0000000..ccff78f --- /dev/null +++ b/R/make_mesh.R @@ -0,0 +1,112 @@ + +#' Make mesh for distances among points +#' +#' \code{make_mesh} builds a tagged list representing distances for isotropic or geometric anisotropic triangulated mesh +#' +#' @param loc_x location (eastings and northings in kilometers, UTM) for each sample or knot +#' @param Method spatial method determines ("Mesh" and "Grid" give +#' @param anisotropic_mesh OPTIONAL, anisotropic mesh (if missing, its recalculated from loc_x) +#' @param ... Arguments passed to \code{INLA::inla.mesh.create} + +#' @return Tagged list containing distance metrics + +#' @export +make_mesh <- +function(loc_x, loc_g, loc_i, Method, Extrapolation_List, anisotropic_mesh=NULL, fine_scale=FALSE, ...){ + + ####################### + # Create the anisotropic SPDE mesh using 2D coordinates + ####################### + + # 2D coordinates SPDE + if( is.null(anisotropic_mesh)){ + if( fine_scale==FALSE ){ + anisotropic_mesh = INLA::inla.mesh.create( loc_x, plot.delay=NULL, ...) + }else{ + loc_z = rbind( loc_x, loc_g, loc_i ) + outer_hull = INLA::inla.nonconvex.hull( as.matrix(loc_z), convex = -0.05, concave = -0.05) + anisotropic_mesh = INLA::inla.mesh.create( loc_x, plot.delay=NULL, boundary=outer_hull, ...) + } + } + + anisotropic_spde = INLA::inla.spde2.matern(anisotropic_mesh, alpha=2) + + # Pre-processing in R for anisotropy + Dset = 1:2 + # Triangle info + TV = anisotropic_mesh$graph$tv # Triangle to vertex indexing + V0 = anisotropic_mesh$loc[TV[,1],Dset] # V = vertices for each triangle + V1 = anisotropic_mesh$loc[TV[,2],Dset] + V2 = anisotropic_mesh$loc[TV[,3],Dset] + E0 = V2 - V1 # E = edge for each triangle + E1 = V0 - V2 + E2 = V1 - V0 + + # Pre-processing for barriers + # Barriers don't affect projection matrix A + # Obtain polygon for water + map_data = rnaturalearth::ne_countries( scale=switch("medium", "low"=110, "medium"=50, "high"=10, 50) ) + attr(map_data,"proj4string") = sp::CRS("+proj=longlat +datum=WGS84") + + # Calculate centroid of each triangle in mesh and convert to SpatialPoints + n_triangles = length(anisotropic_mesh$graph$tv[,1]) + posTri = matrix(NA, nrow=n_triangles, ncol=2) + for(tri_index in 1:n_triangles){ + temp = anisotropic_mesh$loc[ anisotropic_mesh$graph$tv[tri_index,], ] + posTri[tri_index,] = colMeans(temp)[c(1,2)] + } + posTri = sp::SpatialPoints(posTri, proj4string=sp::CRS(Extrapolation_List$projargs) ) + posTri = sp::spTransform(posTri, CRSobj=map_data@proj4string ) + + # Calculate set of triangles barrier.triangles with centroid over land + if( Method == "Barrier" ){ + anisotropic_mesh_triangles_over_land = unlist(sp::over(map_data, posTri, returnList=TRUE)) + }else{ + anisotropic_mesh_triangles_over_land = vector() + } + # + #plot( x=posTri@coords[,1], y=posTri@coords[,2], col=ifelse(1:n_triangles%in%triangles_over_land,"black","red") ) + + # Create Barrier object if requested + # Don't do this unless necessary, because it sometimes throws an error + #Diagnose issues: assign("anisotropic_mesh", anisotropic_mesh, envir = .GlobalEnv) + barrier_finite_elements = INLA:::inla.barrier.fem(mesh=anisotropic_mesh, + barrier.triangles=anisotropic_mesh_triangles_over_land) + barrier_list = list(C0 = barrier_finite_elements$C[[1]], + C1 = barrier_finite_elements$C[[2]], + D0 = barrier_finite_elements$D[[1]], + D1 = barrier_finite_elements$D[[2]], + I = barrier_finite_elements$I ) + # sp::plot( INLA::inla.barrier.polygon(anisotropic_mesh, triangles_over_land) ) + + # Calculate Areas + crossprod_fn = function(Vec1,Vec2) abs(det( rbind(Vec1,Vec2) )) + Tri_Area = rep(NA, nrow(E0)) + for(i in 1:length(Tri_Area)) Tri_Area[i] = crossprod_fn( E0[i,],E1[i,] )/2 # T = area of each triangle + + ################ + # Add the isotropic SPDE mesh for spherical or 2D projection, depending upon `Method` input + ################ + + # Mesh and SPDE for different inputs + if(Method %in% c("Mesh","Grid","Stream_network","Barrier")){ + loc_isotropic_mesh = loc_x + isotropic_mesh = anisotropic_mesh + } + if(Method %in% c("Spherical_mesh")){ + loc_isotropic_mesh = INLA::inla.mesh.map(loc_x, projection="longlat", inverse=TRUE) # Project from lat/long to mesh coordinates + isotropic_mesh = INLA::inla.mesh.create( loc_isotropic_mesh, plot.delay=NULL, ...) + } + isotropic_spde = INLA::inla.spde2.matern(isotropic_mesh, alpha=2) + + #################### + # Return stuff + #################### + #if( isotropic_mesh$n != anisotropic_mesh$n ) stop("Check `Calc_Anisotropic_Mesh` for problem") + + Return = list("loc_x"=loc_x, "loc_isotropic_mesh"=loc_isotropic_mesh, "isotropic_mesh"=isotropic_mesh, + "isotropic_spde"=isotropic_spde, "anisotropic_mesh"=anisotropic_mesh, "anisotropic_spde"=anisotropic_spde, + "Tri_Area"=Tri_Area, "TV"=TV, "E0"=E0, "E1"=E1, "E2"=E2, + "anisotropic_mesh_triangles_over_land"=anisotropic_mesh_triangles_over_land, "barrier_list"=barrier_list ) + return(Return) +} diff --git a/R/make_settings.R b/R/make_settings.R index c3c6a2e..7d39815 100644 --- a/R/make_settings.R +++ b/R/make_settings.R @@ -19,7 +19,7 @@ #' @return Tagged list containing default settings for a given purpose, use \code{names} on output to see list of settings. #' #' @family wrapper functions -#' @seealso \code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots +#' @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 make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE, diff --git a/R/make_spatial_info.R b/R/make_spatial_info.R index 5a71ec5..b81043b 100644 --- a/R/make_spatial_info.R +++ b/R/make_spatial_info.R @@ -14,6 +14,9 @@ #' #' \code{LON_intensity} and \code{LAT_intensity} allow users to specify locations that are used by the k-means algorithm to determine the location of knots, where e.g. users can either hard-code the desired knot locations via these inputs (using \code{n_x} greater than this number of locations), or use the extrapolation-grid to ensure that knots are located proportional to that grid. #' +#' @inheritParams make_kmeans +#' @inheritParams make_mesh +#' #' @param n_x the number of vertices in the SPDE mesh (determines the spatial resolution when Method="Mesh") #' @param Lon_i Longitude for each sample #' @param Lat_i Latitude for each sample @@ -25,7 +28,6 @@ #' @param grid_size_LL the distance between grid cells for the 2D AR1 grid (determines spatial resolution when Method="Grid") when using \code{Method="Spherical_mesh"} #' @param Network_sz_LL data frame with "parent_s", "child_s", "dist_s", "Lat", "Lon", default=NULL only needed with Method == "Stream_network" #' @param ... additional arguments passed to \code{\link[INLA]{inla.mesh.create}} -#' @inheritParams Calc_Kmeans #' @return Tagged list containing objects for running a VAST model #' \describe{ @@ -33,7 +35,7 @@ #' \item{GridList}{A tagged list with inputs related to the 2D AR1 grid} #' \item{a_xl}{A data frame with areas for each knot and each strattum} #' \item{loc_UTM}{A data frame with the converted UTM coordinates for each sample} -#' \item{Kmeans}{Output from \code{Calc_Kmeans} with knots for a triangulated mesh} +#' \item{Kmeans}{Output from \code{make_kmeans} with knots for a triangulated mesh} #' \item{knot_i}{The knot associated with each sample} #' \item{Method}{The Method input (for archival purposes)} #' \item{loc_x}{The UTM location for each knot} @@ -42,9 +44,9 @@ #' @export make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method=NULL, Method="Mesh", - grid_size_km=50, grid_size_LL=1, fine_scale=FALSE, Network_sz_LL=NULL, + anisotropic_mesh=NULL, Kmeans=NULL, grid_size_km=50, grid_size_LL=1, fine_scale=FALSE, Network_sz_LL=NULL, iter.max=1000, randomseed=1, nstart=100, DirPath=paste0(getwd(),"/"), Save_Results=FALSE, - LON_intensity, LAT_intensity, ... ){ + LON_intensity, LAT_intensity, backwards_compatible_kmeans=FALSE, ... ){ # Deprecated options if( Method=="Spherical_mesh" ){ @@ -78,7 +80,8 @@ make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method Grid_bounds = (grid_size_km/110) * apply(loc_e/(grid_size_km/110), MARGIN=2, FUN=function(vec){trunc(range(vec))+c(0,1)}) # Calculate k-means centroids - Kmeans = Calc_Kmeans(n_x=n_x, loc_orig=loc_i[,c("Lon", "Lat")], randomseed=randomseed, ... ) + if(is.null(Kmeans)) Kmeans = make_kmeans(n_x=n_x, loc_orig=loc_i[,c("Lon", "Lat")], randomseed=randomseed, + backwards_compatible_kmeans=backwards_compatible_kmeans, ... ) # Calculate grid for 2D AR1 process loc_grid = expand.grid( 'Lon'=seq(Grid_bounds[1,1],Grid_bounds[2,1],by=grid_size_LL), 'Lat'=seq(Grid_bounds[1,2],Grid_bounds[2,2],by=grid_size_LL) ) @@ -94,7 +97,8 @@ make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method Grid_bounds = grid_size_km * apply(Extrapolation_List$Data_Extrap[,c('E_km','N_km')]/grid_size_km, MARGIN=2, FUN=function(vec){trunc(range(vec))+c(0,1)}) # Calculate k-means centroids - Kmeans = Calc_Kmeans(n_x=n_x, loc_orig=loc_intensity[,c("E_km", "N_km")], randomseed=randomseed, nstart=nstart, DirPath=DirPath, Save_Results=Save_Results ) + if(is.null(Kmeans)) Kmeans = make_kmeans(n_x=n_x, loc_orig=loc_intensity[,c("E_km", "N_km")], randomseed=randomseed, nstart=nstart, + DirPath=DirPath, Save_Results=Save_Results, backwards_compatible_kmeans=backwards_compatible_kmeans ) NN_i = RANN::nn2( data=Kmeans[["centers"]], query=loc_i, k=1)$nn.idx[,1] # Calculate grid for 2D AR1 process @@ -128,7 +132,8 @@ make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method } # Convert loc_x back to location in lat-long coordinates latlon_x - origargs = "+proj=longlat +ellps=WGS84" + #origargs = "+proj=longlat +ellps=WGS84" + origargs = "+proj=longlat +datum=WGS84" latlon_x = project_coordinates( X=loc_x[,"E_km"], Y=loc_x[,"N_km"], projargs=origargs, origargs=Extrapolation_List$projargs )[,c("Y","X")] colnames(latlon_x) = c("Lat", "Lon") @@ -142,9 +147,11 @@ make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method # Make mesh and info for anisotropy SpatialDeltaGLMM:: # Diagnose issues: assign("Kmeans", Kmeans, envir = .GlobalEnv) if(Method != "Stream_network"){ - MeshList = Calc_Anisotropic_Mesh( Method=Method, loc_x=Kmeans$centers, loc_g=loc_g, loc_i=loc_i, Extrapolation_List=Extrapolation_List, fine_scale=fine_scale, ... ) + MeshList = make_mesh( Method=Method, loc_x=Kmeans$centers, loc_g=loc_g, loc_i=loc_i, Extrapolation_List=Extrapolation_List, + fine_scale=fine_scale, anisotropic_mesh=anisotropic_mesh, ... ) }else{ - MeshList = Calc_Anisotropic_Mesh( Method=Method, loc_x=loc_x, loc_g=loc_g, loc_i=loc_i, Extrapolation_List=Extrapolation_List, fine_scale=fine_scale, ... ) + MeshList = make_mesh( Method=Method, loc_x=loc_x, loc_g=loc_g, loc_i=loc_i, Extrapolation_List=Extrapolation_List, + fine_scale=fine_scale, anisotropic_mesh=anisotropic_mesh, ... ) } n_s = switch( tolower(Method), "mesh"=MeshList$anisotropic_spde$n.spde, "grid"=nrow(loc_x), "spherical_mesh"=MeshList$isotropic_spde$n.spde, "stream_network"=nrow(loc_x), "barrier"=MeshList$anisotropic_spde$n.spde, ) diff --git a/R/plot_results.R b/R/plot_results.R index f1964b6..3f4b144 100644 --- a/R/plot_results.R +++ b/R/plot_results.R @@ -28,7 +28,7 @@ #' } #' #' @family wrapper functions -#' @seealso \code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots +#' @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=paste0(getwd(),"/"), diff --git a/R/project_coordinates.R b/R/project_coordinates.R index d7c19d0..37bf819 100644 --- a/R/project_coordinates.R +++ b/R/project_coordinates.R @@ -17,7 +17,7 @@ #' @export project_coordinates <- -function( X, Y, projargs=NA, origargs="+proj=longlat +ellps=WGS84", zone=NA, flip_around_dateline=FALSE, silent=FALSE ){ +function( X, Y, projargs=NA, origargs="+proj=longlat +datum=WGS84", zone=NA, flip_around_dateline=FALSE, silent=FALSE ){ # Original projection origCRS = sp::CRS(origargs) diff --git a/R/zzz.R b/R/zzz.R index adf2035..a2824fb 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,6 +1,6 @@ -#.onLoad <- function(libname, pkgname) { -#} +# .onLoad <- function(libname, pkgname) { +# } .onAttach <- function(libname, pkgname) { if( getOption("repos")["CRAN"] == "@CRAN@" ){ @@ -46,4 +46,27 @@ Rotate_Fn = function( ... ){ FishStatsUtils::rotate_factors( ... ) } +#' Copy of FishStatsUtils::make_mesh +#' +#' Included for continuity when using old scripts +#' +#' Please use \code{?FishStatsUtils::make_mesh} to see list of arguments and usage +#' @export +Calc_Anisotropic_Mesh = function( ... ){ + .Deprecated( new="FishStatsUtils::make_mesh" ) + FishStatsUtils::make_mesh( ... ) +} + +#' Copy of FishStatsUtils::make_kmeans +#' +#' Included for continuity when using old scripts +#' +#' Please use \code{?FishStatsUtils::make_kmeans} to see list of arguments and usage +#' @export +Calc_Kmeans = function( ... ){ + .Deprecated( new="FishStatsUtils::make_kmeans" ) + FishStatsUtils::make_kmeans( ... ) +} + + diff --git a/man/Calc_Anisotropic_Mesh.Rd b/man/Calc_Anisotropic_Mesh.Rd index 9adecfe..90ad1fc 100644 --- a/man/Calc_Anisotropic_Mesh.Rd +++ b/man/Calc_Anisotropic_Mesh.Rd @@ -1,24 +1,30 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Calc_Anisotropic_Mesh.R +% Please edit documentation in R/Calc_Anisotropic_Mesh.R, R/zzz.R \name{Calc_Anisotropic_Mesh} \alias{Calc_Anisotropic_Mesh} \title{Make mesh for distances among points} \usage{ -Calc_Anisotropic_Mesh(loc_x, loc_g, loc_i, Method, Extrapolation_List, - anisotropic_mesh = NULL, fine_scale = FALSE, ...) +Calc_Anisotropic_Mesh(...) + +Calc_Anisotropic_Mesh(...) } \arguments{ +\item{...}{Arguments passed to \code{INLA::inla.mesh.create}} + \item{loc_x}{location (eastings and northings in kilometers, UTM) for each sample or knot} \item{Method}{spatial method determines ("Mesh" and "Grid" give} \item{anisotropic_mesh}{OPTIONAL, anisotropic mesh (if missing, its recalculated from loc_x)} - -\item{...}{Arguments passed to \code{INLA::inla.mesh.create}} } \value{ Tagged list containing distance metrics } \description{ \code{Calc_Anistropic_Mesh} builds a tagged list representing distances for isotropic or geometric anisotropic triangulated mesh + +Included for continuity when using old scripts +} +\details{ +Please use \code{?FishStatsUtils::make_mesh} to see list of arguments and usage } diff --git a/man/Calc_Kmeans.Rd b/man/Calc_Kmeans.Rd index b5c77c1..6bf73fd 100644 --- a/man/Calc_Kmeans.Rd +++ b/man/Calc_Kmeans.Rd @@ -1,35 +1,14 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Calc_Kmeans.R +% Please edit documentation in R/zzz.R \name{Calc_Kmeans} \alias{Calc_Kmeans} -\title{Calculate location for knots approximating spatial variation} +\title{Copy of FishStatsUtils::make_kmeans} \usage{ -Calc_Kmeans(n_x, loc_orig, nstart = 100, randomseed = 1, - iter.max = 1000, DirPath = paste0(getwd(), "/"), - Save_Results = TRUE) -} -\arguments{ -\item{n_x}{the number of knots to select} - -\item{loc_orig}{a matrix with two columns where each row gives the 2-dimensional coordinates to be approximated} - -\item{nstart}{the number of times that the k-means algorithm is run while searching for the best solution (default=100)} - -\item{randomseed}{a random number seed} - -\item{iter.max}{the number of iterations used per k-means algorithm (default=1000)} - -\item{DirPath}{a directory where the algorithm looks for a previously-saved output (default is working directory)} - -\item{Save_Results}{a boolean stating whether to save the output (Default=TRUE)} -} -\value{ -Tagged list containing outputs -\describe{ - \item{centers}{a matrix with 2 columns and n_x rows} - \item{cluster}{A vector with length \code{nrow(loc_orig)} specifying which row of \code{centers} corresponds to each row of loc_orig} -} +Calc_Kmeans(...) } \description{ -\code{Calc_Kmeans} determines the location for a set of knots for approximating spatial variation +Included for continuity when using old scripts +} +\details{ +Please use \code{?FishStatsUtils::make_kmeans} to see list of arguments and usage } diff --git a/man/Calc_Polygon_Areas_and_Polygons_Fn.Rd b/man/Calc_Polygon_Areas_and_Polygons_Fn.Rd index 099faec..fc383c7 100644 --- a/man/Calc_Polygon_Areas_and_Polygons_Fn.Rd +++ b/man/Calc_Polygon_Areas_and_Polygons_Fn.Rd @@ -4,8 +4,12 @@ \alias{Calc_Polygon_Areas_and_Polygons_Fn} \title{Calculate areas and covariates for a given triangulated mesh} \usage{ -Calc_Polygon_Areas_and_Polygons_Fn(loc_x, Data_Extrap, - Covariates = "none", a_el = NULL) +Calc_Polygon_Areas_and_Polygons_Fn( + loc_x, + Data_Extrap, + Covariates = "none", + a_el = NULL +) } \arguments{ \item{loc_x}{location (eastings and northings in kilometers, UTM) for each sample or knot} diff --git a/man/Coherence.Rd b/man/Coherence.Rd index 5a75db1..65892d3 100644 --- a/man/Coherence.Rd +++ b/man/Coherence.Rd @@ -4,8 +4,12 @@ \alias{Coherence} \title{Calculate stability metrics} \usage{ -Coherence(Report, Data, covhat = NULL, yearbounds_zz = matrix(c(1, - Data$n_t), nrow = 1)) +Coherence( + Report, + Data, + covhat = NULL, + yearbounds_zz = matrix(c(1, Data$n_t), nrow = 1) +) } \arguments{ \item{Report}{output report, e.g., from \code{Report <- obj$report()}} diff --git a/man/Crossvalidate_Fn.Rd b/man/Crossvalidate_Fn.Rd index d72bc11..613480a 100644 --- a/man/Crossvalidate_Fn.Rd +++ b/man/Crossvalidate_Fn.Rd @@ -4,8 +4,16 @@ \alias{Crossvalidate_Fn} \title{K-fold crossvalidation} \usage{ -Crossvalidate_Fn(record_dir, parhat, original_data, group_i = NULL, - kfold = 10, newtonsteps = 1, skip_finished = FALSE, ...) +Crossvalidate_Fn( + record_dir, + parhat, + original_data, + group_i = NULL, + kfold = 10, + newtonsteps = 1, + skip_finished = FALSE, + ... +) } \arguments{ \item{record_dir, }{a directory with writing privilege where the runs are stored} @@ -23,8 +31,6 @@ Crossvalidate_Fn(record_dir, parhat, original_data, group_i = NULL, \item{skip_finished}{boolean specifying whether to rerun (skip_finished==FALSE) or skip (skip_finished==TRUE) previously completed runs (Default=FALSE)} \item{...}{Additional arguments to pass to \code{VAST::Build_TMB_Fn}} - -\item{skip_finished}{boolean specifying whether to rerun (skip_finished==FALSE) or skip (skip_finished==TRUE) previously completed runs (Default=FALSE)} } \value{ Results a matrix with total predictive negative log-likelihood for each crossvalidation partition, and number of crossvalidation samples for that partition diff --git a/man/FishStatsUtils.Rd b/man/FishStatsUtils.Rd index c1b7686..61d1edf 100644 --- a/man/FishStatsUtils.Rd +++ b/man/FishStatsUtils.Rd @@ -3,7 +3,6 @@ \docType{package} \name{FishStatsUtils} \alias{FishStatsUtils} -\alias{FishStatsUtils-package} \title{FishStatsUtils} \description{ FishStatsUtils provides utilities (shared code and data) for spatio-temporal models (VAST, SpatialDeltaGLMM, ...) diff --git a/man/Geostat_Sim.Rd b/man/Geostat_Sim.Rd index e6c2fcb..a65d6aa 100644 --- a/man/Geostat_Sim.Rd +++ b/man/Geostat_Sim.Rd @@ -4,9 +4,14 @@ \alias{Geostat_Sim} \title{Simulate data} \usage{ -Geostat_Sim(Sim_Settings, Extrapolation_List, Data_Geostat = NULL, - MakePlot = FALSE, DateFile = paste0(getwd(), "/"), - standardize_fields = FALSE) +Geostat_Sim( + Sim_Settings, + Extrapolation_List, + Data_Geostat = NULL, + MakePlot = FALSE, + DateFile = paste0(getwd(), "/"), + standardize_fields = FALSE +) } \arguments{ \item{Sim_Settings}{an optional tagged list for specifying input parameters (see function code to determine settings)} diff --git a/man/PlotMap_Fn.Rd b/man/PlotMap_Fn.Rd index 7f391c2..2fd29f0 100644 --- a/man/PlotMap_Fn.Rd +++ b/man/PlotMap_Fn.Rd @@ -4,14 +4,35 @@ \alias{PlotMap_Fn} \title{Plot maps with areal results} \usage{ -PlotMap_Fn(MappingDetails, Mat, PlotDF, MapSizeRatio = c(`Width(in)` = 4, - `Height(in)` = 4), Xlim, Ylim, FileName = paste0(getwd(), "/"), - Year_Set, Rescale = FALSE, Rotate = 0, Format = "png", Res = 200, - zone = NA, Cex = 0.01, textmargin = "", add = FALSE, pch = 15, - outermargintext = c("Eastings", "Northings"), zlim = NULL, - Col = NULL, Legend = list(use = FALSE, x = c(10, 30), y = c(10, 30)), - mfrow = c(1, 1), plot_legend_fig = TRUE, land_color = "grey", - ignore.na = FALSE, map_style = "rescale", ...) +PlotMap_Fn( + MappingDetails, + Mat, + PlotDF, + MapSizeRatio = c(`Width(in)` = 4, `Height(in)` = 4), + Xlim, + Ylim, + FileName = paste0(getwd(), "/"), + Year_Set, + Rescale = FALSE, + Rotate = 0, + Format = "png", + Res = 200, + zone = NA, + Cex = 0.01, + textmargin = "", + add = FALSE, + pch = 15, + outermargintext = c("Eastings", "Northings"), + zlim = NULL, + Col = NULL, + Legend = list(use = FALSE, x = c(10, 30), y = c(10, 30)), + mfrow = c(1, 1), + plot_legend_fig = TRUE, + land_color = "grey", + ignore.na = FALSE, + map_style = "rescale", + ... +) } \arguments{ \item{MapSizeRatio}{Default size for each panel} diff --git a/man/Rerun_Fn.Rd b/man/Rerun_Fn.Rd index 1c0537c..c5a1124 100644 --- a/man/Rerun_Fn.Rd +++ b/man/Rerun_Fn.Rd @@ -4,11 +4,19 @@ \alias{Rerun_Fn} \title{Explore counter-factual scenario} \usage{ -Rerun_Fn(parhat0, turnoff_pars, loc_x, +Rerun_Fn( + parhat0, + turnoff_pars, + loc_x, cov_to_turnoff = 1:dim(parhat0[["gamma2_ctp"]])[3], - calculate_COG = TRUE, figname = NULL, Map = "generate", - MapDetails_List = NULL, year_set = 1:ncol(parhat0[["beta1_ct"]]), - c_set = 1:nrow(parhat0[["beta1_ct"]]), ...) + calculate_COG = TRUE, + figname = NULL, + Map = "generate", + MapDetails_List = NULL, + year_set = 1:ncol(parhat0[["beta1_ct"]]), + c_set = 1:nrow(parhat0[["beta1_ct"]]), + ... +) } \arguments{ \item{parhat0}{parameter values to use by default (presumably maximum-likelihood estimates)} diff --git a/man/calc_coherence.Rd b/man/calc_coherence.Rd index 6dc3275..8849c0a 100644 --- a/man/calc_coherence.Rd +++ b/man/calc_coherence.Rd @@ -4,8 +4,12 @@ \alias{calc_coherence} \title{Calculate stability metrics} \usage{ -calc_coherence(Report, Data, covhat = NULL, yearbounds_zz = matrix(c(1, - Data$n_t), nrow = 1)) +calc_coherence( + Report, + Data, + covhat = NULL, + yearbounds_zz = matrix(c(1, Data$n_t), nrow = 1) +) } \arguments{ \item{Report}{output report, e.g., from \code{Report <- obj$report()}} diff --git a/man/calc_synchrony.Rd b/man/calc_synchrony.Rd index d3d6f20..a987461 100644 --- a/man/calc_synchrony.Rd +++ b/man/calc_synchrony.Rd @@ -4,8 +4,7 @@ \alias{calc_synchrony} \title{Calculate synchrony among species or locations} \usage{ -calc_synchrony(Report, Data, yearbounds_zz = matrix(c(1, Data$n_t), nrow - = 1)) +calc_synchrony(Report, Data, yearbounds_zz = matrix(c(1, Data$n_t), nrow = 1)) } \arguments{ \item{Report}{output report, e.g., from \code{Report <- obj$report()}} diff --git a/man/calculate_proportion.Rd b/man/calculate_proportion.Rd index 04b98f5..e102651 100644 --- a/man/calculate_proportion.Rd +++ b/man/calculate_proportion.Rd @@ -4,12 +4,25 @@ \alias{calculate_proportion} \title{Calculate the compositional-expansion} \usage{ -calculate_proportion(TmbData, Index, Expansion_cz = NULL, - Year_Set = NULL, Years2Include = NULL, strata_names = NULL, - category_names = NULL, plot_legend = ifelse(TmbData$n_l > 1, TRUE, - FALSE), DirName = paste0(getwd(), "/"), PlotName = "Proportion.png", - PlotName2 = "Average.png", interval_width = 1, width = 6, - height = 6, xlab = "Category", ylab = "Proportion", ...) +calculate_proportion( + TmbData, + Index, + Expansion_cz = NULL, + Year_Set = NULL, + Years2Include = NULL, + strata_names = NULL, + category_names = NULL, + plot_legend = ifelse(TmbData$n_l > 1, TRUE, FALSE), + DirName = paste0(getwd(), "/"), + PlotName = "Proportion.png", + PlotName2 = "Average.png", + interval_width = 1, + width = 6, + height = 6, + xlab = "Category", + ylab = "Proportion", + ... +) } \arguments{ \item{TmbData}{Formatted data inputs, from `VAST::Data_Fn(...)`} diff --git a/man/combine_lists.Rd b/man/combine_lists.Rd index 4d70611..fa1dd6a 100644 --- a/man/combine_lists.Rd +++ b/man/combine_lists.Rd @@ -4,8 +4,7 @@ \alias{combine_lists} \title{Combine lists} \usage{ -combine_lists(default, input, args_to_use = NULL, - use_partial_matching = FALSE) +combine_lists(default, input, args_to_use = NULL, use_partial_matching = FALSE) } \arguments{ \item{default}{default values for combined list} diff --git a/man/convert_shapefile.Rd b/man/convert_shapefile.Rd index 425533e..171dfa9 100644 --- a/man/convert_shapefile.Rd +++ b/man/convert_shapefile.Rd @@ -4,18 +4,29 @@ \alias{convert_shapefile} \title{Convert shapefile to extrapolation-grid} \usage{ -convert_shapefile(file_path, projargs = NULL, grid_dim_km = c(2, 2), - projargs_for_shapefile = NULL, make_plots = FALSE, quiet = TRUE, - area_tolerance = 0.05, ...) +convert_shapefile( + file_path, + projargs = NULL, + grid_dim_km = c(2, 2), + projargs_for_shapefile = NULL, + make_plots = FALSE, + quiet = TRUE, + area_tolerance = 0.05, + ... +) } \arguments{ \item{file_path}{path for shapefile on harddrive} \item{projargs}{A character string of projection arguments used to project the shapefile prior to construction an extrapolation-grid; the arguments must be entered exactly as in the PROJ.4 documentation; Default \code{projargs=NULL} uses UTM and infers the UTM zone based on the centroid of the shapefile.} +\item{grid_dim_km}{numeric-vector with length two, giving the distance in km between cells in the automatically generated extrapolation grid; only used if \code{Region="other"}} + \item{projargs_for_shapefile}{projection-arguments (e.g., as parsed by \code{sp::CRS}), that are used when reading shapefile and overriding any projection arguments already saved there; Default \code{projargs_for_shapefile=NULL} uses the projection arguments available in the shapefile} \item{make_plots}{Boolean indicating whether to visualize inputs and outputs as maps} + +\item{...}{other objects passed for individual regions (see example script)} } \value{ extrapolation-grid diff --git a/man/fit_model.Rd b/man/fit_model.Rd index ea099f4..9398afe 100644 --- a/man/fit_model.Rd +++ b/man/fit_model.Rd @@ -4,12 +4,33 @@ \alias{fit_model} \title{Fit VAST to data} \usage{ -fit_model(settings, Lat_i, Lon_i, t_i, b_i, a_i, c_iz = rep(0, - length(b_i)), v_i = rep(0, length(b_i)), - working_dir = paste0(getwd(), "/"), Xconfig_zcp = NULL, - covariate_data, formula = ~0, Q_ik = NULL, newtonsteps = 1, - silent = TRUE, build_model = TRUE, run_model = TRUE, - test_fit = TRUE, ...) +fit_model( + settings, + Lat_i, + Lon_i, + t_i, + b_i, + a_i, + c_iz = rep(0, length(b_i)), + v_i = rep(0, length(b_i)), + working_dir = paste0(getwd(), "/"), + X1config_cp = NULL, + X2config_cp = NULL, + covariate_data, + X1_formula = ~0, + X2_formula = ~0, + Q1config_k = NULL, + Q2config_k = NULL, + catchability_data, + Q1_formula = ~0, + Q2_formula = ~0, + newtonsteps = 1, + silent = TRUE, + build_model = TRUE, + run_model = TRUE, + test_fit = TRUE, + ... +) } \arguments{ \item{settings}{Output from \code{\link{make_settings}}} @@ -18,27 +39,47 @@ fit_model(settings, Lat_i, Lon_i, t_i, b_i, a_i, c_iz = rep(0, \item{Lon_i}{Longitude for each sample} -\item{b_i}{Sampled biomass for each observation i} +\item{t_i}{Time for each observation i} -\item{a_i}{Sampled area for each observation i} +\item{b_i}{Sampled value (biomass, counts, etc.) for each observation i} + +\item{a_i}{Sampled area for each observation i; use \code{a_i=1} for observations without a natural area measurement, while +noting that resulting densities no longer have interpretable units in that case)} \item{c_iz}{Category (e.g., species, length-bin) for each observation i} -\item{v_i}{sampling category (e.g., vessel or tow) associated with overdispersed variation for each observation i (by default \code{v_i=0} for all samples, which will not affect things given the default values for \code{OverdispersionConfig})} +\item{v_i}{sampling category (e.g., vessel or tow) associated with overdispersed variation for each observation i +(by default \code{v_i=0} for all samples, which will not affect things given the default values for \code{OverdispersionConfig})} -\item{Xconfig_zcp}{OPTIONAL, 3D array of settings for each dynamic density covariate, where the first dimension corresponds to 1st or 2nd linear predictors, second dimension corresponds to model category, and third dimension corresponds to each density covariate +\item{X1config_cp}{matrix of settings for each density covariate for the 1st lienar predictor, + where the row corresponds to model category, and column corresponds to each density covariate \describe{ - \item{Xconfig_zcp[z,c,p]=0}{\code{X_itp[,,p]} has no effect on linear predictor z for category c} - \item{Xconfig_zcp[z,c,p]=1}{\code{X_itp[,,p]} has a linear effect on linear predictor z for category c} - \item{Xconfig_zcp[z,c,p]=2}{\code{X_itp[,,p]} has a spatially varying, zero-centered linear effect on linear predictor z for category c} - \item{Xconfig_zcp[z,c,p]=3}{\code{X_itp[,,p]} has a spatially varying linear effect on linear predictor z for category c} + \item{X1config_cp[c,p]=0}{\code{X_ip[,p]} has no effect on the 1st linear predictor for category c} + \item{X1config_cp[c,p]=1}{\code{X_ip[,p]} has a linear effect on 1st linear predictor for category c} + \item{X1config_cp[c,p]=2}{\code{X_ip[,p]} has a spatially varying, zero-centered linear effect on 1st linear predictor for category c} + \item{X1config_cp[c,p]=3}{\code{X_ip[,p]} has a spatially varying linear effect on 1st linear predictor for category c} + \item{X1config_cp[c,p]=-1}{\code{X1config_cp[c,p]=-1} is the same as \code{X1config_cp[c,p]=2}, but without including the log-likelihood term; this is useful in special cases when carefully mirroring spatially varying coefficients, e.g., to use cohort effects in a age-structured spatio-temporal model} }} +\item{X2config_cp}{Same as argument \code{X1config_cp} but for 2nd linear predictor} + \item{covariate_data}{data frame of covariate values with columns \code{Lat}, \code{Lon}, and \code{Year}, and other columns matching names in \code{formula}; \code{Year=NA} can be used for covariates that do not change among years (e.g., depth)} -\item{formula}{an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. Similar specification to \code{\link{stats::lm}}} +\item{X1_formula}{right-sided formula affecting the 1st linear predictor for density, e.g., +\code{X1_formula=~BOT_DEPTH+BOT_DEPTH^2} for a quadratic effect of variable \code{BOT_DEPTH}} + +\item{X2_formula}{same as \code{X1_formula} but affecting the 2nd linear predictor for density} + +\item{Q1config_k}{Same as argument \code{X1config_cp} but affecting affecting the 1st linear predictor for catchability, +and note that it is a vector (instead of matrix) given the more flexible form for catchability covariates} + +\item{Q2config_k}{Same as argument \code{Q1config_cp} but affecting affecting the 2nd linear predictor for catchability} + +\item{catchability_data}{data-frame of covariates for use when specifying \code{Q1_formula} and \code{Q2_formula}} + +\item{Q1_formula}{same as \code{X1_formula} but affecting the 1st linear predictor for catchability (a.k.a. detectability)} -\item{Q_ik}{matrix of catchability covariates (e.g., measured variables affecting catch rates but not caused by variation in species density) for each observation i} +\item{Q2_formula}{same as \code{X1_formula} but affecting the 2nd linear predictor for catchability (a.k.a. detectability)} \item{newtonsteps}{number of extra newton steps to take after optimization (alternative to \code{loopnum})} @@ -107,14 +148,20 @@ library(VAST) example = load_example( data_set="EBS_pollock" ) # Make settings -settings = make_settings( n_x=50, Region=example$Region, purpose="index", - strata.limits=example$strata.limits ) +settings = make_settings( n_x=50, + Region=example$Region, + purpose="index", + strata.limits=example$strata.limits ) # Run model -fit = fit_model( "settings"=settings, "Lat_i"=example$sampling_data[,'Lat'], - "Lon_i"=example$sampling_data[,'Lon'], "t_i"=example$sampling_data[,'Year'], - "c_i"=rep(0,nrow(example$sampling_data)), "b_i"=example$sampling_data[,'Catch_KG'], - "a_i"=example$sampling_data[,'AreaSwept_km2'], "v_i"=example$sampling_data[,'Vessel'] ) +fit = fit_model( "settings"=settings, + "Lat_i"=example$sampling_data[,'Lat'], + "Lon_i"=example$sampling_data[,'Lon'], + "t_i"=example$sampling_data[,'Year'], + "c_i"=rep(0,nrow(example$sampling_data)), + "b_i"=example$sampling_data[,'Catch_KG'], + "a_i"=example$sampling_data[,'AreaSwept_km2'], + "v_i"=example$sampling_data[,'Vessel'] ) # Plot results plot_results( settings=settings, fit=fit ) @@ -122,11 +169,12 @@ plot_results( settings=settings, fit=fit ) } \seealso{ -\code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots +\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 \code{\link{summary.fit_model}} for methods to summarize output, including obtain a dataframe of estimated densities -Other wrapper functions: \code{\link{make_settings}}, - \code{\link{plot_results}} +Other wrapper functions: +\code{\link{make_settings}()}, +\code{\link{plot_results}()} } \concept{wrapper functions} diff --git a/man/format_covariates.Rd b/man/format_covariates.Rd index 18d54a8..6c2e0cb 100644 --- a/man/format_covariates.Rd +++ b/man/format_covariates.Rd @@ -4,9 +4,17 @@ \alias{format_covariates} \title{Format habitat covariate matrix} \usage{ -format_covariates(Lat_e, Lon_e, t_e, Cov_ep, Extrapolation_List, - Spatial_List, FUN = mean, Year_Set = min(t_e):max(t_e), - na.omit = "error") +format_covariates( + Lat_e, + Lon_e, + t_e, + Cov_ep, + Extrapolation_List, + Spatial_List, + FUN = mean, + Year_Set = min(t_e):max(t_e), + na.omit = "error" +) } \arguments{ \item{Lat_e, }{Latitude for covariate sample e} diff --git a/man/make_covariates.Rd b/man/make_covariates.Rd index 128333e..87d549c 100644 --- a/man/make_covariates.Rd +++ b/man/make_covariates.Rd @@ -4,8 +4,13 @@ \alias{make_covariates} \title{Format habitat covariate matrix} \usage{ -make_covariates(formula, covariate_data, Year_i, spatial_list, - extrapolation_list) +make_covariates( + formula, + covariate_data, + Year_i, + spatial_list, + extrapolation_list +) } \arguments{ \item{formula}{an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. Similar specification to \code{\link{stats::lm}}} diff --git a/man/make_extrapolation_info.Rd b/man/make_extrapolation_info.Rd index c0dbc16..1b45572 100644 --- a/man/make_extrapolation_info.Rd +++ b/man/make_extrapolation_info.Rd @@ -4,17 +4,31 @@ \alias{make_extrapolation_info} \title{Build extrapolation grid} \usage{ -make_extrapolation_info(Region, projargs = NA, zone = NA, +make_extrapolation_info( + Region, + projargs = NA, + zone = NA, strata.limits = data.frame(STRATA = "All_areas"), - create_strata_per_region = FALSE, max_cells = NULL, - input_grid = NULL, observations_LL = NULL, grid_dim_km = c(2, 2), - maximum_distance_from_sample = NULL, grid_in_UTM = TRUE, - grid_dim_LL = c(0.1, 0.1), region = c("south_coast", "west_coast"), + create_strata_per_region = FALSE, + max_cells = NULL, + input_grid = NULL, + observations_LL = NULL, + grid_dim_km = c(2, 2), + maximum_distance_from_sample = NULL, + grid_in_UTM = TRUE, + grid_dim_LL = c(0.1, 0.1), + region = c("south_coast", "west_coast"), strata_to_use = c("SOG", "WCVI", "QCS", "HS", "WCHG"), - epu_to_use = c("All", "Georges_Bank", "Mid_Atlantic_Bight", - "Scotian_Shelf", "Gulf_of_Maine", "Other")[1], survey = "Chatham_rise", - surveyname = "propInWCGBTS", flip_around_dateline, nstart = 100, - area_tolerance = 0.05, ...) + epu_to_use = c("All", "Georges_Bank", "Mid_Atlantic_Bight", "Scotian_Shelf", + "Gulf_of_Maine", "Other")[1], + survey = "Chatham_rise", + surveyname = "propInWCGBTS", + flip_around_dateline, + nstart = 100, + area_tolerance = 0.05, + backwards_compatible_kmeans = FALSE, + ... +) } \arguments{ \item{Region}{a character vector, where each element that is matched against potential values to determine the region for the extrapolation grid. Current options are "california_current", "west_coast_hook_and_line", "british_columbia", "eastern_bering_sea", "northern_bering_sea", "bering_sea_slope", "st_matthews_island", "aleutian_islands", "gulf_of_alaska", "northwest_atlantic", "south_africa", "gulf_of_st_lawrence", "new_zealand", "habcam", "gulf_of_mexico", "ATL-IBTS-Q1", "ATL-IBTS-Q4", "BITS", "BTS", "BTS-VIIA", "EVHOE", "IE-IGFS", "NIGFS", "NS_IBTS", "PT-IBTS", "SP-ARSA", "SP-NORTH", "SP-PORC", "stream_network", "user", "other", or the absolute path and file name for a GIS shapefile} @@ -55,6 +69,10 @@ make_extrapolation_info(Region, projargs = NA, zone = NA, \item{nstart}{the number of times that the k-means algorithm is run while searching for the best solution (default=100)} +\item{backwards_compatible_kmeans}{a boolean stating how to deal with changes in the kmeans algorithm implemented in R version 3.6.0, +where \code{backwards_compatible_kmeans==TRUE} modifies the default algorithm to maintain backwards compatibility, and +where \code{backwards_compatible_kmeans==FALSE} breaks backwards compatibility between R versions prior to and after R 3.6.0.} + \item{...}{other objects passed for individual regions (see example script)} } \value{ diff --git a/man/make_kmeans.Rd b/man/make_kmeans.Rd new file mode 100644 index 0000000..9b3bca1 --- /dev/null +++ b/man/make_kmeans.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/make_kmeans.R +\name{make_kmeans} +\alias{make_kmeans} +\title{Calculate location for knots approximating spatial variation} +\usage{ +make_kmeans( + n_x, + loc_orig, + nstart = 100, + randomseed = 1, + iter.max = 1000, + DirPath = paste0(getwd(), "/"), + Save_Results = TRUE, + backwards_compatible_kmeans = FALSE +) +} +\arguments{ +\item{n_x}{the number of knots to select} + +\item{loc_orig}{a matrix with two columns where each row gives the 2-dimensional coordinates to be approximated} + +\item{nstart}{the number of times that the k-means algorithm is run while searching for the best solution (default=100)} + +\item{randomseed}{a random number seed} + +\item{iter.max}{the number of iterations used per k-means algorithm (default=1000)} + +\item{DirPath}{a directory where the algorithm looks for a previously-saved output (default is working directory)} + +\item{Save_Results}{a boolean stating whether to save the output (Default=TRUE)} + +\item{backwards_compatible_kmeans}{a boolean stating how to deal with changes in the kmeans algorithm implemented in R version 3.6.0, +where \code{backwards_compatible_kmeans==TRUE} modifies the default algorithm to maintain backwards compatibility, and +where \code{backwards_compatible_kmeans==FALSE} breaks backwards compatibility between R versions prior to and after R 3.6.0.} +} +\value{ +Tagged list containing outputs +\describe{ + \item{centers}{a matrix with 2 columns and n_x rows} + \item{cluster}{A vector with length \code{nrow(loc_orig)} specifying which row of \code{centers} corresponds to each row of loc_orig} +} +} +\description{ +\code{make_kmeans} determines the location for a set of knots for approximating spatial variation +} diff --git a/man/make_map_info.Rd b/man/make_map_info.Rd index 34b6171..f8d3c45 100644 --- a/man/make_map_info.Rd +++ b/man/make_map_info.Rd @@ -4,9 +4,14 @@ \alias{make_map_info} \title{Get information for adequate mapping for regions} \usage{ -make_map_info(Region, Extrapolation_List, spatial_list = NULL, +make_map_info( + Region, + Extrapolation_List, + spatial_list = NULL, NN_Extrap = spatial_list$PolygonList$NN_Extrap, - fine_scale = spatial_list$fine_scale, Include) + fine_scale = spatial_list$fine_scale, + Include +) } \description{ Function provides information to be used when plotting results diff --git a/man/make_mesh.Rd b/man/make_mesh.Rd new file mode 100644 index 0000000..1511154 --- /dev/null +++ b/man/make_mesh.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/make_mesh.R +\name{make_mesh} +\alias{make_mesh} +\title{Make mesh for distances among points} +\usage{ +make_mesh( + loc_x, + loc_g, + loc_i, + Method, + Extrapolation_List, + anisotropic_mesh = NULL, + fine_scale = FALSE, + ... +) +} +\arguments{ +\item{loc_x}{location (eastings and northings in kilometers, UTM) for each sample or knot} + +\item{Method}{spatial method determines ("Mesh" and "Grid" give} + +\item{anisotropic_mesh}{OPTIONAL, anisotropic mesh (if missing, its recalculated from loc_x)} + +\item{...}{Arguments passed to \code{INLA::inla.mesh.create}} +} +\value{ +Tagged list containing distance metrics +} +\description{ +\code{make_mesh} builds a tagged list representing distances for isotropic or geometric anisotropic triangulated mesh +} diff --git a/man/make_settings.Rd b/man/make_settings.Rd index 9908ed3..304391a 100644 --- a/man/make_settings.Rd +++ b/man/make_settings.Rd @@ -4,12 +4,28 @@ \alias{make_settings} \title{Make list of settings} \usage{ -make_settings(n_x, Region, purpose = "index", fine_scale = TRUE, - strata.limits = data.frame(STRATA = "All_areas"), zone = NA, - FieldConfig, RhoConfig, OverdispersionConfig, ObsModel, bias.correct, - Options, use_anisotropy, vars_to_correct, Version, - treat_nonencounter_as_zero, n_categories, VamConfig, max_cells, - knot_method) +make_settings( + n_x, + Region, + purpose = "index", + fine_scale = TRUE, + strata.limits = data.frame(STRATA = "All_areas"), + zone = NA, + FieldConfig, + RhoConfig, + OverdispersionConfig, + ObsModel, + bias.correct, + Options, + use_anisotropy, + vars_to_correct, + Version, + treat_nonencounter_as_zero, + n_categories, + VamConfig, + max_cells, + knot_method +) } \arguments{ \item{n_x}{the number of vertices in the SPDE mesh (determines the spatial resolution when Method="Mesh")} @@ -24,7 +40,10 @@ make_settings(n_x, Region, purpose = "index", fine_scale = TRUE, \item{zone}{UTM zone used for projecting Lat-Lon to km distances; use \code{zone=NA} by default to automatically detect UTM zone from the location of extrapolation-grid samples} -\item{FieldConfig}{a vector of format c("Omega1"=0, "Epsilon1"=10, "Omega2"="AR1", "Epsilon2"=10), where Omega refers to spatial variation, Epsilon refers to spatio-temporal variation, Omega1 refers to variation in encounter probability, and Omega2 refers to variation in positive catch rates, where 0 is off, "AR1" is an AR1 process, and >0 is the number of elements in a factor-analysis covariance} +\item{FieldConfig}{a vector of format c("Omega1"=0, "Epsilon1"=10, "Omega2"="AR1", "Epsilon2"=10), where Omega refers to spatial variation, +Epsilon refers to spatio-temporal variation, Omega1 refers to variation in encounter probability, and +Omega2 refers to variation in positive catch rates, where 0 is off, "AR1" is an AR1 process, +and >0 is the number of elements in a factor-analysis covariance} \item{RhoConfig}{vector of form c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsilon2"=0) specifying whether either intercepts (Beta1 and Beta2) or spatio-temporal variation (Epsilon1 and Epsilon2) is structured among time intervals (0: each year as fixed effect; 1: each year as random following IID distribution; 2: each year as random following a random walk; 3: constant among years as fixed effect; 4: each year as random following AR1 process); If missing, assumed to be zero for each element} @@ -63,9 +82,10 @@ Tagged list containing default settings for a given purpose, use \code{names} on This function assembles a default set of user-decisions for a specified modelling purpose. The default settings are guessed based on generic guidance, and should be carefully reviewed for real-world purposes. If the user supplies values for individual settings e.g. \code{FieldConfig}, then these values override the defaults that are provided by interpreting \code{purpose} } \seealso{ -\code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots +\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 -Other wrapper functions: \code{\link{fit_model}}, - \code{\link{plot_results}} +Other wrapper functions: +\code{\link{fit_model}()}, +\code{\link{plot_results}()} } \concept{wrapper functions} diff --git a/man/make_spatial_info.Rd b/man/make_spatial_info.Rd index a10cff8..54c1830 100644 --- a/man/make_spatial_info.Rd +++ b/man/make_spatial_info.Rd @@ -4,12 +4,29 @@ \alias{make_spatial_info} \title{Build objects related to spatial information} \usage{ -make_spatial_info(n_x, Lon_i, Lat_i, Extrapolation_List, - knot_method = NULL, Method = "Mesh", grid_size_km = 50, - grid_size_LL = 1, fine_scale = FALSE, Network_sz_LL = NULL, - iter.max = 1000, randomseed = 1, nstart = 100, - DirPath = paste0(getwd(), "/"), Save_Results = FALSE, LON_intensity, - LAT_intensity, ...) +make_spatial_info( + n_x, + Lon_i, + Lat_i, + Extrapolation_List, + knot_method = NULL, + Method = "Mesh", + anisotropic_mesh = NULL, + Kmeans = NULL, + grid_size_km = 50, + grid_size_LL = 1, + fine_scale = FALSE, + Network_sz_LL = NULL, + iter.max = 1000, + randomseed = 1, + nstart = 100, + DirPath = paste0(getwd(), "/"), + Save_Results = FALSE, + LON_intensity, + LAT_intensity, + backwards_compatible_kmeans = FALSE, + ... +) } \arguments{ \item{n_x}{the number of vertices in the SPDE mesh (determines the spatial resolution when Method="Mesh")} @@ -24,6 +41,8 @@ make_spatial_info(n_x, Lon_i, Lat_i, Extrapolation_List, \item{Method}{a character of either "Grid" or "Mesh" where "Grid" is a 2D AR1 process, and "Mesh" is the SPDE method with geometric anisotropy} +\item{anisotropic_mesh}{OPTIONAL, anisotropic mesh (if missing, its recalculated from loc_x)} + \item{grid_size_km}{the distance between grid cells for the 2D AR1 grid (determines spatial resolution when Method="Grid") when not using \code{Method="Spherical_mesh"}} \item{grid_size_LL}{the distance between grid cells for the 2D AR1 grid (determines spatial resolution when Method="Grid") when using \code{Method="Spherical_mesh"}} @@ -42,6 +61,10 @@ make_spatial_info(n_x, Lon_i, Lat_i, Extrapolation_List, \item{Save_Results}{a boolean stating whether to save the output (Default=TRUE)} +\item{backwards_compatible_kmeans}{a boolean stating how to deal with changes in the kmeans algorithm implemented in R version 3.6.0, +where \code{backwards_compatible_kmeans==TRUE} modifies the default algorithm to maintain backwards compatibility, and +where \code{backwards_compatible_kmeans==FALSE} breaks backwards compatibility between R versions prior to and after R 3.6.0.} + \item{...}{additional arguments passed to \code{\link[INLA]{inla.mesh.create}}} } \value{ @@ -51,7 +74,7 @@ Tagged list containing objects for running a VAST model \item{GridList}{A tagged list with inputs related to the 2D AR1 grid} \item{a_xl}{A data frame with areas for each knot and each strattum} \item{loc_UTM}{A data frame with the converted UTM coordinates for each sample} - \item{Kmeans}{Output from \code{Calc_Kmeans} with knots for a triangulated mesh} + \item{Kmeans}{Output from \code{make_kmeans} with knots for a triangulated mesh} \item{knot_i}{The knot associated with each sample} \item{Method}{The Method input (for archival purposes)} \item{loc_x}{The UTM location for each knot} diff --git a/man/plot.make_extrapolation_info.Rd b/man/plot.make_extrapolation_info.Rd index 1ab31d2..3f37fdd 100644 --- a/man/plot.make_extrapolation_info.Rd +++ b/man/plot.make_extrapolation_info.Rd @@ -4,8 +4,7 @@ \alias{plot.make_extrapolation_info} \title{Plot extrapolation-grid} \usage{ -\method{plot}{make_extrapolation_info}(x, cex = 0.01, - land_color = "grey", map_resolution = "medium", ...) +\method{plot}{make_extrapolation_info}(x, cex = 0.01, land_color = "grey", map_resolution = "medium", ...) } \arguments{ \item{x}{Output from \code{\link{make_extrapolation_info}}} diff --git a/man/plot_biomass_index.Rd b/man/plot_biomass_index.Rd index 5d566dc..1ced5b6 100644 --- a/man/plot_biomass_index.Rd +++ b/man/plot_biomass_index.Rd @@ -4,12 +4,26 @@ \alias{plot_biomass_index} \title{Plot index of abundance} \usage{ -plot_biomass_index(TmbData, Sdreport, Year_Set = NULL, - Years2Include = NULL, DirName = paste0(getwd(), "/"), - PlotName = "Index", interval_width = 1, strata_names = NULL, - category_names = NULL, use_biascorr = TRUE, plot_legend = TRUE, - total_area_km2 = NULL, plot_log = FALSE, width = NULL, - height = NULL, create_covariance_table = FALSE, ...) +plot_biomass_index( + TmbData, + Sdreport, + Year_Set = NULL, + Years2Include = NULL, + DirName = paste0(getwd(), "/"), + PlotName = "Index", + interval_width = 1, + strata_names = NULL, + category_names = NULL, + use_biascorr = TRUE, + plot_legend = TRUE, + total_area_km2 = NULL, + plot_log = FALSE, + width = NULL, + height = NULL, + create_covariance_table = FALSE, + Yrange = c(ifelse(plot_log == TRUE, NA, 0), NA), + ... +) } \arguments{ \item{TmbData}{Formatted data inputs, from `VAST::Data_Fn(...)`} diff --git a/man/plot_cov.Rd b/man/plot_cov.Rd index 18934a6..3278e9b 100644 --- a/man/plot_cov.Rd +++ b/man/plot_cov.Rd @@ -4,8 +4,15 @@ \alias{plot_cov} \title{Plot covariance matrix} \usage{ -plot_cov(Cov, zlim = NULL, names = 1:nrow(Cov), names2 = names, - ncolors = 21, digits = 2, ...) +plot_cov( + Cov, + zlim = NULL, + names = 1:nrow(Cov), + names2 = names, + ncolors = 21, + digits = 2, + ... +) } \arguments{ \item{Cov}{matrix (covariance or correlation) used for plotting} diff --git a/man/plot_data.Rd b/man/plot_data.Rd index 98c6fa8..7939a0a 100644 --- a/man/plot_data.Rd +++ b/man/plot_data.Rd @@ -4,13 +4,26 @@ \alias{plot_data} \title{Plot location of sampling data} \usage{ -plot_data(Extrapolation_List, Spatial_List, Data_Geostat = NULL, - Lat_i = Data_Geostat[, "Lat"], Lon_i = Data_Geostat[, "Lon"], - Year_i = Data_Geostat[, "Year"], PlotDir = paste0(getwd(), "/"), - Plot1_name = "Data_and_knots.png", Plot2_name = "Data_by_year.png", - col = "red", cex = 0.01, pch = 19, Year_Set, - projargs = "+proj=longlat", map_resolution = "medium", - land_color = "grey", country = NULL, ...) +plot_data( + Extrapolation_List, + Spatial_List, + Data_Geostat = NULL, + Lat_i = Data_Geostat[, "Lat"], + Lon_i = Data_Geostat[, "Lon"], + Year_i = Data_Geostat[, "Year"], + PlotDir = paste0(getwd(), "/"), + Plot1_name = "Data_and_knots.png", + Plot2_name = "Data_by_year.png", + col = "red", + cex = 0.01, + pch = 19, + Year_Set, + projargs = "+proj=longlat", + map_resolution = "medium", + land_color = "grey", + country = NULL, + ... +) } \arguments{ \item{Extrapolation_List}{Output from \code{Prepare_Extrapolation_Data_Fn}} diff --git a/man/plot_encounter_diagnostic.Rd b/man/plot_encounter_diagnostic.Rd index 86a936e..3e2adb5 100644 --- a/man/plot_encounter_diagnostic.Rd +++ b/man/plot_encounter_diagnostic.Rd @@ -4,9 +4,15 @@ \alias{plot_encounter_diagnostic} \title{Check predicted encounter probability against observed encounter frequency} \usage{ -plot_encounter_diagnostic(Report, Data_Geostat, cutpoints_z = seq(0, 1, - length = 21), interval_width = 1.96, DirName = paste0(getwd(), "/"), - PlotName = "Diag--Encounter_prob.png", ...) +plot_encounter_diagnostic( + Report, + Data_Geostat, + cutpoints_z = seq(0, 1, length = 21), + interval_width = 1.96, + DirName = paste0(getwd(), "/"), + PlotName = "Diag--Encounter_prob.png", + ... +) } \arguments{ \item{Report}{tagged list of outputs from TMB model via \code{Obj$report()}} diff --git a/man/plot_factors.Rd b/man/plot_factors.Rd index 53dc573..52b345e 100644 --- a/man/plot_factors.Rd +++ b/man/plot_factors.Rd @@ -4,11 +4,23 @@ \alias{plot_factors} \title{Plot factor-decomposition of covariance} \usage{ -plot_factors(Report, ParHat, Data, SD = NULL, Year_Set = NULL, - category_names = NULL, RotationMethod = "PCA", - mapdetails_list = NULL, Dim_year = NULL, Dim_species = NULL, - plotdir = paste0(getwd(), "/"), land_color = "grey", zlim = NA, - testcutoff = 1e-04, ...) +plot_factors( + Report, + ParHat, + Data, + SD = NULL, + Year_Set = NULL, + category_names = NULL, + RotationMethod = "PCA", + mapdetails_list = NULL, + Dim_year = NULL, + Dim_species = NULL, + plotdir = paste0(getwd(), "/"), + land_color = "grey", + zlim = NA, + testcutoff = 1e-04, + ... +) } \arguments{ \item{Report}{output report, e.g., from \code{Report <- obj$report()}} diff --git a/man/plot_index.Rd b/man/plot_index.Rd index 02db05c..0cd25e1 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -4,15 +4,34 @@ \alias{plot_index} \title{Plot an index} \usage{ -plot_index(Index_ctl, sd_Index_ctl = array(0, dim(Index_ctl)), - Year_Set = NULL, Years2Include = NULL, strata_names = NULL, - category_names = NULL, scale = "uniform", plot_legend = NULL, - DirName = paste0(getwd(), "/"), PlotName = "Index.png", - interval_width = 1, width = NULL, height = NULL, xlab = "Year", - ylab = "Index", bounds_type = "whiskers", col = NULL, - col_bounds = NULL, Yrange = c(0, NA), type = "b", - plot_lines_args = list(), plot_args = list(), - SampleSize_ctz = NULL, Y2range = c(0, NA), y2lab = "", ...) +plot_index( + Index_ctl, + sd_Index_ctl = array(0, dim(Index_ctl)), + Year_Set = NULL, + Years2Include = NULL, + strata_names = NULL, + category_names = NULL, + scale = "uniform", + plot_legend = NULL, + DirName = paste0(getwd(), "/"), + PlotName = "Index.png", + interval_width = 1, + width = NULL, + height = NULL, + xlab = "Year", + ylab = "Index", + bounds_type = "whiskers", + col = NULL, + col_bounds = NULL, + Yrange = c(0, NA), + type = "b", + plot_lines_args = list(), + plot_args = list(), + SampleSize_ctz = NULL, + Y2range = c(0, NA), + y2lab = "", + ... +) } \arguments{ \item{Index_ctl}{A matrix or array of time-series estimates for multiple categories \code{c}, years \code{t}, and strata \code{l}} diff --git a/man/plot_loadings.Rd b/man/plot_loadings.Rd index 859051c..376f7c1 100644 --- a/man/plot_loadings.Rd +++ b/man/plot_loadings.Rd @@ -4,10 +4,18 @@ \alias{plot_loadings} \title{Plotting loadings matrix} \usage{ -plot_loadings(L_pj, At = 1:nrow(L_pj), whichfactor = 1, - addtitle = TRUE, LabelPosition = "Right", Buffer = c(0, 0.1), - Labels = rownames(L_pj), Cex = 1.2, - legend_text = "Proportion of explained variance", ...) +plot_loadings( + L_pj, + At = 1:nrow(L_pj), + whichfactor = 1, + addtitle = TRUE, + LabelPosition = "Right", + Buffer = c(0, 0.1), + Labels = rownames(L_pj), + Cex = 1.2, + legend_text = "Proportion of explained variance", + ... +) } \arguments{ \item{L_pj}{Loadings matrix for `p` categories and `j` factors} diff --git a/man/plot_maps.Rd b/man/plot_maps.Rd index 22b31f1..9c7cfd1 100644 --- a/man/plot_maps.Rd +++ b/man/plot_maps.Rd @@ -4,12 +4,29 @@ \alias{plot_maps} \title{Plot standard maps} \usage{ -plot_maps(plot_set = 3, Obj = NULL, PlotDF, Sdreport = NULL, - projargs = "+proj=longlat", Panel = "Category", Year_Set = NULL, - Years2Include = NULL, category_names = NULL, quiet = FALSE, - working_dir = paste0(getwd(), "/"), MapSizeRatio, n_cells, - plot_value = "estimate", n_samples = 100, Report, TmbData, - zlim = NULL, country = NULL, sample_fixed = TRUE, ...) +plot_maps( + plot_set = 3, + Obj = NULL, + PlotDF, + Sdreport = NULL, + projargs = "+proj=longlat", + Panel = "Category", + Year_Set = NULL, + Years2Include = NULL, + category_names = NULL, + quiet = FALSE, + working_dir = paste0(getwd(), "/"), + MapSizeRatio, + n_cells, + plot_value = "estimate", + n_samples = 100, + Report, + TmbData, + zlim = NULL, + country = NULL, + sample_fixed = TRUE, + ... +) } \arguments{ \item{plot_set}{integer-vector defining plots to create diff --git a/man/plot_overdispersion.Rd b/man/plot_overdispersion.Rd index 15def59..dd7b2d0 100644 --- a/man/plot_overdispersion.Rd +++ b/man/plot_overdispersion.Rd @@ -4,9 +4,17 @@ \alias{plot_overdispersion} \title{Plot estimated overdispersion} \usage{ -plot_overdispersion(filename1, filename2, Data, ParHat, Report, - SD = NULL, Map = NULL, ControlList1 = list(Width = 8, Height = 4, - Res = 200, Units = "in"), ControlList2 = ControlList1) +plot_overdispersion( + filename1, + filename2, + Data, + ParHat, + Report, + SD = NULL, + Map = NULL, + ControlList1 = list(Width = 8, Height = 4, Res = 200, Units = "in"), + ControlList2 = ControlList1 +) } \arguments{ \item{filename1}{filename (including absolute path) for overdispersion for component #1} diff --git a/man/plot_quantile_diagnostic.Rd b/man/plot_quantile_diagnostic.Rd index 14126bb..b08ed71 100644 --- a/man/plot_quantile_diagnostic.Rd +++ b/man/plot_quantile_diagnostic.Rd @@ -4,11 +4,16 @@ \alias{plot_quantile_diagnostic} \title{Diagnostic QQ function} \usage{ -plot_quantile_diagnostic(TmbData, Report, DateFile = paste0(getwd(), - "/"), save_dir = paste0(DateFile, "/QQ_Fn/"), +plot_quantile_diagnostic( + TmbData, + Report, + DateFile = paste0(getwd(), "/"), + save_dir = paste0(DateFile, "/QQ_Fn/"), FileName_PP = "Posterior_Predictive", FileName_Phist = "Posterior_Predictive-Histogram", - FileName_QQ = "Q-Q_plot", FileName_Qhist = "Q-Q_hist") + FileName_QQ = "Q-Q_plot", + FileName_Qhist = "Q-Q_hist" +) } \arguments{ \item{TmbData}{TMB Model input data list} diff --git a/man/plot_quantile_residuals.Rd b/man/plot_quantile_residuals.Rd index 0035673..22a0360 100644 --- a/man/plot_quantile_residuals.Rd +++ b/man/plot_quantile_residuals.Rd @@ -4,9 +4,14 @@ \alias{plot_quantile_residuals} \title{Plot quantile residuals on map} \usage{ -plot_quantile_residuals(dharmaRes, fit, - file_name = "quantile_residuals_on_map", Year_Set = NULL, - Years2Include = NULL, ...) +plot_quantile_residuals( + dharmaRes, + fit, + file_name = "quantile_residuals_on_map", + Year_Set = NULL, + Years2Include = NULL, + ... +) } \arguments{ \item{...}{arguments passed to \code{FisshStatsUtils::plot_variable}} diff --git a/man/plot_range_edge.Rd b/man/plot_range_edge.Rd index caa9b9a..4cdb13d 100644 --- a/man/plot_range_edge.Rd +++ b/man/plot_range_edge.Rd @@ -4,11 +4,23 @@ \alias{plot_range_edge} \title{Plot shifts in range edges} \usage{ -plot_range_edge(Sdreport, Obj, Year_Set = NULL, Years2Include = NULL, - strata_names = NULL, category_names = NULL, - working_dir = paste0(getwd(), "/"), quantiles = c(0.05, 0.95), - n_samples = 100, interval_width = 1, width = NULL, height = NULL, - calculate_relative_to_average = FALSE, seed = 123456, ...) +plot_range_edge( + Sdreport, + Obj, + Year_Set = NULL, + Years2Include = NULL, + strata_names = NULL, + category_names = NULL, + working_dir = paste0(getwd(), "/"), + quantiles = c(0.05, 0.95), + n_samples = 100, + interval_width = 1, + width = NULL, + height = NULL, + calculate_relative_to_average = FALSE, + seed = 123456, + ... +) } \arguments{ \item{Sdreport}{Standard deviation outputs from TMB model via \code{sdreport(Obj)}} diff --git a/man/plot_range_index.Rd b/man/plot_range_index.Rd index 3003b8e..9e6acdc 100644 --- a/man/plot_range_index.Rd +++ b/man/plot_range_index.Rd @@ -4,14 +4,23 @@ \alias{plot_range_index} \title{Plot shifts in distribution and area occupied} \usage{ -plot_range_index(Sdreport, Report, TmbData, Year_Set = NULL, - Years2Include = NULL, strata_names = NULL, - PlotDir = paste0(getwd(), "/"), FileName_COG = paste0(PlotDir, - "/center_of_gravity.png"), FileName_Area = paste0(PlotDir, - "/Area.png"), FileName_EffArea = paste0(PlotDir, - "/Effective_Area.png"), Znames = rep("", ncol(TmbData$Z_xm)), - use_biascorr = TRUE, category_names = NULL, interval_width = 1, - ...) +plot_range_index( + Sdreport, + Report, + TmbData, + Year_Set = NULL, + Years2Include = NULL, + strata_names = NULL, + PlotDir = paste0(getwd(), "/"), + FileName_COG = paste0(PlotDir, "/center_of_gravity.png"), + FileName_Area = paste0(PlotDir, "/Area.png"), + FileName_EffArea = paste0(PlotDir, "/Effective_Area.png"), + Znames = rep("", ncol(TmbData$Z_xm)), + use_biascorr = TRUE, + category_names = NULL, + interval_width = 1, + ... +) } \arguments{ \item{Sdreport}{TMB output from `TMB::sdreport(Obj)`} diff --git a/man/plot_residuals.Rd b/man/plot_residuals.Rd index 8e85163..f3e72e3 100644 --- a/man/plot_residuals.Rd +++ b/man/plot_residuals.Rd @@ -4,10 +4,21 @@ \alias{plot_residuals} \title{Plot Pearson residuals on map} \usage{ -plot_residuals(Lat_i, Lon_i, TmbData, Report, Q, - projargs = "+proj=longlat", working_dir = paste0(getwd(), "/"), - spatial_list, extrapolation_list, Year_Set = NULL, - Years2Include = NULL, zrange, ...) +plot_residuals( + Lat_i, + Lon_i, + TmbData, + Report, + Q, + projargs = "+proj=longlat", + working_dir = paste0(getwd(), "/"), + spatial_list, + extrapolation_list, + Year_Set = NULL, + Years2Include = NULL, + zrange, + ... +) } \arguments{ \item{Lat_i}{Latitude for every observation \code{i}} diff --git a/man/plot_results.Rd b/man/plot_results.Rd index 94bfe53..8eadab3 100644 --- a/man/plot_results.Rd +++ b/man/plot_results.Rd @@ -4,12 +4,24 @@ \alias{plot_results} \title{Plot results} \usage{ -plot_results(fit, settings = fit$settings, plot_set = 3, - working_dir = paste0(getwd(), "/"), year_labels = fit$year_labels, - years_to_plot = fit$years_to_plot, use_biascorr = TRUE, map_list, - category_names, check_residuals = TRUE, projargs = "+proj=longlat", - zrange, n_samples = 100, calculate_relative_to_average = FALSE, - type = 1, ...) +plot_results( + fit, + settings = fit$settings, + plot_set = 3, + working_dir = paste0(getwd(), "/"), + year_labels = fit$year_labels, + years_to_plot = fit$years_to_plot, + use_biascorr = TRUE, + map_list, + category_names, + check_residuals = TRUE, + projargs = "+proj=longlat", + zrange, + n_samples = 100, + calculate_relative_to_average = FALSE, + type = 1, + ... +) } \arguments{ \item{fit}{Output from \code{fit_model}} @@ -74,9 +86,10 @@ This function takes a fitted VAST model and generates a standard set of diagnost It does this by calling a series of mid-level plotting functions; see list of functions in Value section of documentation. } \seealso{ -\code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots +\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 -Other wrapper functions: \code{\link{fit_model}}, - \code{\link{make_settings}} +Other wrapper functions: +\code{\link{fit_model}()}, +\code{\link{make_settings}()} } \concept{wrapper functions} diff --git a/man/plot_variable.Rd b/man/plot_variable.Rd index 9d4c782..95eb3d5 100644 --- a/man/plot_variable.Rd +++ b/man/plot_variable.Rd @@ -4,14 +4,35 @@ \alias{plot_variable} \title{Plot results on a multi-panel figure of maps} \usage{ -plot_variable(Y_gt, map_list, panel_labels, projargs = "+proj=longlat", - map_resolution = "medium", file_name = "density", - working_dir = paste0(getwd(), "/"), Format = "png", Res = 200, - add = FALSE, outermargintext = c("Eastings", "Northings"), - zlim = NULL, col, mar = c(0, 0, 2, 0), oma = c(4, 4, 0, 0), - legend_x = c(0, 0.05), legend_y = c(0.05, 0.45), cex.legend = 1, - mfrow, land_color = "grey", n_cells, xlim, ylim, country = NULL, - contour_nlevels = 0, fun = mean, ...) +plot_variable( + Y_gt, + map_list, + panel_labels, + projargs = "+proj=longlat", + map_resolution = "medium", + file_name = "density", + working_dir = paste0(getwd(), "/"), + Format = "png", + Res = 200, + add = FALSE, + outermargintext = c("Eastings", "Northings"), + zlim = NULL, + col, + mar = c(0, 0, 2, 0), + oma = c(4, 4, 0, 0), + legend_x = c(0, 0.05), + legend_y = c(0.05, 0.45), + cex.legend = 1, + mfrow, + land_color = "grey", + n_cells, + xlim, + ylim, + country = NULL, + contour_nlevels = 0, + fun = mean, + ... +) } \arguments{ \item{Y_gt}{matrix where values for every column are plotted as a map} diff --git a/man/predict.fit_model.Rd b/man/predict.fit_model.Rd new file mode 100644 index 0000000..48432ea --- /dev/null +++ b/man/predict.fit_model.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fit_model.R +\name{predict.fit_model} +\alias{predict.fit_model} +\title{Predict density for new samples (\emph{Beta version; may change without notice})} +\usage{ +\method{predict}{fit_model}( + x, + what = "D_i", + Lat_i, + Lon_i, + t_i, + a_i, + c_iz = rep(0, length(t_i)), + v_i = rep(0, length(t_i)), + covariate_data = NULL, + catchability_data = NULL, + working_dir = paste0(getwd(), "/") +) +} +\arguments{ +\item{x}{Output from \code{\link{fit_model}}} + +\item{what}{Which output from \code{fit$Report} should be extracted; default is predicted density} + +\item{t_i}{Time for each observation i} + +\item{a_i}{Sampled area for each observation i; use \code{a_i=1} for observations without a natural area measurement, while +noting that resulting densities no longer have interpretable units in that case)} + +\item{c_iz}{Category (e.g., species, length-bin) for each observation i} + +\item{v_i}{sampling category (e.g., vessel or tow) associated with overdispersed variation for each observation i +(by default \code{v_i=0} for all samples, which will not affect things given the default values for \code{OverdispersionConfig})} + +\item{covariate_data}{data frame of covariate values with columns \code{Lat}, \code{Lon}, and \code{Year}, and other columns matching names in \code{formula}; \code{Year=NA} can be used for covariates that do not change among years (e.g., depth)} + +\item{catchability_data}{data-frame of covariates for use when specifying \code{Q1_formula} and \code{Q2_formula}} +} +\description{ +\code{predict.fit_model} calculates predictions given new data +} +\details{ +When the user does not supply \code{covariate_data}, predictions are based upon covariate values +interpolated from covariates supplied when fitting the model. When the user does supply \code{covariate_data}, +predictions are created when interpolating based on those new supplied values. +} diff --git a/man/project_coordinates.Rd b/man/project_coordinates.Rd index ec83c52..4d1e507 100644 --- a/man/project_coordinates.Rd +++ b/man/project_coordinates.Rd @@ -4,9 +4,15 @@ \alias{project_coordinates} \title{Convert from Lat-Long to UTM} \usage{ -project_coordinates(X, Y, projargs = NA, - origargs = "+proj=longlat +ellps=WGS84", zone = NA, - flip_around_dateline = FALSE, silent = FALSE) +project_coordinates( + X, + Y, + projargs = NA, + origargs = "+proj=longlat +datum=WGS84", + zone = NA, + flip_around_dateline = FALSE, + silent = FALSE +) } \arguments{ \item{X}{vector of horizontal-coordinates (e.g., longitude by default)} diff --git a/man/rotate_factors.Rd b/man/rotate_factors.Rd index 34b697c..f70443b 100644 --- a/man/rotate_factors.Rd +++ b/man/rotate_factors.Rd @@ -4,8 +4,14 @@ \alias{rotate_factors} \title{Rotate results} \usage{ -rotate_factors(Cov_jj = NULL, L_pj = NULL, Psi_sjt = NULL, - RotationMethod = "PCA", testcutoff = 1e-10, quiet = FALSE) +rotate_factors( + Cov_jj = NULL, + L_pj = NULL, + Psi_sjt = NULL, + RotationMethod = "PCA", + testcutoff = 1e-10, + quiet = FALSE +) } \arguments{ \item{Cov_jj}{Covariance calculated from loadings matrix} diff --git a/man/sample_variable.Rd b/man/sample_variable.Rd index 0ab8c4b..c00da15 100644 --- a/man/sample_variable.Rd +++ b/man/sample_variable.Rd @@ -4,8 +4,14 @@ \alias{sample_variable} \title{Sample from predictive distribution of a variable} \usage{ -sample_variable(Sdreport, Obj, variable_name, n_samples = 100, - sample_fixed = TRUE, seed = 123456) +sample_variable( + Sdreport, + Obj, + variable_name, + n_samples = 100, + sample_fixed = TRUE, + seed = 123456 +) } \arguments{ \item{Sdreport}{TMB output from `\code{TMB::sdreport(Obj)}`} diff --git a/man/smallPlot.Rd b/man/smallPlot.Rd index b540df5..be4dfc5 100644 --- a/man/smallPlot.Rd +++ b/man/smallPlot.Rd @@ -4,9 +4,22 @@ \alias{smallPlot} \title{Inset small plot within figure} \usage{ -smallPlot(expr, x = c(5, 70), y = c(50, 100), x1, y1, x2, y2, - mar = c(12, 14, 3, 3), mgp = c(1.8, 0.8, 0), bg = par("bg"), - border = par("fg"), las = 1, resetfocus = TRUE, ...) +smallPlot( + expr, + x = c(5, 70), + y = c(50, 100), + x1, + y1, + x2, + y2, + mar = c(12, 14, 3, 3), + mgp = c(1.8, 0.8, 0), + bg = par("bg"), + border = par("fg"), + las = 1, + resetfocus = TRUE, + ... +) } \value{ parameters of small plot, invisible. diff --git a/man/summarize_covariance.Rd b/man/summarize_covariance.Rd index 8f48130..6a63efd 100644 --- a/man/summarize_covariance.Rd +++ b/man/summarize_covariance.Rd @@ -4,11 +4,22 @@ \alias{summarize_covariance} \title{Explore spatio-temporal covariance} \usage{ -summarize_covariance(Report, Data, ParHat, SD = NULL, - category_order = 1:Data$n_c, category_names = 1:Data$n_c, - plotdir = paste0(getwd(), "/"), figname = "Cov", plotTF = NULL, - plot_cor = TRUE, mgp = c(2, 0.5, 0), tck = -0.02, oma = c(0, 5, - 2, 0), ...) +summarize_covariance( + Report, + Data, + ParHat, + SD = NULL, + category_order = 1:Data$n_c, + category_names = 1:Data$n_c, + plotdir = paste0(getwd(), "/"), + figname = "Cov", + plotTF = NULL, + plot_cor = TRUE, + mgp = c(2, 0.5, 0), + tck = -0.02, + oma = c(0, 5, 2, 0), + ... +) } \arguments{ \item{Report}{output report, e.g., from \code{Report <- obj$report()}} diff --git a/man/summary.fit_model.Rd b/man/summary.fit_model.Rd index c53fd79..7a4cbfe 100644 --- a/man/summary.fit_model.Rd +++ b/man/summary.fit_model.Rd @@ -4,8 +4,14 @@ \alias{summary.fit_model} \title{Extract summary of spatial estimates} \usage{ -\method{summary}{fit_model}(x, what = "density", n_samples = 250, - working_dir = NULL, type = 1, ...) +\method{summary}{fit_model}( + x, + what = "density", + n_samples = 250, + working_dir = NULL, + type = 1, + ... +) } \arguments{ \item{x}{Output from \code{\link{fit_model}}} From 4a89be53668145b9ddf27cf45c0d0389a1b08189 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Fri, 7 Aug 2020 18:01:54 -0700 Subject: [PATCH 02/15] small fix... ... to eliminate warning during devtools::document(.) --- NAMESPACE | 1 - R/deprecated.R | 1 + R/fit_model.R | 1 - R/plot_index.R | 2 +- 4 files changed, 2 insertions(+), 3 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index e83f2d8..e3b3b78 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -52,7 +52,6 @@ export(make_settings) export(make_spatial_info) export(map_hypervariance) export(match_strata_fn) -export(plot) export(plot_anisotropy) export(plot_biomass_index) export(plot_cov) diff --git a/R/deprecated.R b/R/deprecated.R index 96051f2..ba90f58 100644 --- a/R/deprecated.R +++ b/R/deprecated.R @@ -1022,6 +1022,7 @@ function( colvec, heatrange, textmargin=NULL, labeltransform="uniform", dopar=TR plot_lines = function( x, y, ybounds, fn=lines, col_bounds="black", bounds_type="whiskers", border=NA, border_lty="solid", lwd_bounds=1, ... ){ + # Function still used in plot_index #warning( "`plot_lines` is soft-deprecated" ) fn( y=y, x=x, ... ) diff --git a/R/fit_model.R b/R/fit_model.R index 7d54cd1..3a7d32c 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -260,7 +260,6 @@ print.fit_model <- function(x, ...) #' @param ... Arguments passed to \code{\link{plot_results}} #' @return NULL #' @method plot fit_model -#' @export plot #' @export plot.fit_model <- function(x, what="results", ...) { diff --git a/R/plot_index.R b/R/plot_index.R index c94dd7c..895ce99 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -6,7 +6,7 @@ #' @param Index_ctl A matrix or array of time-series estimates for multiple categories \code{c}, years \code{t}, and strata \code{l} #' @param sd_Index_ctl A matrix or array of variances for each estimate #' @inheritParams plot_biomass_index -#' @inheritParams plot_lines + #' @param Yrange lower and upper bound for left-hand y-axis, corresponding to input \code{Index_ctl} (use \code{Yrange[1]=NA} and/or \code{Yrange[2]=NA} for using the lower and upper bound of estimate intervals) #' @param Y2range lower and upper bound for right-hand y-axis, corresponding to input \code{SampleSize_ctz} (see Yrange for more info) #' @param SampleSize_ctz optional array of sample sizes for each category and year to be plotted on each panel From d3295eaf50eb90638ba1bd04183b4fcf19772785 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Fri, 7 Aug 2020 21:00:26 -0700 Subject: [PATCH 03/15] eliminate warning during document(.) --- DESCRIPTION | 1 + R/Calc_Polygon_Areas_and_Polygons_Fn.R | 1 + man/Coherence.Rd | 3 ++- man/calc_coherence.Rd | 3 ++- man/calc_synchrony.Rd | 3 ++- man/calculate_proportion.Rd | 8 +++++++- man/make_settings.Rd | 17 +++++++++++++---- 7 files changed, 28 insertions(+), 8 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 15f82b8..5d09ee1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -53,6 +53,7 @@ Remotes: License: GPL-3 LazyData: yes BuildVignettes: yes +Encoding: UTF-8 RoxygenNote: 7.1.1 URL: http://github.com/james-thorson/FishStatsUtils BugReports: http://github.com/james-thorson/FishStatsUtils/issues diff --git a/R/Calc_Polygon_Areas_and_Polygons_Fn.R b/R/Calc_Polygon_Areas_and_Polygons_Fn.R index 6e44531..87b997d 100644 --- a/R/Calc_Polygon_Areas_and_Polygons_Fn.R +++ b/R/Calc_Polygon_Areas_and_Polygons_Fn.R @@ -4,6 +4,7 @@ #' \code{Calc_Polygon_Areas_and_Polygons_Fn} builds outputs for a given triangulated mesh used for approximating spatial variation #' #' @inheritParams Calc_Anisotropic_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) #' @param a_el a matrix with \code{nrow(Data_Extrap)} rows and l columns, giving the area within the l-th stratum for the e-th row of Data_Extrap diff --git a/man/Coherence.Rd b/man/Coherence.Rd index 65892d3..7be0c8b 100644 --- a/man/Coherence.Rd +++ b/man/Coherence.Rd @@ -18,7 +18,8 @@ Coherence( \item{covhat}{estimated covariance used for calculating coherence} -\item{yearbounds_zz}{matrix with two columns, giving first and last years for defining one or more periods (rows) used to calculate changes in synchrony over time (only used if \code{Options['Calculate_Synchrony']=1})} +\item{yearbounds_zz}{matrix with two columns, giving first and last years for defining one or more periods (rows) used to +calculate changes in synchrony over time (only used if \code{Options['Calculate_Synchrony']=1})} } \value{ Tagged list containing measures of synchrony diff --git a/man/calc_coherence.Rd b/man/calc_coherence.Rd index 8849c0a..86be873 100644 --- a/man/calc_coherence.Rd +++ b/man/calc_coherence.Rd @@ -18,7 +18,8 @@ calc_coherence( \item{covhat}{estimated covariance used for calculating coherence} -\item{yearbounds_zz}{matrix with two columns, giving first and last years for defining one or more periods (rows) used to calculate changes in synchrony over time (only used if \code{Options['Calculate_Synchrony']=1})} +\item{yearbounds_zz}{matrix with two columns, giving first and last years for defining one or more periods (rows) used to +calculate changes in synchrony over time (only used if \code{Options['Calculate_Synchrony']=1})} } \value{ Tagged list containing measures of synchrony diff --git a/man/calc_synchrony.Rd b/man/calc_synchrony.Rd index a987461..43e1f37 100644 --- a/man/calc_synchrony.Rd +++ b/man/calc_synchrony.Rd @@ -11,7 +11,8 @@ calc_synchrony(Report, Data, yearbounds_zz = matrix(c(1, Data$n_t), nrow = 1)) \item{Data}{tagged list of input data} -\item{yearbounds_zz}{matrix with two columns, giving first and last years for defining one or more periods (rows) used to calculate changes in synchrony over time (only used if \code{Options['Calculate_Synchrony']=1})} +\item{yearbounds_zz}{matrix with two columns, giving first and last years for defining one or more periods (rows) used to +calculate changes in synchrony over time (only used if \code{Options['Calculate_Synchrony']=1})} } \value{ Tagged list containing measures of synchrony diff --git a/man/calculate_proportion.Rd b/man/calculate_proportion.Rd index e102651..592e2ae 100644 --- a/man/calculate_proportion.Rd +++ b/man/calculate_proportion.Rd @@ -29,7 +29,13 @@ calculate_proportion( \item{Index}{output from \code{FishStatsUtils::plot_biomass_index}} -\item{Expansion_cz}{matrix specifying how densities are expanded when calculating annual indices, with a row for each category \code{c} and two columns. The first column specifies whether to calculate annual index for category \code{c} as the weighted-sum across density estimates, where density is weighted by area ("area-weighted expansion", \code{Expansion[c,1]=0}, the default), where density is weighted by the expanded value for another category ("abundance weighted expansion" \code{Expansion[c1,1]=1}), or the index is calculated as the weighted average of density weighted by the expanded value for another category ("abundance weighted-average expansion" \code{Expansion[c1,1]=2}). The 2nd column is only used when \code{Expansion[c1,1]=1} or \code{Expansion[c1,1]=2} , and specifies the category to use for abundance-weighted expansion, where \code{Expansion[c1,2]=c2} and \code{c2} must be lower than \code{c1}.} +\item{Expansion_cz}{matrix specifying how densities are expanded when calculating annual indices, with a row for each category \code{c} and two columns. +The first column specifies whether to calculate annual index for category \code{c} as the weighted-sum across density estimates, +where density is weighted by area ("area-weighted expansion", \code{Expansion[c,1]=0}, the default), +where density is weighted by the expanded value for another category ("abundance weighted expansion" \code{Expansion[c1,1]=1}), +or the index is calculated as the weighted average of density weighted by the expanded value for another category +("abundance weighted-average expansion" \code{Expansion[c1,1]=2}). The 2nd column is only used when \code{Expansion[c1,1]=1} or \code{Expansion[c1,1]=2}, +and specifies the category to use for abundance-weighted expansion, where \code{Expansion[c1,2]=c2} and \code{c2} must be lower than \code{c1}.} \item{Year_Set}{Year names for labeling panels} diff --git a/man/make_settings.Rd b/man/make_settings.Rd index 304391a..01689c3 100644 --- a/man/make_settings.Rd +++ b/man/make_settings.Rd @@ -45,17 +45,26 @@ Epsilon refers to spatio-temporal variation, Omega1 refers to variation in encou Omega2 refers to variation in positive catch rates, where 0 is off, "AR1" is an AR1 process, and >0 is the number of elements in a factor-analysis covariance} -\item{RhoConfig}{vector of form c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsilon2"=0) specifying whether either intercepts (Beta1 and Beta2) or spatio-temporal variation (Epsilon1 and Epsilon2) is structured among time intervals (0: each year as fixed effect; 1: each year as random following IID distribution; 2: each year as random following a random walk; 3: constant among years as fixed effect; 4: each year as random following AR1 process); If missing, assumed to be zero for each element} +\item{RhoConfig}{vector of form c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsilon2"=0) specifying whether either intercepts (Beta1 and Beta2) +or spatio-temporal variation (Epsilon1 and Epsilon2) is structured among time intervals (0: each year as fixed effect; +1: each year as random following IID distribution; 2: each year as random following a random walk; +3: constant among years as fixed effect; 4: each year as random following AR1 process); If missing, assumed to be zero for each element} -\item{OverdispersionConfig}{a vector of format c("eta1"=0, "eta2"="AR1") governing any correlated overdispersion among categories for each level of v_i, where eta1 is for encounter probability, and eta2 is for positive catch rates, where 0 is off, "AR1" is an AR1 process, and >0 is the number of elements in a factor-analysis covariance (by default, c("eta1"=0, "eta2"=0) and this turns off overdispersion)} +\item{OverdispersionConfig}{a vector of format c("eta1"=0, "eta2"="AR1") governing any correlated overdispersion +among categories for each level of v_i, where eta1 is for encounter probability, and eta2 is for positive catch rates, +where 0 is off, "AR1" is an AR1 process, and >0 is the number of elements in a factor-analysis covariance (by default, +c("eta1"=0, "eta2"=0) and this turns off overdispersion)} -\item{Options}{a vector of form c('SD_site_logdensity'=FALSE,'Calculate_Range'=FALSE,'Calculate_effective_area'=FALSE,'Calculate_Cov_SE'=FALSE,'Calculate_Synchrony'=FALSE,'Calculate_proportion'=FALSE), where Calculate_Range=1 turns on calculation of center of gravity, and Calculate_effective_area=1 turns on calculation of effective area occupied} +\item{Options}{a vector of form c('SD_site_logdensity'=FALSE,'Calculate_Range'=FALSE,'Calculate_effective_area'=FALSE, +'Calculate_Cov_SE'=FALSE,'Calculate_Synchrony'=FALSE,'Calculate_proportion'=FALSE), where \code{Calculate_Range=1} +turns on calculation of center of gravity, and Calculate_effective_area=1 turns on calculation of effective area occupied} \item{use_anisotropy}{Boolean indicating whether to estimate two additional parameters representing geometric anisotropy} \item{vars_to_correct}{a character-vector listing which parameters to include for bias-correction, as passed to \code{TMBhelper::Optimize}} -\item{Version}{a version number; If missing, defaults to latest version using \code{FishStatsUtils::get_latest_version(package="VAST")}} +\item{Version}{Which CPP version to use. If missing, defaults to latest version using \code{FishStatsUtils::get_latest_version(package="VAST")}. +Can be used to specify using an older CPP, to maintain backwards compatibility.} \item{treat_nonencounter_as_zero}{Boolean indicating whether to treat any year-category combination as having zero biomass when generating abundance indices and resulting compositional estimates} From ae27c8a463421ff5a43bd7d95647ac7f5388bf1d Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Sun, 9 Aug 2020 15:08:53 -0700 Subject: [PATCH 04/15] fix make_settings(.)... ... to track change from Index_cyl to Index_ctl --- R/make_settings.R | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/R/make_settings.R b/R/make_settings.R index 7d39815..81851cf 100644 --- a/R/make_settings.R +++ b/R/make_settings.R @@ -52,7 +52,7 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE, if(missing(bias.correct)) bias.correct = TRUE if(missing(treat_nonencounter_as_zero)) treat_nonencounter_as_zero = TRUE if(missing(Options)) Options = c("SD_site_logdensity"=FALSE, "Calculate_Range"=TRUE, "Calculate_effective_area"=TRUE, "treat_nonencounter_as_zero"=treat_nonencounter_as_zero ) - if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" ) + if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Index_ctl" ) if(missing(knot_method)) knot_method = "samples" if(missing(max_cells)) max_cells = Inf } @@ -72,7 +72,7 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE, if(missing(bias.correct)) bias.correct = TRUE if(missing(treat_nonencounter_as_zero)) treat_nonencounter_as_zero = TRUE if(missing(Options)) Options = c("SD_site_logdensity"=FALSE, "Calculate_Range"=TRUE, "Calculate_effective_area"=TRUE, "treat_nonencounter_as_zero"=treat_nonencounter_as_zero ) - if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" ) + if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Index_ctl" ) if(missing(knot_method)) knot_method = "grid" if(missing(max_cells)) max_cells = max( 2000, n_x*10 ) } @@ -95,7 +95,7 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE, if(missing(treat_nonencounter_as_zero)) treat_nonencounter_as_zero = FALSE if(missing(Options)) Options = c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=TRUE, "treat_nonencounter_as_zero"=treat_nonencounter_as_zero, "Project_factors"=TRUE ) - if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" ) + if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Index_ctl" ) if(missing(knot_method)) knot_method = "samples" if(missing(max_cells)) max_cells = Inf } @@ -118,7 +118,7 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE, if(missing(bias.correct)) bias.correct = TRUE if(missing(Options)) Options = c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=FALSE, "Calculate_Fratio"=TRUE, "Estimate_B0"=TRUE ) - if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Bratio_cyl" ) + if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Bratio_cyl", "Index_ctl", "Bratio_ctl" ) if(missing(knot_method)) knot_method = "samples" if(missing(max_cells)) max_cells = Inf } @@ -140,7 +140,7 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE, if(missing(ObsModel)) ObsModel = c(1,1) if(missing(bias.correct)) bias.correct = FALSE if(missing(Options)) Options = c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=TRUE, "Project_factors"=TRUE ) - if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" ) + if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Index_ctl" ) if(missing(knot_method)) knot_method = "samples" if(missing(max_cells)) max_cells = Inf } @@ -166,7 +166,7 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE, if(missing(ObsModel)) ObsModel = c(1,1) if(missing(bias.correct)) bias.correct = FALSE if(missing(Options)) Options = c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=TRUE, "Project_factors"=TRUE ) - if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" ) + if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Index_ctl" ) if(missing(knot_method)) knot_method = "samples" if(missing(max_cells)) max_cells = Inf } @@ -183,7 +183,7 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE, if(missing(ObsModel)) ObsModel = c(2,1) if(missing(bias.correct)) bias.correct = FALSE if(missing(Options)) Options = c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=TRUE, "Project_factors"=TRUE ) - if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" ) + if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Index_ctl" ) if(missing(knot_method)) knot_method = "grid" if(missing(max_cells)) max_cells = max( 2000, n_x*10 ) } From 3fd2dd333592e1ca6f49b1ef811b7083b1325b17 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Mon, 17 Aug 2020 14:30:34 -0700 Subject: [PATCH 05/15] fix small bug... ... wherein DHARMa residuals calculation was changing n_g in the original object outside of function due to way TMB passes objects --- R/fit_model.R | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/R/fit_model.R b/R/fit_model.R index 3a7d32c..81e33b6 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -356,6 +356,11 @@ summary.fit_model <- function(x, what="density", n_samples=250, if( tolower(what) == "residuals" ){ # extract objects Obj = x$tmb_list$Obj + + # Change n_g + # Must change back explicitly because TMB appears to pass env as a pointer, so changes in copy affect original x outside of function! + n_g_orig = Obj$env$data$n_g + on.exit( Obj$env$data$n_g = n_g_orig ) Obj$env$data$n_g = 0 # check for issues From e025a81ad1c5ea1618cf0b61f3511f08fd86ac8d Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Mon, 17 Aug 2020 15:04:24 -0700 Subject: [PATCH 06/15] fix bug in last commit --- R/fit_model.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/fit_model.R b/R/fit_model.R index 81e33b6..55e3b36 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -360,7 +360,8 @@ summary.fit_model <- function(x, what="density", n_samples=250, # Change n_g # Must change back explicitly because TMB appears to pass env as a pointer, so changes in copy affect original x outside of function! n_g_orig = Obj$env$data$n_g - on.exit( Obj$env$data$n_g = n_g_orig ) + Tmp = function(n_g){Obj$env$data$n_g = n_g} + on.exit( Tmp(n_g_orig) ) Obj$env$data$n_g = 0 # check for issues From cc718be2c0f7483091910d7206d3f44f67847ad0 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Wed, 19 Aug 2020 08:12:37 -0700 Subject: [PATCH 07/15] revert settings during simulate_data(.) ... ... to deal with way TMB uses pointers within TMB object, such that changes in local copy within function change original object outside of function --- R/fit_model.R | 4 ++-- R/simulate_data.R | 6 ++++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/R/fit_model.R b/R/fit_model.R index 55e3b36..7618641 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -360,8 +360,8 @@ summary.fit_model <- function(x, what="density", n_samples=250, # Change n_g # Must change back explicitly because TMB appears to pass env as a pointer, so changes in copy affect original x outside of function! n_g_orig = Obj$env$data$n_g - Tmp = function(n_g){Obj$env$data$n_g = n_g} - on.exit( Tmp(n_g_orig) ) + revert_settings = function(n_g){Obj$env$data$n_g = n_g} + on.exit( revert_settings(n_g_orig) ) Obj$env$data$n_g = 0 # check for issues diff --git a/R/simulate_data.R b/R/simulate_data.R index 3119304..7b9e35d 100644 --- a/R/simulate_data.R +++ b/R/simulate_data.R @@ -53,9 +53,15 @@ simulate_data = function( fit, type=1, random_seed=NULL ){ # Extract stuff Obj = fit$tmb_list$Obj + simulate_random_effects_orig = Obj$env$data$Options_list$Options['simulate_random_effects'] + + # Revert settings when done + revert_settings = function(simulate_random_effects){Obj$env$data$Options_list$Options['simulate_random_effects'] = simulate_random_effects} + on.exit( revert_settings(simulate_random_effects_orig) ) # Simulate conditional upon fixed and random effect estimates if( type==1 ){ + # Change and revert settings Obj$env$data$Options_list$Options['simulate_random_effects'] = FALSE set.seed(random_seed) Return = Obj$simulate( complete=TRUE ) From 31bdf5adb1ff646c8a1350294ac0624690c824b8 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Fri, 28 Aug 2020 21:04:29 -0700 Subject: [PATCH 08/15] correct .prj projection for ICES shapefiles --- inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.prj | 2 +- inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.prj | 2 +- inst/region_shapefiles/BITS/Shapefile.prj | 2 +- inst/region_shapefiles/BTS-VIIa/Shapefile.prj | 2 +- inst/region_shapefiles/BTS/Shapefile.prj | 2 +- inst/region_shapefiles/EVHOE/Shapefile.prj | 2 +- inst/region_shapefiles/IE-IGFS/Shapefile.prj | 2 +- inst/region_shapefiles/NIGFS/Shapefile.prj | 2 +- inst/region_shapefiles/NS_IBTS/Shapefile.prj | 2 +- inst/region_shapefiles/PT-IBTS/Shapefile.prj | 2 +- inst/region_shapefiles/SP-ARSA/Shapefile.prj | 2 +- inst/region_shapefiles/SP-NORTH/Shapefile.prj | 2 +- inst/region_shapefiles/SP-PORC/Shapefile.prj | 2 +- 13 files changed, 13 insertions(+), 13 deletions(-) diff --git a/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.prj b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.prj +++ b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.prj b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.prj +++ b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/BITS/Shapefile.prj b/inst/region_shapefiles/BITS/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/BITS/Shapefile.prj +++ b/inst/region_shapefiles/BITS/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/BTS-VIIa/Shapefile.prj b/inst/region_shapefiles/BTS-VIIa/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/BTS-VIIa/Shapefile.prj +++ b/inst/region_shapefiles/BTS-VIIa/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/BTS/Shapefile.prj b/inst/region_shapefiles/BTS/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/BTS/Shapefile.prj +++ b/inst/region_shapefiles/BTS/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/EVHOE/Shapefile.prj b/inst/region_shapefiles/EVHOE/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/EVHOE/Shapefile.prj +++ b/inst/region_shapefiles/EVHOE/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/IE-IGFS/Shapefile.prj b/inst/region_shapefiles/IE-IGFS/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/IE-IGFS/Shapefile.prj +++ b/inst/region_shapefiles/IE-IGFS/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/NIGFS/Shapefile.prj b/inst/region_shapefiles/NIGFS/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/NIGFS/Shapefile.prj +++ b/inst/region_shapefiles/NIGFS/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/NS_IBTS/Shapefile.prj b/inst/region_shapefiles/NS_IBTS/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/NS_IBTS/Shapefile.prj +++ b/inst/region_shapefiles/NS_IBTS/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/PT-IBTS/Shapefile.prj b/inst/region_shapefiles/PT-IBTS/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/PT-IBTS/Shapefile.prj +++ b/inst/region_shapefiles/PT-IBTS/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/SP-ARSA/Shapefile.prj b/inst/region_shapefiles/SP-ARSA/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/SP-ARSA/Shapefile.prj +++ b/inst/region_shapefiles/SP-ARSA/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/SP-NORTH/Shapefile.prj b/inst/region_shapefiles/SP-NORTH/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/SP-NORTH/Shapefile.prj +++ b/inst/region_shapefiles/SP-NORTH/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/inst/region_shapefiles/SP-PORC/Shapefile.prj b/inst/region_shapefiles/SP-PORC/Shapefile.prj index 6e3eaed..a41f629 100644 --- a/inst/region_shapefiles/SP-PORC/Shapefile.prj +++ b/inst/region_shapefiles/SP-PORC/Shapefile.prj @@ -1 +1 @@ -PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] \ No newline at end of file +GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file From 0d64326ae888c2f3be32890b2c6ee9c252441c30 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Mon, 14 Sep 2020 12:05:11 -0700 Subject: [PATCH 09/15] add message to fit_model explaining location of files --- R/convert_shapefile.R | 1 + R/fit_model.R | 1 + 2 files changed, 2 insertions(+) diff --git a/R/convert_shapefile.R b/R/convert_shapefile.R index 6c29477..3214753 100644 --- a/R/convert_shapefile.R +++ b/R/convert_shapefile.R @@ -24,6 +24,7 @@ convert_shapefile = function( file_path, projargs=NULL, grid_dim_km=c(2,2), proj make_plots=FALSE, quiet=TRUE, area_tolerance=0.05, ... ){ shapefile_orig = rgdal::readOGR( file_path, verbose=FALSE, p4s=projargs_for_shapefile ) + # raster::shapefile(.) has simplified read-write interface for future reference if( !(shapefile_orig@class %in% c("SpatialPolygonsDataFrame","SpatialPolygons")) ){ warning( "object at `file_path` doesn't appear to be a shapefile") } diff --git a/R/fit_model.R b/R/fit_model.R index 7618641..d7f686e 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -110,6 +110,7 @@ fit_model = function( settings, Lat_i, Lon_i, t_i, b_i, a_i, c_iz=rep(0,length(b years_to_plot = which( year_labels %in% t_i ) # Save record + message("\n### Writing output from `fit_model` in directory: ", working_dir) dir.create(working_dir, showWarnings=FALSE, recursive=TRUE) #save( settings, file=file.path(working_dir,"Record.RData")) capture.output( settings, file=file.path(working_dir,"settings.txt")) From 963d300cb87873c3c14670fd0c09d8c2524bbb00 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Mon, 14 Sep 2020 12:19:12 -0700 Subject: [PATCH 10/15] provide example-usage for sample_variable --- R/fit_model.R | 1 - R/sample_variable.R | 14 ++++++++++++-- man/fit_model.Rd | 1 - man/sample_variable.Rd | 12 ++++++++++++ 4 files changed, 24 insertions(+), 4 deletions(-) diff --git a/R/fit_model.R b/R/fit_model.R index d7f686e..f6ec923 100644 --- a/R/fit_model.R +++ b/R/fit_model.R @@ -61,7 +61,6 @@ #' @examples #' \dontrun{ #' # Load packages -#' library(TMB) #' library(VAST) #' #' # load data set diff --git a/R/sample_variable.R b/R/sample_variable.R index 7919cff..64f0a2b 100644 --- a/R/sample_variable.R +++ b/R/sample_variable.R @@ -14,7 +14,17 @@ #' @param seed integer used to set random-number seed when sampling variables, as passed to \code{set.seed(.)} #' @param sample_fixed whether to sample fixed and random effects, \code{sample_fixed=TRUE} as by default, or just sample random effects, \code{sample_fixed=FALSE} #' - +#' @examples +#' \dontrun{ +#' # Run model using selected inputs, but also with getJointPrecision=TRUE +#' fit = fit_model( ..., +#' getJointPrecision=TRUE ) +#' +#' # Run sample_variable +#' sample = sample_variable( Sdreport=fit$parameter_estimates$SD, +#' Obj=fit$tmb_list$Obj, variable_name="D_gct" ) +#' } +#' #' @export sample_variable = function( Sdreport, Obj, variable_name, n_samples=100, sample_fixed=TRUE, seed=123456 ){ @@ -24,7 +34,7 @@ sample_variable = function( Sdreport, Obj, variable_name, n_samples=100, stop("jointPrecision not present in Sdreport; please re-run with `getJointPrecision=TRUE`") } if( !(variable_name %in% names(Obj$report())) ){ - stop( variable_name, " not found in report(.); please choose check your requested variable name" ) + stop( variable_name, " not found in report(.); please choose check your requested variable name from available list: ", paste(names(Obj$report()),collapse=", ") ) } #### Local function diff --git a/man/fit_model.Rd b/man/fit_model.Rd index 9398afe..cb32e7a 100644 --- a/man/fit_model.Rd +++ b/man/fit_model.Rd @@ -139,7 +139,6 @@ Category \code{g} corresponds to values in \code{fit$data_list$g_i}. \examples{ \dontrun{ # Load packages -library(TMB) library(VAST) # load data set diff --git a/man/sample_variable.Rd b/man/sample_variable.Rd index c00da15..5bff83d 100644 --- a/man/sample_variable.Rd +++ b/man/sample_variable.Rd @@ -31,3 +31,15 @@ sample_variable( Exploratory testing suggests that calling \code{Obj$report(newvalues)} within a function like \code{sample_variable(.)} doesn't affect the results of subsequent calls of \code{Obj$report()} so this function does not appear to change anything in the global environment } +\examples{ +\dontrun{ +# Run model using selected inputs, but also with getJointPrecision=TRUE +fit = fit_model( ..., + getJointPrecision=TRUE ) + +# Run sample_variable +sample = sample_variable( Sdreport=fit$parameter_estimates$SD, + Obj=fit$tmb_list$Obj, variable_name="D_gct" ) +} + +} From 405b20cce22303d7176b011cf8d1d5770b4c73ac Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Sun, 20 Sep 2020 14:41:19 -0700 Subject: [PATCH 11/15] improved docs --- man/fit_model.Rd | 13 ++++++++----- man/make_settings.Rd | 5 +---- man/predict.fit_model.Rd | 11 +++++++---- 3 files changed, 16 insertions(+), 13 deletions(-) diff --git a/man/fit_model.Rd b/man/fit_model.Rd index cb32e7a..4c0b569 100644 --- a/man/fit_model.Rd +++ b/man/fit_model.Rd @@ -39,16 +39,19 @@ fit_model( \item{Lon_i}{Longitude for each sample} -\item{t_i}{Time for each observation i} +\item{t_i}{Vector of integers, providing the time (e.g., calendar year) for each observation i} -\item{b_i}{Sampled value (biomass, counts, etc.) for each observation i} +\item{b_i}{Numeric vector, providing sampled value (biomass, counts, etc.) for each observation i} -\item{a_i}{Sampled area for each observation i; use \code{a_i=1} for observations without a natural area measurement, while +\item{a_i}{Numeric vector containing values greater than zero, providing sampled area for each +observation i; use \code{a_i=1} for observations without a natural area measurement, while noting that resulting densities no longer have interpretable units in that case)} -\item{c_iz}{Category (e.g., species, length-bin) for each observation i} +\item{c_iz}{Vector of integers ranging from 0 to the number of variables minus 1, providing the +category (e.g., species, length-bin) for each observation i} -\item{v_i}{sampling category (e.g., vessel or tow) associated with overdispersed variation for each observation i +\item{v_i}{Vector of integers ranging from 0 to the number of vessels minus 1, +providing sampling category (e.g., vessel or tow) associated with overdispersed variation for each observation i (by default \code{v_i=0} for all samples, which will not affect things given the default values for \code{OverdispersionConfig})} \item{X1config_cp}{matrix of settings for each density covariate for the 1st lienar predictor, diff --git a/man/make_settings.Rd b/man/make_settings.Rd index 01689c3..e9b583f 100644 --- a/man/make_settings.Rd +++ b/man/make_settings.Rd @@ -40,10 +40,7 @@ make_settings( \item{zone}{UTM zone used for projecting Lat-Lon to km distances; use \code{zone=NA} by default to automatically detect UTM zone from the location of extrapolation-grid samples} -\item{FieldConfig}{a vector of format c("Omega1"=0, "Epsilon1"=10, "Omega2"="AR1", "Epsilon2"=10), where Omega refers to spatial variation, -Epsilon refers to spatio-temporal variation, Omega1 refers to variation in encounter probability, and -Omega2 refers to variation in positive catch rates, where 0 is off, "AR1" is an AR1 process, -and >0 is the number of elements in a factor-analysis covariance} +\item{FieldConfig}{See Details section of \code{\link[VAST]{make_data}} for details} \item{RhoConfig}{vector of form c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsilon2"=0) specifying whether either intercepts (Beta1 and Beta2) or spatio-temporal variation (Epsilon1 and Epsilon2) is structured among time intervals (0: each year as fixed effect; diff --git a/man/predict.fit_model.Rd b/man/predict.fit_model.Rd index 48432ea..95ad39d 100644 --- a/man/predict.fit_model.Rd +++ b/man/predict.fit_model.Rd @@ -23,14 +23,17 @@ \item{what}{Which output from \code{fit$Report} should be extracted; default is predicted density} -\item{t_i}{Time for each observation i} +\item{t_i}{Vector of integers, providing the time (e.g., calendar year) for each observation i} -\item{a_i}{Sampled area for each observation i; use \code{a_i=1} for observations without a natural area measurement, while +\item{a_i}{Numeric vector containing values greater than zero, providing sampled area for each +observation i; use \code{a_i=1} for observations without a natural area measurement, while noting that resulting densities no longer have interpretable units in that case)} -\item{c_iz}{Category (e.g., species, length-bin) for each observation i} +\item{c_iz}{Vector of integers ranging from 0 to the number of variables minus 1, providing the +category (e.g., species, length-bin) for each observation i} -\item{v_i}{sampling category (e.g., vessel or tow) associated with overdispersed variation for each observation i +\item{v_i}{Vector of integers ranging from 0 to the number of vessels minus 1, +providing sampling category (e.g., vessel or tow) associated with overdispersed variation for each observation i (by default \code{v_i=0} for all samples, which will not affect things given the default values for \code{OverdispersionConfig})} \item{covariate_data}{data frame of covariate values with columns \code{Lat}, \code{Lon}, and \code{Year}, and other columns matching names in \code{formula}; \code{Year=NA} can be used for covariates that do not change among years (e.g., depth)} From 899ed37e4ace674e4ed07f2caa4a8fc21e93b04a Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Mon, 21 Sep 2020 13:02:12 -0700 Subject: [PATCH 12/15] fix Region="User" to allow using depth-based post-strata --- R/Prepare_XXXX_Extrapolation_Data_Fn.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/Prepare_XXXX_Extrapolation_Data_Fn.R b/R/Prepare_XXXX_Extrapolation_Data_Fn.R index bd62074..c7215ae 100644 --- a/R/Prepare_XXXX_Extrapolation_Data_Fn.R +++ b/R/Prepare_XXXX_Extrapolation_Data_Fn.R @@ -642,7 +642,7 @@ function( input_grid, strata.limits=NULL, projargs=NA, zone=NA, flip_around_date # Augment with strata for each extrapolation cell Tmp = cbind("BEST_LAT_DD"=Data_Extrap[,'Lat'], "BEST_LON_DD"=Data_Extrap[,'Lon']) if( "Depth" %in% colnames(Data_Extrap) ){ - Tmp = cbind( Tmp, Data_Extrap[,'Depth'] ) + Tmp = cbind( Tmp, "BEST_DEPTH_M" = Data_Extrap[,'Depth'] ) } a_el = as.data.frame(matrix(NA, nrow=nrow(Data_Extrap), ncol=nrow(strata.limits), dimnames=list(NULL,strata.limits[,'STRATA']))) for(l in 1:ncol(a_el)){ From 15f597cf16bc4ae4ebfa108b54e5d53fa4b52821 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Mon, 21 Sep 2020 21:38:19 -0700 Subject: [PATCH 13/15] improve stability/readability for make_covariates... ... by first processing formula and then assembling X_ip and X_gtp (rather than the opposite order previously) --- R/make_covariates.R | 91 +++++++++++++++++++++------------------------ 1 file changed, 42 insertions(+), 49 deletions(-) diff --git a/R/make_covariates.R b/R/make_covariates.R index 2e03b9d..b2df52e 100644 --- a/R/make_covariates.R +++ b/R/make_covariates.R @@ -28,83 +28,76 @@ #' @export make_covariates = function( formula, covariate_data, Year_i, spatial_list, extrapolation_list ){ - # Errors + # Check for bad entries if( !is.data.frame(covariate_data) ) stop("Please ensure that `covariate_data` is a data frame") if( !all(c("Lat","Lon","Year") %in% names(covariate_data)) ){ stop( "`data` in `make_covariates(.)` must include columns `Lat`, `Lon`, and `Year`" ) } - # transform data inputs - sample_data = data.frame( "Year"=Year_i, "Lat"=spatial_list$latlon_i[,'Lat'], "Lon"=spatial_list$latlon_i[,'Lon'] ) - covariate_names = setdiff( names(covariate_data), names(sample_data) ) - # set of years needed Year_Set = min(Year_i):max(Year_i) + # Make `covariate_df` by expanding for rows with Year=NA + covariate_df = covariate_data[ which(!is.na(covariate_data[,'Year'])), ] + for( tI in seq_along(Year_Set) ){ + newrows = covariate_data[ which(is.na(covariate_data[,'Year'])), ] + newrows[,"Year"] = Year_Set[tI] + covariate_df = rbind( covariate_df, newrows ) + } + + # Make model.matrix + # To ensure identifiability given betas (intercepts), add intercept to formula + # and then remove that term from model.matrix. This will not fix identifiability + # issues arising when both conditions are met: + # factor(Year) has an interaction with another factor, and + # betas vary among years (are not constant) + Model_matrix = model.matrix( update.formula(formula, ~.+1), data=covariate_df ) + Columns_to_keep = which( attr(Model_matrix,"assign") != 0 ) + coefficient_names = attr(Model_matrix,"dimnames")[[2]][Columns_to_keep] + X = Model_matrix[,Columns_to_keep,drop=FALSE] + dimnames(X) = list(NULL, coefficient_names) + + # transform data inputs + sample_i = data.frame( "Year"=Year_i, "Lat"=spatial_list$latlon_i[,'Lat'], "Lon"=spatial_list$latlon_i[,'Lon'] ) + #covariate_names = setdiff( names(covariate_data), names(sample_data) ) + # extract latitude and longitude for extrapolation grid latlon_g = spatial_list$latlon_g # Create data frame of necessary size - DF_zp = NULL - #DF_ip = cbind( sample_data, covariate_data[rep(1,nrow(sample_data)),covariate_names] ) - DF_ip = data.frame( sample_data, covariate_data[rep(1,nrow(sample_data)),covariate_names] ) - colnames(DF_ip) = c( names(sample_data) ,covariate_names ) - #DF_ip[,covariate_names] = NA + X_gtp = array( NA, dim=c(nrow(latlon_g),length(Year_Set),ncol(X)), dimnames=list(NULL,Year_Set,colnames(X)) ) + X_ip = array( NA, dim=c(nrow(sample_i),ncol(X)), dimnames=list(NULL,colnames(X)) ) # Loop through data and extrapolation-grid for( tI in seq_along(Year_Set) ){ # Subset to same year - tmp_covariate_data = covariate_data[ which(Year_Set[tI]==covariate_data[,'Year'] | is.na(covariate_data[,'Year'])), , drop=FALSE] - if( nrow(tmp_covariate_data)==0 ){ + tmp_covariate_df = covariate_df[ which(Year_Set[tI]==covariate_df[,'Year']), , drop=FALSE] + tmp_X = X[ which(Year_Set[tI]==covariate_df[,'Year']), , drop=FALSE] + if( nrow(tmp_covariate_df)==0 ){ stop("Year ", Year_Set[tI], " not found in `covariate_data` please specify covariate values for all years" ) } - # - Which = which(Year_Set[tI]==sample_data[,'Year']) + + # Fill in values in X_ip + Which = which(Year_Set[tI]==sample_i[,'Year']) # Do nearest neighbors to define covariates for observations, skipping years without observations if( length(Which) > 0 ){ - NN = RANN::nn2( data=tmp_covariate_data[,c("Lat","Lon")], query=sample_data[Which,c("Lat","Lon")], k=1 ) - # Add to data-frame - nearest_covariates = tmp_covariate_data[ NN$nn.idx[,1], covariate_names, drop=FALSE ] - DF_ip[Which, covariate_names] = nearest_covariates + NN = RANN::nn2( data=tmp_covariate_df[,c("Lat","Lon")], query=sample_i[Which,c("Lat","Lon")], k=1 ) + # Fill in values + X_ip[Which, ] = tmp_X[ NN$nn.idx[,1], , drop=FALSE] } # Do nearest neighbors to define covariates for extrapolation grid, including years without observations - NN = RANN::nn2( data=tmp_covariate_data[,c("Lat","Lon")], query=latlon_g[,c("Lat","Lon")], k=1 ) + NN = RANN::nn2( data=tmp_covariate_df[,c("Lat","Lon")], query=latlon_g[,c("Lat","Lon")], k=1 ) # Add rows - nearest_covariates = tmp_covariate_data[ NN$nn.idx[,1], covariate_names, drop=FALSE ] - newrows = cbind("Year"=Year_Set[tI], latlon_g, nearest_covariates ) - DF_zp = rbind( DF_zp, newrows ) + X_gtp[,tI,] = tmp_X[ NN$nn.idx[,1], , drop=FALSE ] } - # Convert to dimensions requested - DF = rbind( DF_ip, DF_zp ) - - # Make model.matrix - # To ensure identifiability given betas (intercepts), add intercept to formula - # and then remove that term from model.matrix. This will not fix identifiability - # issues arising when both conditions are met: - # factor(Year) has an interaction with another factor, and - # betas vary among years (are not constant) - Model_matrix = model.matrix( update.formula(formula, ~.+1), data=DF ) - Columns_to_keep = which( attr(Model_matrix,"assign") != 0 ) - coefficient_names = attr(Model_matrix,"dimnames")[[2]][Columns_to_keep] - X = Model_matrix[,Columns_to_keep,drop=FALSE] - - # Make X_ip - X_ip = X[ 1:nrow(DF_ip), , drop=FALSE ] + # Make X_itp X_itp = aperm( X_ip %o% rep(1,length(Year_Set)), perm=c(1,3,2) ) - if( any(is.na(X_itp)) ) stop("Problem with `X_itp` in `make_covariates(.)") - # Make X_gpt and then permute dimensions - X_gpt = NULL - indices = nrow(X_ip) - for( tI in seq_along(Year_Set) ){ - indices = max(indices) + 1:nrow(latlon_g) - if( max(indices)>nrow(X) ) stop("Check problem in `make_covariates`") - X_gpt = abind::abind( X_gpt, X[ indices, , drop=FALSE ], along=3 ) - } - X_gtp = aperm( X_gpt, perm=c(1,3,2) ) + # Check for obvious problems + if( any(is.na(X_itp)) ) stop("Problem with `X_itp` in `make_covariates(.)") if( any(is.na(X_gtp)) ) stop("Problem with `X_gtp` in `make_covariates(.)") # warnings @@ -113,7 +106,7 @@ make_covariates = function( formula, covariate_data, Year_i, spatial_list, extra } # return stuff - Return = list( "X_gtp"=X_gtp, "X_itp"=X_itp, "covariate_names"=covariate_names, + Return = list( "X_gtp"=X_gtp, "X_itp"=X_itp, # "covariate_names"=covariate_names, "coefficient_names"=coefficient_names ) return( Return ) } From a518a2db49048d035748b4bf93c5954739a9a220 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Tue, 22 Sep 2020 09:23:38 -0700 Subject: [PATCH 14/15] fix bug in revised make_covariates(.) --- R/make_covariates.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/make_covariates.R b/R/make_covariates.R index b2df52e..95cf933 100644 --- a/R/make_covariates.R +++ b/R/make_covariates.R @@ -41,7 +41,7 @@ make_covariates = function( formula, covariate_data, Year_i, spatial_list, extra covariate_df = covariate_data[ which(!is.na(covariate_data[,'Year'])), ] for( tI in seq_along(Year_Set) ){ newrows = covariate_data[ which(is.na(covariate_data[,'Year'])), ] - newrows[,"Year"] = Year_Set[tI] + newrows[,"Year"] = rep( Year_Set[tI], nrow(newrows) ) covariate_df = rbind( covariate_df, newrows ) } From e6022c115d07fd1e300f90f0d3bdfda21305eac1 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Tue, 22 Sep 2020 13:24:15 -0700 Subject: [PATCH 15/15] Update DESCRIPTION --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5d09ee1..a286f91 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: FishStatsUtils Type: Package Title: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox Version: 2.8.0 -Date: 2020-08-07 +Date: 2020-09-22 Authors@R: c(person(given = "James", family = "Thorson",