diff --git a/NAMESPACE b/NAMESPACE index 90cea17..c451746 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,6 +19,7 @@ export(get_surveyidx_aic) export(get_surveyidx_resid) export(get_surveyidx_sim) export(get_surveyidx_stratmean) +export(plot_simulation_list) export(plot_surveyidx) export(plot_surveyidx_grid) export(qres_tweedie) diff --git a/R/functions.R b/R/functions.R index 15ece16..15d71c6 100644 --- a/R/functions.R +++ b/R/functions.R @@ -1127,6 +1127,97 @@ get_surveyidx_sim <- function( prob0=out_mu0)) } + +##' Plot survey index list (e.g. retrospective analysis) +##' +##' @title Plot survey index list (e.g. retrospective analysis) +##' @param x (named) list of "surveyIdx" objects for example from "retro.surveyIdx" or "leaveout.surveyIdx" +##' @param base Either index of x that should considered the "base run" (integer), OR object of class "surveyIdx". Confidence bounds will be shown for this model only. +##' @param rescale Should indices be rescaled to have mean 1 (over the set of intersecting years)? Default: FALSE +##' @param lwd line width argument to plot +##' @param main if not NULL override main plotting default title of "Age group a" +##' @param allCI show 95\% confidence lines for all indices? Default FALSE. +##' @param includeCI Show confidence intervals? Default TRUE. +##' @param ylim Y axis range. If NULL (default) then determine automatically. +##' @return nothing +##' @export +plot_simulation_list<-function(x, base=1, rescale=FALSE,lwd=1.5,main=NULL,allCI=FALSE,includeCI=TRUE,ylim=NULL){ + if(class(base)=="surveyIdx"){ + x = c( list(base), x) + base = 1 + } + stopifnot(is.numeric(base)) + nx = length(x) + mainwasnull = is.null(main) + n = ncol(x[[base]]$idx) + if(n>1){ + op <- par(mfrow=n2mfrow(n)) + on.exit(par(op)) + } + + cols = rainbow(nx) + if(nx==2) cols = 2:3 + cols[base] = "black" + allyears = lapply(x, function(x) rownames(x$idx)) + rsidx = 1:nrow(x[[nx]]$idx) + if(rescale){ + commonyears = allyears[[1]] + if(nx>1){ + for(i in 2:nx){ + commonyears = intersect(commonyears,allyears[[i]]) + } + if(length(commonyears)==0) stop("rescaling not possible because the set of common years is empty") + } + } + + ss <- ssbase <- 1 + + for(aa in 1:n){ + + rangevec = x[[1]]$idx[,aa] + for(xx in 2:nx) rangevec = c(rangevec,x[[xx]]$idx[,aa]) + if(includeCI){ + for(xx in 1:nx) rangevec = c(rangevec,x[[xx]]$lo[,aa],x[[xx]]$up[,aa]) + } + + yl = range(rangevec) + if(rescale){ + rsidx = which(rownames(x[[base]]$idx) %in% commonyears ) + ssbase = mean( x[[base]]$idx[rsidx,aa], na.rm=TRUE) + yl = yl/ssbase + } + if(!is.null(ylim)) yl = ylim + + if(mainwasnull) main <- paste("Age group", colnames(x[[base]]$idx)[aa]) + y = as.numeric(rownames(x[[base]]$idx)) + plot(y,x[[base]]$idx[,aa]/ssbase,type="b",ylim=yl,main=main,xlab="Year",ylab="Index") + + if(includeCI) + polygon(c(y, rev(y)), c(x[[base]]$lo[,aa], rev(x[[base]]$up[,aa]))/ssbase, col = "lightgrey", border = NA) + + for(i in 1:length(x)){ + y = as.numeric(rownames(x[[i]]$idx)) + if(rescale){ + rsidx = which(rownames(x[[i]]$idx) %in% commonyears ) + ss = mean( x[[i]]$idx[rsidx,aa], na.rm=TRUE) + } + lines(y,x[[i]]$idx[,aa]/ss,col=cols[i],type="b", lwd=lwd) + + if(includeCI && allCI && i!=base){ + lines(y,x[[i]]$lo[,aa]/ss,col=cols[i],lwd=lwd*0.6,lty=2) + lines(y,x[[i]]$up[,aa]/ss,col=cols[i],lwd=lwd*0.6,lty=2) + } + + } + y = as.numeric(rownames(x[[base]]$idx)) + lines(y,x[[base]]$idx[,aa]/ssbase,type="b",lwd=lwd) + + } + if(!is.null(names(x))){ + legend("topleft",legend=names(x),col=cols,lty=1,lwd=lwd,pch=1) + } +} + # Randomized quantile residuals ------------------------------------------------ #' Randomized quantile residuals for class 'surveyIndex' diff --git a/README.Rmd b/README.Rmd index 9df471e..caa183c 100644 --- a/README.Rmd +++ b/README.Rmd @@ -6,7 +6,6 @@ output: ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) -link_foss <- "https://github.com/EmilyMarkowitz-NOAA/surveyIndex" ``` @@ -36,13 +35,19 @@ cat( ) ``` -> Code is still in development +> Code is still in development at https://github.com/EmilyMarkowitz-NOAA/surveyIndex -Code was originally developed by: **Casper W. Berg** (@casperwberg) +*Code was originally developed by:* -> Berg et al. (2014): "Evaluation of alternative age-based methods for estimating relative abundance from survey data in relation to assessment models", Fisheries Research 151(2014) 91-99. +**Casper W. Berg** (@casperwberg) -And then modified and adapted for the AFSC by: +National Institute of Aquatic Resources, + +Technical University of Denmark + +[**Berg et al. (2014): "Evaluation of alternative age-based methods for estimating relative abundance from survey data in relation to assessment models", Fisheries Research 151(2014) 91-99.**](https://doi.org/10.1016/j.fishres.2013.10.005) + +*And then modified and adapted for the AFSC by:* **Emily Markowitz** (@EmilyMarkowitz-noaa; Emily.Markowitz AT noaa.gov) diff --git a/README.md b/README.md index a22aad6..d9a4ec2 100644 --- a/README.md +++ b/README.md @@ -7,15 +7,23 @@ data. [![](https://img.shields.io/github/last-commit/EmilyMarkowitz-NOAA/surveyIndex.svg)](https://github.com/EmilyMarkowitz-NOAA/surveyIndex/commits/main) -> Code is still in development +> Code is still in development at +> -Code was originally developed by: **Casper W. Berg** (@casperwberg) +*Code was originally developed by:* -> Berg et al. (2014): “Evaluation of alternative age-based methods for -> estimating relative abundance from survey data in relation to -> assessment models”, Fisheries Research 151(2014) 91-99. +**Casper W. Berg** (@casperwberg) -And then modified and adapted for the AFSC by: +National Institute of Aquatic Resources, + +Technical University of Denmark + +[**Berg et al. (2014): “Evaluation of alternative age-based methods for +estimating relative abundance from survey data in relation to assessment +models”, Fisheries Research 151(2014) +91-99.**](https://doi.org/10.1016/j.fishres.2013.10.005) + +*And then modified and adapted for the AFSC by:* **Emily Markowitz** (@EmilyMarkowitz-noaa; Emily.Markowitz AT noaa.gov) @@ -65,20 +73,20 @@ This package was last produced using: ## [1] parallel stats graphics grDevices utils datasets methods base ## ## other attached packages: - ## [1] surveyIndex_0.1.0 marmap_1.0.6 RANN_2.6.1 maptools_1.1-4 sp_1.5-0 mapdata_2.3.0 - ## [7] maps_3.4.0 mgcv_1.8-40 nlme_3.1-157 DATRAS_1.01 RODBC_1.3-19 roxygen2_7.2.1 - ## [13] devtools_2.4.4 knitr_1.40 badger_0.2.1 pkgdown_2.0.6 usethis_2.1.6 here_1.0.1 + ## [1] surveyIndex_0.1.0 pkgdown_2.0.6.9000 knitr_1.40 badger_0.2.1 marmap_1.0.6 RANN_2.6.1 + ## [7] maptools_1.1-4 sp_1.5-0 mapdata_2.3.0 maps_3.4.0 mgcv_1.8-40 nlme_3.1-157 + ## [13] DATRAS_1.01 RODBC_1.3-19 roxygen2_7.2.1 devtools_2.4.4 usethis_2.1.6 here_1.0.1 ## ## loaded via a namespace (and not attached): - ## [1] colorspace_2.0-3 ellipsis_0.3.2 gitcreds_0.1.2 rprojroot_2.0.3 fs_1.5.2 + ## [1] colorspace_2.0-3 gitcreds_0.1.2 ellipsis_0.3.2 rprojroot_2.0.3 fs_1.5.2 ## [6] rstudioapi_0.14 remotes_2.4.2 gh_1.3.1 bit64_4.0.5 fansi_1.0.3 ## [11] xml2_1.3.3 codetools_0.2-18 splines_4.2.1 ncdf4_1.19 cachem_1.0.6 - ## [16] pkgload_1.3.0 jsonlite_1.8.0 icesDatras_1.4.0 shiny_1.7.2 BiocManager_1.30.18 - ## [21] compiler_4.2.1 httr_1.4.4 rvcheck_0.2.1 assertthat_0.2.1 Matrix_1.5-0 + ## [16] pkgload_1.3.0 jsonlite_1.8.0 icesDatras_1.4.0 shiny_1.7.2 httr_1.4.4 + ## [21] BiocManager_1.30.18 compiler_4.2.1 rvcheck_0.2.1 assertthat_0.2.1 Matrix_1.4-1 ## [26] fastmap_1.1.0 cli_3.4.0 later_1.3.0 htmltools_0.5.3 prettyunits_1.1.1 ## [31] tools_4.2.1 gtable_0.3.1 glue_1.6.2 reshape2_1.4.4 dplyr_1.0.10 ## [36] Rcpp_1.0.9 raster_3.5-29 vctrs_0.4.1 xfun_0.32 stringr_1.4.1 - ## [41] brio_1.1.3 ps_1.7.1 testthat_3.1.4 mime_0.12 miniUI_0.1.1.1 + ## [41] ps_1.7.1 brio_1.1.3 testthat_3.1.4 mime_0.12 miniUI_0.1.1.1 ## [46] lifecycle_1.0.2 dlstats_0.1.5 sys_3.4 terra_1.6-17 MASS_7.3-57 ## [51] scales_1.2.1 promises_1.2.0.1 credentials_1.3.2 RColorBrewer_1.1-3 gert_1.8.0 ## [56] yaml_2.3.5 curl_4.3.2 memoise_2.0.1 ggplot2_3.3.6 yulab.utils_0.0.5 diff --git a/docs/404.html b/docs/404.html index 053d505..a69fb64 100644 --- a/docs/404.html +++ b/docs/404.html @@ -96,7 +96,7 @@

Page not found (404)

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/LICENSE-text.html b/docs/LICENSE-text.html index 6742769..9aea796 100644 --- a/docs/LICENSE-text.html +++ b/docs/LICENSE-text.html @@ -744,7 +744,7 @@

License

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/articles/A-data-prep.html b/docs/articles/A-data-prep.html index b14ea8a..febf395 100644 --- a/docs/articles/A-data-prep.html +++ b/docs/articles/A-data-prep.html @@ -227,8 +227,8 @@

2. Prepare the data for s Quarter = "2") %>% dplyr::mutate(across((c("Year", "Ship", "COMMON_NAME")), as.factor)) %>% dplyr::select(wCPUE, GEAR_TEMPERATURE, BOTTOM_DEPTH, COMMON_NAME, EFFORT, - Year, Ship, Lon, Lat, lat, lon, sx, sy, ctime, - TimeShotHour, timeOfYear, Gear, Quarter, HaulDur, hauljoin) + Year, Ship, Lon, Lat, lat, lon, sx, sy, + ctime, TimeShotHour, timeOfYear, Gear, Quarter, HaulDur, hauljoin) head(dat_wrangled) #> # A tibble: 6 × 20 @@ -291,7 +291,7 @@

4. Prepare covariate data dat_cov <- surveyIndex::pred_grid_ebs %>% dplyr::select(-Shape_Area) %>% - dplyr::mutate( + dplyr::mutate( sx = ((lon - mean(lon, na.rm = TRUE))/1000), sy = ((lat - mean(lat, na.rm = TRUE))/1000)) @@ -335,7 +335,7 @@

4a. Data that varies over onl # to make future runs of this faster: save(extrap_data0, extrap_data, file = paste0("../inst/VigA_bottom_depth_raster_", - min(YEARS),"-",max(YEARS), ".rdata")) + min(YEARS),"-",max(YEARS), ".rdata")) 
 # Just so we can see what we are looking at:
 plot(extrap_data0, main = "Interpolated Bottom Depths")
@@ -343,7 +343,7 @@

4a. Data that varies over onl 
 
 dat_cov <- cbind.data.frame(dat_cov, 
-                      "BOTTOM_DEPTH" = extrap_data$var1.pred) %>%
+                            "BOTTOM_DEPTH" = extrap_data$var1.pred) %>%
   stats::na.omit()
 
 head(dat_cov)
@@ -379,12 +379,12 @@ 
4b. Data tha
 }
 
 extrap_data0 <- coldpool::ebs_bottom_temperature[[tmp]] %>% 
-    as(., Class = "SpatialPointsDataFrame") %>%
-      sf::st_as_sf() %>%
+  as(., Class = "SpatialPointsDataFrame") %>%
+  sf::st_as_sf() %>%
   sf::st_transform(crs = crs_latlon)  %>%
   stars::st_rasterize() %>% 
-    stars::st_extract(x = .,
-                      at = as.matrix(dat_cov[,c("lon", "lat")]))
+  stars::st_extract(x = .,
+                    at = as.matrix(dat_cov[,c("lon", "lat")]))
 
 names(extrap_data0) <- paste0("GEAR_TEMPERATURE", YEARS)
 
@@ -424,59 +424,60 @@ 
5b. Catch DataNow, we need to fill in the data with the zeros!

 
 
-  # Identify vars that will be used --------------------------------------------
+# Identify vars that will be used --------------------------------------------
 
-  varsbyyr <- unique( # c("GEAR_TEMPERATURE", "cpi")
-    gsub(pattern = "[0-9]+", 
-         replacement = "", 
-         x = names(dat_cov)[grepl(names(dat_cov), 
-                                   pattern = "[0-9]+")]))
-  
-  vars <- unique( # c("BOTTOM_DEPTH")
-    names(dat_cov)[!grepl(names(dat_cov), 
-                           pattern = "[0-9]+")])
-  vars <- vars[!(vars %in% c("LONG", "LAT", "lon", "lat", "sx", "sy"))]
-  
-  ## Fill catch data with zeros ---------------------------------------------------------
-  
-  data_hauls <- dat_wrangled %>% 
-    dplyr::select(Year, sx, sy, 
-                  dplyr::all_of(varsbyyr), dplyr::all_of(vars), 
-                  Ship, hauljoin, 
-                  lat, lon, 
-                  ctime, TimeShotHour, timeOfYear, Gear, Quarter, 
-                  EFFORT, HaulDur)  %>% 
-    # dplyr::filter(!is.na(GEAR_TEMPERATURE)) %>% 
-    na.omit() %>%
-    dplyr::distinct()
-  
-  data_catch <- dat_wrangled %>% 
-    dplyr::select(COMMON_NAME, wCPUE, hauljoin)
-  
-  dat_catch_haul <- dplyr::left_join(x = data_hauls, 
-                        y = data_catch, 
-                        by = c("hauljoin")) %>% 
-    dplyr::mutate(wCPUE = ifelse(is.na(wCPUE), 0, wCPUE))
-  
-  head(dat_catch_haul)
-#> # A tibble: 6 × 18
-#>   Year       sx       sy GEAR_TEMPE…¹ BOTTO…² Ship  haulj…³    lat     lon ctime
+varsbyyr <- unique( # c("GEAR_TEMPERATURE", "cpi")
+  gsub(pattern = "[0-9]+", 
+       replacement = "", 
+       x = names(dat_cov)[grepl(names(dat_cov), 
+                                pattern = "[0-9]+")]))
+
+vars <- unique( # c("BOTTOM_DEPTH")
+  names(dat_cov)[!grepl(names(dat_cov), 
+                        pattern = "[0-9]+")])
+vars <- vars[!(vars %in% c("LONG", "LAT", "lon", "lat", "sx", "sy"))]
+
+## Fill catch data with zeros ---------------------------------------------------------
+
+data_hauls <- dat_wrangled %>% 
+  dplyr::select(Year, sx, sy, 
+                dplyr::all_of(varsbyyr), dplyr::all_of(vars), 
+                Ship, hauljoin, 
+                lat, lon, Lat, Lon, 
+                ctime, TimeShotHour, timeOfYear, Gear, Quarter, 
+                EFFORT, HaulDur)  %>% 
+  # dplyr::filter(!is.na(GEAR_TEMPERATURE)) %>% 
+  na.omit() %>%
+  dplyr::distinct()
+
+data_catch <- dat_wrangled %>% 
+  dplyr::select(COMMON_NAME, wCPUE, hauljoin)
+
+dat_catch_haul <- dplyr::left_join(x = data_hauls, 
+                                   y = data_catch, 
+                                   by = c("hauljoin")) %>% 
+  dplyr::mutate(wCPUE = ifelse(is.na(wCPUE), 0, wCPUE))
+
+head(dat_catch_haul)
+#> # A tibble: 6 × 20
+#>   Year       sx       sy GEAR_TEMPE…¹ BOTTO…² Ship  haulj…³    lat     lon   Lat
 #>   <fct>   <dbl>    <dbl>        <dbl>   <dbl> <fct> <chr>    <dbl>   <dbl> <dbl>
-#> 1 2016  0.00647 -0.00183          6.1      54 94    10_E-1… 7.26e5 -4.32e5  2016
-#> 2 2016  0.00647 -0.00183          6.1      54 94    10_E-1… 7.26e5 -4.32e5  2016
-#> 3 2016  0.00647 -0.00183          6.1      54 94    10_E-1… 7.26e5 -4.32e5  2016
-#> 4 2016  0.00518 -0.00217          5.2      73 162   31_D-1… 6.98e5 -5.16e5  2016
-#> 5 2016  0.00518 -0.00217          5.2      73 162   31_D-1… 6.98e5 -5.16e5  2016
-#> 6 2016  0.00518 -0.00217          5.2      73 162   31_D-1… 6.98e5 -5.16e5  2016
-#> # … with 8 more variables: TimeShotHour <dbl>, timeOfYear <dbl>, Gear <chr>,
-#> #   Quarter <chr>, EFFORT <dbl>, HaulDur <dbl>, COMMON_NAME <fct>, wCPUE <dbl>,
-#> #   and abbreviated variable names ¹​GEAR_TEMPERATURE, ²​BOTTOM_DEPTH, ³​hauljoin

+#> 1 2016  0.00647 -0.00183          6.1      54 94    10_E-1… 7.26e5 -4.32e5  56.3
+#> 2 2016  0.00647 -0.00183          6.1      54 94    10_E-1… 7.26e5 -4.32e5  56.3
+#> 3 2016  0.00647 -0.00183          6.1      54 94    10_E-1… 7.26e5 -4.32e5  56.3
+#> 4 2016  0.00518 -0.00217          5.2      73 162   31_D-1… 6.98e5 -5.16e5  56.0
+#> 5 2016  0.00518 -0.00217          5.2      73 162   31_D-1… 6.98e5 -5.16e5  56.0
+#> 6 2016  0.00518 -0.00217          5.2      73 162   31_D-1… 6.98e5 -5.16e5  56.0
+#> # … with 10 more variables: Lon <dbl>, ctime <dbl>, TimeShotHour <dbl>,
+#> #   timeOfYear <dbl>, Gear <chr>, Quarter <chr>, EFFORT <dbl>, HaulDur <dbl>,
+#> #   COMMON_NAME <fct>, wCPUE <dbl>, and abbreviated variable names
+#> #   ¹​GEAR_TEMPERATURE, ²​BOTTOM_DEPTH, ³​hauljoin
- allpd <- lapply(YEARS, FUN = surveyIndex::get_prediction_grid, x = dat_cov, 
-                 vars = vars, varsbyyr = varsbyyr)
-  names(allpd) <- as.character(YEARS)
-  
-  head(allpd[1][[1]])
+allpd <- lapply(YEARS, FUN = surveyIndex::get_prediction_grid, x = dat_cov, 
+                vars = vars, varsbyyr = varsbyyr)
+names(allpd) <- as.character(YEARS)
+
+head(allpd[1][[1]])
 #>          lon      lat           sx          sy BOTTOM_DEPTH GEAR_TEMPERATURE
 #> 7  -176.2507 62.09301 -0.007524943 0.003800355     92.30377       0.06603578
 #> 11 -175.9816 62.13632 -0.007255819 0.003843669     91.79176      -0.19408995
@@ -497,12 +498,12 @@ 
5a. Covariate Data
 
-  ## split data by species, make into DATRASraw + Nage matrix
-  ds <- split(dat_catch_haul,dat_catch_haul$COMMON_NAME)
-  ds <- lapply(ds, surveyIndex::get_datrasraw)
-  ## OBS, response is added here in "Nage" matrix -- use wCPUE
-  ds <- lapply(ds,function(x) { x[[2]]$Nage <- matrix(x$wCPUE,ncol=1); colnames(x[[2]]$Nage)<-1; x } )
-  
+## split data by species, make into DATRASraw + Nage matrix
+ds <- split(dat_catch_haul,dat_catch_haul$COMMON_NAME)
+ds <- lapply(ds, surveyIndex::get_datrasraw)
+## OBS, response is added here in "Nage" matrix -- use wCPUE
+ds <- lapply(ds,function(x) { x[[2]]$Nage <- matrix(x$wCPUE,ncol=1); colnames(x[[2]]$Nage)<-1; x } )
+
 ds
 #> $`red king crab`
 #> Object of class 'DATRASraw'
@@ -547,10 +548,10 @@ 
6. Formulas"fm_1_s_t_st" = "Year +
     s(sx,sy,bs=c('ts'),k=376) + 
     s(sx,sy,bs=c('ts'),k=10,by=Year)",
-
+  
   # Mdoel with simple covariates
   "fm_2_cov" = 
-"s(BOTTOM_DEPTH,bs='ts',k=10) +
+    "s(BOTTOM_DEPTH,bs='ts',k=10) +
 s(log(GEAR_TEMPERATURE+3),bs='ts',k=10)"
 )
@@ -559,16 +560,16 @@

7. Fit the Model

Here are all of the models we want to try fitting:

-  comb <- tidyr::crossing(
-    "SPECIES" = SPECIES, 
-    "fm_name" = gsub(pattern = " ", replacement = "_", x = names(fm))) %>% 
-    dplyr::left_join(
-      x = ., 
-      y = data.frame("fm" = gsub(pattern = "\n", replacement = "", 
-                                 x = unlist(fm), fixed = TRUE), 
-                     "fm_name" = gsub(pattern = " ", replacement = "_", 
-                                      x = names(fm))), 
-                     by = "fm_name")
+comb <- tidyr::crossing(
+  "SPECIES" = SPECIES, 
+  "fm_name" = gsub(pattern = " ", replacement = "_", x = names(fm))) %>% 
+  dplyr::left_join(
+    x = ., 
+    y = data.frame("fm" = gsub(pattern = "\n", replacement = "", 
+                               x = unlist(fm), fixed = TRUE), 
+                   "fm_name" = gsub(pattern = " ", replacement = "_", 
+                                    x = names(fm))), 
+    by = "fm_name")
 
 comb
 #> # A tibble: 6 × 3
@@ -581,36 +582,36 @@ 
7. Fit the Model#> 5 yellowfin sole  fm_1_s_t_st Year +    s(sx,sy,bs=c('ts'),k=376) +     s(sx,sy…
 #> 6 yellowfin sole  fm_2_cov    s(BOTTOM_DEPTH,bs='ts',k=10) +s(log(GEAR_TEMPERAT…
-  
-  models <- fittimes <- list()
-
-  for(i in 1:nrow(comb)){
-    cat("Fitting ",comb$SPECIES[i],"\n", comb$fm_name[i], ": ", comb$fm[i], "\n")
-
-    temp <- paste0(comb$SPECIES[i], " ", comb$fm_name[i])
-
-    fittimes[[ temp ]] <-
-      system.time ( models[[ temp ]] <-
-                     surveyIndex::get_surveyidx(
-                       x = ds[[ comb$SPECIES[i] ]],
-                                   ages = 1,
-                                   myids = NULL,
-                                   predD = allpd,
-                                   cutOff = 0,
-                                   fam = "Tweedie",
-                                   modelP = comb$fm[i],
-                                   gamma = 1,
-                                   control = list(trace = TRUE,
-                                                  maxit = 20))  )
+
+models <- fittimes <- list()
 
-  }
+for(i in 1:nrow(comb)){
+  cat("Fitting ",comb$SPECIES[i],"\n", comb$fm_name[i], ": ", comb$fm[i], "\n")
+  
+  temp <- paste0(comb$SPECIES[i], " ", comb$fm_name[i])
+  
+  fittimes[[ temp ]] <-
+    system.time ( models[[ temp ]] <-
+                    surveyIndex::get_surveyidx(
+                      x = ds[[ comb$SPECIES[i] ]],
+                      ages = 1,
+                      myids = NULL,
+                      predD = allpd,
+                      cutOff = 0,
+                      fam = "Tweedie",
+                      modelP = comb$fm[i],
+                      gamma = 1,
+                      control = list(trace = TRUE,
+                                     maxit = 20))  )
+  
+}
 
-  save(models, fittimes, file = paste0("../inst/VigA_model_fits.Rdata"))
+save(models, fittimes, file = paste0("../inst/VigA_model_fits.Rdata")) 
-  ## Check basis dimensions splines (spatial resolution)
-  # sink(paste0(".", dir_out, "gamcheck.txt"))
+## Check basis dimensions splines (spatial resolution)
+# sink(paste0(".", dir_out, "gamcheck.txt"))
 par(mfrow = c(2,2))
-  lapply(models,function(x) gam.check(x$pModels[[1]]))
+lapply(models,function(x) gam.check(x$pModels[[1]])) 

#> 
 #> Method: ML   Optimizer: outer newton
@@ -645,7 +646,7 @@ 
7. Fit the Model#> 
 #>                                 k'   edf k-index p-value    
 #> s(BOTTOM_DEPTH)              9.000 4.388    0.76  <2e-16 ***
-#> s(log(GEAR_TEMPERATURE + 3)) 9.000 0.944    0.77  <2e-16 ***
+#> s(log(GEAR_TEMPERATURE + 3)) 9.000 0.944    0.77    0.01 ** 
 #> ---
 #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

@@ -660,16 +661,14 @@

7. Fit the Model#> Basis dimension (k) checking results. Low p-value (k-index<1) may #> indicate that k is too low, especially if edf is close to k'. #> -#> k' edf k-index p-value -#> s(sx,sy) 3.75e+02 1.26e+02 0.87 0.165 -#> s(sx,sy):Year2015 9.00e+00 6.23e+00 0.87 0.090 . -#> s(sx,sy):Year2016 9.00e+00 1.29e-02 0.87 0.145 -#> s(sx,sy):Year2017 9.00e+00 1.94e+00 0.87 0.155 -#> s(sx,sy):Year2018 9.00e+00 8.12e+00 0.87 0.140 -#> s(sx,sy):Year2019 9.00e+00 2.01e+00 0.87 0.085 . -#> s(sx,sy):Year2021 9.00e+00 6.41e-03 0.87 0.115 -#> --- -#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +#> k' edf k-index p-value +#> s(sx,sy) 3.75e+02 1.26e+02 0.87 0.11 +#> s(sx,sy):Year2015 9.00e+00 6.23e+00 0.87 0.14 +#> s(sx,sy):Year2016 9.00e+00 1.29e-02 0.87 0.16 +#> s(sx,sy):Year2017 9.00e+00 1.94e+00 0.87 0.14 +#> s(sx,sy):Year2018 9.00e+00 8.12e+00 0.87 0.13 +#> s(sx,sy):Year2019 9.00e+00 2.01e+00 0.87 0.11 +#> s(sx,sy):Year2021 9.00e+00 6.41e-03 0.87 0.17 

#> 
 #> Method: ML   Optimizer: outer newton
@@ -702,10 +701,10 @@ 
7. Fit the Model#>                         k'      edf k-index p-value
 #> s(sx,sy)          3.75e+02 1.52e+02       1    0.99
 #> s(sx,sy):Year2015 9.00e+00 1.80e+00       1    0.99
-#> s(sx,sy):Year2016 9.00e+00 1.15e-03       1    0.99
-#> s(sx,sy):Year2017 9.00e+00 2.08e-03       1    0.99
-#> s(sx,sy):Year2018 9.00e+00 1.87e+00       1    1.00
-#> s(sx,sy):Year2019 9.00e+00 6.64e+00       1    0.99
+#> s(sx,sy):Year2016 9.00e+00 1.15e-03       1    1.00
+#> s(sx,sy):Year2017 9.00e+00 2.08e-03       1    0.98
+#> s(sx,sy):Year2018 9.00e+00 1.87e+00       1    0.99
+#> s(sx,sy):Year2019 9.00e+00 6.64e+00       1    0.98
 #> s(sx,sy):Year2021 9.00e+00 1.95e+00       1    0.99

#> 
@@ -752,11 +751,11 @@ 
7. Fit the Model#> $`yellowfin sole fm_2_cov`
 #> $`yellowfin sole fm_2_cov`$mfrow
 #> [1] 2 2
-  # sink()
+# sink() 
-  ## Model summaries
-  # sink(paste0(".", dir_out, "summaries.txt"))
-  lapply(models,function(x) summary(x$pModels[[1]]))
+## Model summaries
+# sink(paste0(".", dir_out, "summaries.txt"))
+lapply(models,function(x) summary(x$pModels[[1]]))
 #> $`red king crab fm_1_s_t_st`
 #> 
 #> Family: Tweedie(p=1.99) 
@@ -936,15 +935,15 @@ 
7. Fit the Model#> 
 #> R-sq.(adj) =  0.152   Deviance explained = 35.1%
 #> -ML = 7295.9  Scale est. = 1.4146    n = 1497
-  # sink()
+# sink()
 
-  surveyIndex::get_surveyidx_aic(x = models)
+surveyIndex::get_surveyidx_aic(x = models)
 #> numeric(0)
-  temp <- sapply(models, `[`, "pModels")
-  mods <- sapply(temp, `[`, 1)
+temp <- sapply(models, `[`, "pModels")
+mods <- sapply(temp, `[`, 1)
 
-  lapply(X = mods, FUN = AIC)
+lapply(X = mods, FUN = AIC)
 #> $`red king crab fm_1_s_t_st.pModels`
 #> [1] 2296.156
 #> 
@@ -967,40 +966,124 @@ 
7. Fit the Model8. Indicies of Abundance
 

 
-dat <- data.frame()
+
+dat <- data.frame()
 for (i in 1:length(models)){
   temp <- models[[i]]
-  dat <- dplyr::bind_rows(dat, 
-                data.frame(idx = temp$idx[,1], 
-                    Year = rownames(temp$idx), 
-                    group = names(models)[i]))
+  dat0 <- data.frame(idx = temp$idx[,1], 
+                     lo = temp$lo[,1], 
+                     up = temp$up[,1],
+                     Year = rownames(temp$idx), 
+                     group = names(models)[i], 
+                     formula = paste0("cpue_kgha ~ ", 
+                                      as.character(temp$pModels[[1]]$formula)[[3]]))
+  
+  dat <- dplyr::bind_rows(dat, dat0) 
 }
 
-dat$facet_group <- sapply(X = strsplit(x = dat$group, split = " "), 
-                          `[`, 1)
+dat$facet_group <- paste0(sapply(X = strsplit(x = dat$group, split = " fm"), `[`, 1))
+# dat$model <- paste0(sapply(X = strsplit(x = dat$group, split = " fm"), `[`, 2))
 
-dat$idx[dat$Year == 2020] <- NA
+dat[dat$Year == 2020, c("idx", "up", "lo")] <- NA
 
 ggplot2::ggplot(data = dat, 
-            mapping = aes(x = Year, 
-                          y = idx, 
-                          group = group, 
-                          color = group)) +
+                mapping = aes(x = Year, 
+                              y = idx, 
+                              group = formula, 
+                              color = formula)) +
   geom_line(size = 1.5) + 
-  geom_point(size = 2) +
+  geom_point(size = 2)  + 
+  geom_ribbon(aes(ymin = lo, ymax = up, fill = formula), 
+              alpha=0.1, 
+              linetype="dashed",
+              color="grey") + 
+  ggtitle("Annual Index Model Results") +
   facet_wrap(vars(facet_group), scales = "free", ncol = 1) +
-  theme(legend.position = "bottom")

+  theme(legend.position = "bottom", 
+        legend.direction = "vertical")

-

8. Predict and plot +

9. Predict and plot 

 
-# predict.gam(object = models$`red king crab fm_2_cov`, 
-#             newdata = dat_catch_haul %>% 
-#               dplyr::filter(Year == 2021) %>%
-#               dplyr::select(Year, sx, sy, GEAR_TEMPERATURE, BOTTOM_DEPTH))
+dat_pred <- dat_catch_haul %>% + dplyr::select(Year, sx, sy, Lon, Lat, GEAR_TEMPERATURE, BOTTOM_DEPTH) + +dat <- data.frame() +for (i in 1:length(models)) { + temp <- models[[i]] + dat0 <- data.frame(idx = + predict.gam( + object = temp$pModels[[1]], + newdata = dat_pred), + group = names(models)[i], + formula = paste0("cpue_kgha ~ ", + as.character(temp$pModels[[1]]$formula)[[3]]) + ) + dat00 <- dplyr::bind_cols(dat0, dat_pred) + dat <- dplyr::bind_rows(dat, dat00) + +# dat_r <- raster::rasterFromXYZ(xyz = dat00[,c("lon", "lat", "idx")]) + +} + +dat$facet_group <- paste0(sapply(X = strsplit(x = dat$group, split = " fm"), `[`, 1)) + +for (i in 1:length(unique(dat$facet_group))){ + + ggplot2::ggplot(data = dat %>% + dplyr::filter(facet_group == unique(dat$facet_group)[i]), + mapping = aes(x = Lon, + y = Lat, + group = group, + color = idx)) + + geom_point() + + ggtitle(paste0("Annual Index Model Results for ", unique(dat$facet_group)[i])) + + facet_grid(cols = vars(group), + rows = vars(Year)) + + theme_bw() + +}
+ +
+

10. Simulations +

+
+
+sims <- fittimes_sims <- list()
+for(i in 1:nrow(comb)){
+  
+  cat("Simulating ",comb$SPECIES[i],"\n", comb$fm_name[i], ": ", comb$fm[i], "\n")
+  
+  temp <- paste0(comb$SPECIES[i], " ", comb$fm_name[i])
+  
+  fittimes[[ temp ]] <-
+    system.time ( sims[[ temp ]] <-
+                    surveyIndex::get_surveyidx_sim(
+                      model = models[[i]], 
+                      d = ds[[ comb$SPECIES[i] ]]) )
+}
+#> Simulating  red king crab 
+#>  fm_1_s_t_st :  Year +    s(sx,sy,bs=c('ts'),k=376) +     s(sx,sy,bs=c('ts'),k=10,by=Year) 
+#> Simulating  red king crab 
+#>  fm_2_cov :  s(BOTTOM_DEPTH,bs='ts',k=10) +s(log(GEAR_TEMPERATURE+3),bs='ts',k=10) 
+#> Simulating  walleye pollock 
+#>  fm_1_s_t_st :  Year +    s(sx,sy,bs=c('ts'),k=376) +     s(sx,sy,bs=c('ts'),k=10,by=Year) 
+#> Simulating  walleye pollock 
+#>  fm_2_cov :  s(BOTTOM_DEPTH,bs='ts',k=10) +s(log(GEAR_TEMPERATURE+3),bs='ts',k=10) 
+#> Simulating  yellowfin sole 
+#>  fm_1_s_t_st :  Year +    s(sx,sy,bs=c('ts'),k=376) +     s(sx,sy,bs=c('ts'),k=10,by=Year) 
+#> Simulating  yellowfin sole 
+#>  fm_2_cov :  s(BOTTOM_DEPTH,bs='ts',k=10) +s(log(GEAR_TEMPERATURE+3),bs='ts',k=10)
+
+par(mfrow = c(2, 2)) # Create a 2 x 2 plotting matrix
+for(i in 1:nrow(comb)){
+  plot(sims[[i]]$sim, main = paste0(names(sims)[i], " sims"))
+  plot(sims[[i]]$mu[[1]], main = paste0(names(sims)[i], " mu"))
+}
+

@@ -1021,7 +1104,7 @@

8. Predict and plot

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/articles/A-data-prep_files/figure-html/indicie_abund-1.png b/docs/articles/A-data-prep_files/figure-html/indicie_abund-1.png index e7eae2b..4d0dcd3 100644 Binary files a/docs/articles/A-data-prep_files/figure-html/indicie_abund-1.png and b/docs/articles/A-data-prep_files/figure-html/indicie_abund-1.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/model_check1-1.png b/docs/articles/A-data-prep_files/figure-html/model_check1-1.png index 8155b4b..a1e2cc1 100644 Binary files a/docs/articles/A-data-prep_files/figure-html/model_check1-1.png and b/docs/articles/A-data-prep_files/figure-html/model_check1-1.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/model_check1-2.png b/docs/articles/A-data-prep_files/figure-html/model_check1-2.png index 3234042..d1e4a55 100644 Binary files a/docs/articles/A-data-prep_files/figure-html/model_check1-2.png and b/docs/articles/A-data-prep_files/figure-html/model_check1-2.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/model_check1-3.png b/docs/articles/A-data-prep_files/figure-html/model_check1-3.png index ec2380d..a1d4991 100644 Binary files a/docs/articles/A-data-prep_files/figure-html/model_check1-3.png and b/docs/articles/A-data-prep_files/figure-html/model_check1-3.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/model_check1-4.png b/docs/articles/A-data-prep_files/figure-html/model_check1-4.png index b098da2..6d243e4 100644 Binary files a/docs/articles/A-data-prep_files/figure-html/model_check1-4.png and b/docs/articles/A-data-prep_files/figure-html/model_check1-4.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/model_check1-5.png b/docs/articles/A-data-prep_files/figure-html/model_check1-5.png index 54e47f7..f0486dc 100644 Binary files a/docs/articles/A-data-prep_files/figure-html/model_check1-5.png and b/docs/articles/A-data-prep_files/figure-html/model_check1-5.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/model_check1-6.png b/docs/articles/A-data-prep_files/figure-html/model_check1-6.png index ee4e708..8f30410 100644 Binary files a/docs/articles/A-data-prep_files/figure-html/model_check1-6.png and b/docs/articles/A-data-prep_files/figure-html/model_check1-6.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/sim_gam-1.png b/docs/articles/A-data-prep_files/figure-html/sim_gam-1.png new file mode 100644 index 0000000..4667c47 Binary files /dev/null and b/docs/articles/A-data-prep_files/figure-html/sim_gam-1.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/sim_gam-2.png b/docs/articles/A-data-prep_files/figure-html/sim_gam-2.png new file mode 100644 index 0000000..68b8783 Binary files /dev/null and b/docs/articles/A-data-prep_files/figure-html/sim_gam-2.png differ diff --git a/docs/articles/A-data-prep_files/figure-html/sim_gam-3.png b/docs/articles/A-data-prep_files/figure-html/sim_gam-3.png new file mode 100644 index 0000000..b6478ea Binary files /dev/null and b/docs/articles/A-data-prep_files/figure-html/sim_gam-3.png differ diff --git a/docs/articles/B-model-tuning.html b/docs/articles/B-model-tuning.html index 33e019f..3f781a1 100644 --- a/docs/articles/B-model-tuning.html +++ b/docs/articles/B-model-tuning.html @@ -103,7 +103,7 @@

Nuances in model choices

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/articles/C-model-comparisons.html b/docs/articles/C-model-comparisons.html index 038fe81..37e3387 100644 --- a/docs/articles/C-model-comparisons.html +++ b/docs/articles/C-model-comparisons.html @@ -241,7 +241,7 @@

4. Assess the model#> indicate that k is too low, especially if edf is close to k'. #> #> k' edf k-index p-value -#> s(sx,sy) 375.00 95.70 0.94 0.57 +#> s(sx,sy) 375.00 95.70 0.94 0.54 #> s(sx,sy):Year 10.00 3.08 0.94 0.60 vis.gam(bb) @@ -307,7 +307,7 @@

4. Assess the model

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/articles/C-model-comparisons_files/figure-html/model_plots_1-1.png b/docs/articles/C-model-comparisons_files/figure-html/model_plots_1-1.png index dd4f437..4e2fca1 100644 Binary files a/docs/articles/C-model-comparisons_files/figure-html/model_plots_1-1.png and b/docs/articles/C-model-comparisons_files/figure-html/model_plots_1-1.png differ diff --git a/docs/articles/C-model-comparisons_files/figure-html/model_plots_2-1.png b/docs/articles/C-model-comparisons_files/figure-html/model_plots_2-1.png index 00b6c35..dc23777 100644 Binary files a/docs/articles/C-model-comparisons_files/figure-html/model_plots_2-1.png and b/docs/articles/C-model-comparisons_files/figure-html/model_plots_2-1.png differ diff --git a/docs/articles/index.html b/docs/articles/index.html index bf9f1b8..d4e502f 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -73,7 +73,7 @@

All vignettes

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/authors.html b/docs/authors.html index 1e6ceca..6b2c181 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -98,7 +98,7 @@

Citation

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/index.html b/docs/index.html index 426cf5d..22a0b63 100644 --- a/docs/index.html +++ b/docs/index.html @@ -81,13 +81,14 @@

surveyIndex

-

Code is still in development

+

Code is still in development at https://github.com/EmilyMarkowitz-NOAA/surveyIndex

-

Code was originally developed by: Casper W. Berg (@casperwberg)

-
-

Berg et al. (2014): “Evaluation of alternative age-based methods for estimating relative abundance from survey data in relation to assessment models”, Fisheries Research 151(2014) 91-99.

-
-

And then modified and adapted for the AFSC by:

+

Code was originally developed by:

+

Casper W. Berg (@casperwberg)

+

National Institute of Aquatic Resources,

+

Technical University of Denmark

+

Berg et al. (2014): “Evaluation of alternative age-based methods for estimating relative abundance from survey data in relation to assessment models”, Fisheries Research 151(2014) 91-99.

+

And then modified and adapted for the AFSC by:

Emily Markowitz (@EmilyMarkowitz-noaa; Emily.Markowitz AT noaa.gov)

Margaret Siple (@MargaretSiple-noaa; Margaret.Siple AT noaa.gov)

Alaska Fisheries Science Center

@@ -128,20 +129,20 @@

Metadata## [1] parallel stats graphics grDevices utils datasets methods base ## ## other attached packages: -## [1] surveyIndex_0.1.0 marmap_1.0.6 RANN_2.6.1 maptools_1.1-4 sp_1.5-0 mapdata_2.3.0 -## [7] maps_3.4.0 mgcv_1.8-40 nlme_3.1-157 DATRAS_1.01 RODBC_1.3-19 roxygen2_7.2.1 -## [13] devtools_2.4.4 knitr_1.40 badger_0.2.1 pkgdown_2.0.6 usethis_2.1.6 here_1.0.1 +## [1] surveyIndex_0.1.0 pkgdown_2.0.6.9000 knitr_1.40 badger_0.2.1 marmap_1.0.6 RANN_2.6.1 +## [7] maptools_1.1-4 sp_1.5-0 mapdata_2.3.0 maps_3.4.0 mgcv_1.8-40 nlme_3.1-157 +## [13] DATRAS_1.01 RODBC_1.3-19 roxygen2_7.2.1 devtools_2.4.4 usethis_2.1.6 here_1.0.1 ## ## loaded via a namespace (and not attached): -## [1] colorspace_2.0-3 ellipsis_0.3.2 gitcreds_0.1.2 rprojroot_2.0.3 fs_1.5.2 +## [1] colorspace_2.0-3 gitcreds_0.1.2 ellipsis_0.3.2 rprojroot_2.0.3 fs_1.5.2 ## [6] rstudioapi_0.14 remotes_2.4.2 gh_1.3.1 bit64_4.0.5 fansi_1.0.3 ## [11] xml2_1.3.3 codetools_0.2-18 splines_4.2.1 ncdf4_1.19 cachem_1.0.6 -## [16] pkgload_1.3.0 jsonlite_1.8.0 icesDatras_1.4.0 shiny_1.7.2 BiocManager_1.30.18 -## [21] compiler_4.2.1 httr_1.4.4 rvcheck_0.2.1 assertthat_0.2.1 Matrix_1.5-0 +## [16] pkgload_1.3.0 jsonlite_1.8.0 icesDatras_1.4.0 shiny_1.7.2 httr_1.4.4 +## [21] BiocManager_1.30.18 compiler_4.2.1 rvcheck_0.2.1 assertthat_0.2.1 Matrix_1.4-1 ## [26] fastmap_1.1.0 cli_3.4.0 later_1.3.0 htmltools_0.5.3 prettyunits_1.1.1 ## [31] tools_4.2.1 gtable_0.3.1 glue_1.6.2 reshape2_1.4.4 dplyr_1.0.10 ## [36] Rcpp_1.0.9 raster_3.5-29 vctrs_0.4.1 xfun_0.32 stringr_1.4.1 -## [41] brio_1.1.3 ps_1.7.1 testthat_3.1.4 mime_0.12 miniUI_0.1.1.1 +## [41] ps_1.7.1 brio_1.1.3 testthat_3.1.4 mime_0.12 miniUI_0.1.1.1 ## [46] lifecycle_1.0.2 dlstats_0.1.5 sys_3.4 terra_1.6-17 MASS_7.3-57 ## [51] scales_1.2.1 promises_1.2.0.1 credentials_1.3.2 RColorBrewer_1.1-3 gert_1.8.0 ## [56] yaml_2.3.5 curl_4.3.2 memoise_2.0.1 ggplot2_3.3.6 yulab.utils_0.0.5 @@ -209,7 +210,7 @@

Developers

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 9b03ac0..21860eb 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -1,9 +1,9 @@ pandoc: '2.18' -pkgdown: 2.0.6 -pkgdown_sha: ~ +pkgdown: 2.0.6.9000 +pkgdown_sha: 956f076f661eb2de9413f64ea24a580cf444466d articles: A-data-prep: A-data-prep.html B-model-tuning: B-model-tuning.html C-model-comparisons: C-model-comparisons.html -last_built: 2022-09-13T23:03Z +last_built: 2022-09-14T04:19Z diff --git a/docs/reference/Rplot002.png b/docs/reference/Rplot002.png index 671f1cf..e984268 100644 Binary files a/docs/reference/Rplot002.png and b/docs/reference/Rplot002.png differ diff --git a/docs/reference/anova_likelihood_ratio_test.html b/docs/reference/anova_likelihood_ratio_test.html index 29ee133..44e6660 100644 --- a/docs/reference/anova_likelihood_ratio_test.html +++ b/docs/reference/anova_likelihood_ratio_test.html @@ -125,7 +125,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/calc_distance.html b/docs/reference/calc_distance.html index 8a02fea..b58fb7c 100644 --- a/docs/reference/calc_distance.html +++ b/docs/reference/calc_distance.html @@ -117,7 +117,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/concentration_transform-2.png b/docs/reference/concentration_transform-2.png index 74f2cc0..349120c 100644 Binary files a/docs/reference/concentration_transform-2.png and b/docs/reference/concentration_transform-2.png differ diff --git a/docs/reference/concentration_transform.html b/docs/reference/concentration_transform.html index 4dfdb64..02987ac 100644 --- a/docs/reference/concentration_transform.html +++ b/docs/reference/concentration_transform.html @@ -108,7 +108,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/consistency_external.html b/docs/reference/consistency_external.html index 8424f58..ca10fa8 100644 --- a/docs/reference/consistency_external.html +++ b/docs/reference/consistency_external.html @@ -101,7 +101,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/consistency_internal.html b/docs/reference/consistency_internal.html index 9750516..2ae2708 100644 --- a/docs/reference/consistency_internal.html +++ b/docs/reference/consistency_internal.html @@ -96,7 +96,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/convert_crs.html b/docs/reference/convert_crs.html index ac260e6..f8fc73d 100644 --- a/docs/reference/convert_crs.html +++ b/docs/reference/convert_crs.html @@ -131,7 +131,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/convert_deg_rad.html b/docs/reference/convert_deg_rad.html index f6a6d15..76a08de 100644 --- a/docs/reference/convert_deg_rad.html +++ b/docs/reference/convert_deg_rad.html @@ -101,7 +101,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/export_surveyidx.html b/docs/reference/export_surveyidx.html index 5105fe6..5aa0908 100644 --- a/docs/reference/export_surveyidx.html +++ b/docs/reference/export_surveyidx.html @@ -112,7 +112,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/fix_age_group.html b/docs/reference/fix_age_group.html index d66467c..5beb311 100644 --- a/docs/reference/fix_age_group.html +++ b/docs/reference/fix_age_group.html @@ -109,7 +109,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_bathy_grid.html b/docs/reference/get_bathy_grid.html index 7fb9a90..e8affe5 100644 --- a/docs/reference/get_bathy_grid.html +++ b/docs/reference/get_bathy_grid.html @@ -125,7 +125,7 @@

Value

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_datrasraw.html b/docs/reference/get_datrasraw.html index 8a95796..a981ea6 100644 --- a/docs/reference/get_datrasraw.html +++ b/docs/reference/get_datrasraw.html @@ -125,7 +125,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_effect.html b/docs/reference/get_effect.html index f2ffede..c74d77d 100644 --- a/docs/reference/get_effect.html +++ b/docs/reference/get_effect.html @@ -108,7 +108,7 @@

Value

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_grid.html b/docs/reference/get_grid.html index 7bddf13..a422769 100644 --- a/docs/reference/get_grid.html +++ b/docs/reference/get_grid.html @@ -92,7 +92,7 @@

Value

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_prediction_grid.html b/docs/reference/get_prediction_grid.html index dc6939e..c57ef6d 100644 --- a/docs/reference/get_prediction_grid.html +++ b/docs/reference/get_prediction_grid.html @@ -133,7 +133,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_surveyidx.html b/docs/reference/get_surveyidx.html index 38e8e4c..722b765 100644 --- a/docs/reference/get_surveyidx.html +++ b/docs/reference/get_surveyidx.html @@ -262,7 +262,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_surveyidx_aic.html b/docs/reference/get_surveyidx_aic.html index efb6b58..00b5b74 100644 --- a/docs/reference/get_surveyidx_aic.html +++ b/docs/reference/get_surveyidx_aic.html @@ -96,7 +96,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_surveyidx_resid.html b/docs/reference/get_surveyidx_resid.html index 9845e85..c362ccb 100644 --- a/docs/reference/get_surveyidx_resid.html +++ b/docs/reference/get_surveyidx_resid.html @@ -96,7 +96,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_surveyidx_sim.html b/docs/reference/get_surveyidx_sim.html index 83a1334..323a8df 100644 --- a/docs/reference/get_surveyidx_sim.html +++ b/docs/reference/get_surveyidx_sim.html @@ -104,7 +104,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/get_surveyidx_stratmean.html b/docs/reference/get_surveyidx_stratmean.html index b46586d..eeae8aa 100644 --- a/docs/reference/get_surveyidx_stratmean.html +++ b/docs/reference/get_surveyidx_stratmean.html @@ -100,7 +100,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/index.html b/docs/reference/index.html index 5be9a54..97c15c8 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -139,6 +139,10 @@

All functions noaa_afsc_public_foss

Public data from FOSS

+ +

plot_simulation_list()

+ +

Plot survey index list (e.g. retrospective analysis)

plot_surveyidx()

@@ -172,7 +176,7 @@

All functions
-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/noaa_afsc_public_foss.html b/docs/reference/noaa_afsc_public_foss.html index b593271..a377900 100644 --- a/docs/reference/noaa_afsc_public_foss.html +++ b/docs/reference/noaa_afsc_public_foss.html @@ -196,7 +196,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/plot_simulation_list.html b/docs/reference/plot_simulation_list.html new file mode 100644 index 0000000..03baf64 --- /dev/null +++ b/docs/reference/plot_simulation_list.html @@ -0,0 +1,139 @@ + +Plot survey index list (e.g. retrospective analysis) — plot_simulation_list • surveyIndex + + +
+
+ + + +
+
+
+

Plot survey index list (e.g. retrospective analysis)

+ + +
+ +
+

Plot survey index list (e.g. retrospective analysis)

+
+ +
+ 
plot_simulation_list(
+  x,
+  base = 1,
+  rescale = FALSE,
+  lwd = 1.5,
+  main = NULL,
+  allCI = FALSE,
+  includeCI = TRUE,
+  ylim = NULL
+)
+
+ +
+

Arguments

+
x
+

(named) list of "surveyIdx" objects for example from "retro.surveyIdx" or "leaveout.surveyIdx"

+ + +
base
+

Either index of x that should considered the "base run" (integer), OR object of class "surveyIdx". Confidence bounds will be shown for this model only.

+ + +
rescale
+

Should indices be rescaled to have mean 1 (over the set of intersecting years)? Default: FALSE

+ + +
lwd
+

line width argument to plot

+ + +
main
+

if not NULL override main plotting default title of "Age group a"

+ + +
allCI
+

show 95% confidence lines for all indices? Default FALSE.

+ + +
includeCI
+

Show confidence intervals? Default TRUE.

+ + +
ylim
+

Y axis range. If NULL (default) then determine automatically.

+ +
+
+

Value

+ + +

nothing

+
+ +
+ +
+ + +
+

Developed by Emily Markowitz, Margaret Siple, Casper Berg.

+
+ +
+

Site built with pkgdown 2.0.6.9000.

+
+ +
+ + + + + + + + diff --git a/docs/reference/plot_surveyidx.html b/docs/reference/plot_surveyidx.html index 77a3dc6..c168259 100644 --- a/docs/reference/plot_surveyidx.html +++ b/docs/reference/plot_surveyidx.html @@ -191,7 +191,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/plot_surveyidx_grid.html b/docs/reference/plot_surveyidx_grid.html index d1f458c..7481593 100644 --- a/docs/reference/plot_surveyidx_grid.html +++ b/docs/reference/plot_surveyidx_grid.html @@ -102,7 +102,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/pred_grid_ebs.html b/docs/reference/pred_grid_ebs.html index 5632d56..8f8a26d 100644 --- a/docs/reference/pred_grid_ebs.html +++ b/docs/reference/pred_grid_ebs.html @@ -106,7 +106,7 @@

Examples

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/qres_tweedie.html b/docs/reference/qres_tweedie.html index ec57f33..44b9f53 100644 --- a/docs/reference/qres_tweedie.html +++ b/docs/reference/qres_tweedie.html @@ -92,7 +92,7 @@

Details

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/docs/reference/redo_surveyidx.html b/docs/reference/redo_surveyidx.html index 0c9f659..04f263d 100644 --- a/docs/reference/redo_surveyidx.html +++ b/docs/reference/redo_surveyidx.html @@ -120,7 +120,7 @@

Value

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.6.9000.

diff --git a/man/plot_simulation_list.Rd b/man/plot_simulation_list.Rd new file mode 100644 index 0000000..32c3357 --- /dev/null +++ b/man/plot_simulation_list.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/functions.R +\name{plot_simulation_list} +\alias{plot_simulation_list} +\title{Plot survey index list (e.g. retrospective analysis)} +\usage{ +plot_simulation_list( + x, + base = 1, + rescale = FALSE, + lwd = 1.5, + main = NULL, + allCI = FALSE, + includeCI = TRUE, + ylim = NULL +) +} +\arguments{ +\item{x}{(named) list of "surveyIdx" objects for example from "retro.surveyIdx" or "leaveout.surveyIdx"} + +\item{base}{Either index of x that should considered the "base run" (integer), OR object of class "surveyIdx". Confidence bounds will be shown for this model only.} + +\item{rescale}{Should indices be rescaled to have mean 1 (over the set of intersecting years)? Default: FALSE} + +\item{lwd}{line width argument to plot} + +\item{main}{if not NULL override main plotting default title of "Age group a"} + +\item{allCI}{show 95\% confidence lines for all indices? Default FALSE.} + +\item{includeCI}{Show confidence intervals? Default TRUE.} + +\item{ylim}{Y axis range. If NULL (default) then determine automatically.} +} +\value{ +nothing +} +\description{ +Plot survey index list (e.g. retrospective analysis) +} diff --git a/notes.R b/notes.R index 845ed28..665fae3 100644 --- a/notes.R +++ b/notes.R @@ -35,13 +35,14 @@ library(here) library(usethis) library(pkgdown) -# devtools::install_github("r-lib/pkgdown") -# pkgdown::build_favicons() -# devtools::build_vignettes() -# usethis::use_pkgdown(config_file = "./docs/pkgdown.yml") +# Run once to configure package to use pkgdown +# usethis::use_pkgdown() +# Run to build the website +pkgdown::build_site() -pkgdown::build_site(pkg = here::here()) -# usethis::use_github_action("pkgdown") + +# If you’re using GitHub, we also recommend setting up GitHub actions to automatically build and publish your site: +# usethis::use_pkgdown_github_pages() # Save Package tar.gz version0 <- "0.1.0" diff --git a/vignettes/A-data-prep.Rmd b/vignettes/A-data-prep.Rmd index 04e5fa0..31614f9 100644 --- a/vignettes/A-data-prep.Rmd +++ b/vignettes/A-data-prep.Rmd @@ -27,6 +27,8 @@ PKG <- c( "sf", "gstat", "magrittr", + "raster", + # RACE-GAP Specific R packages "akgfmaps", # devtools::install_github("sean-rohan-noaa/akgfmaps", build_vignettes = TRUE) "coldpool" # devtools::install_github("sean-rohan-noaa/coldpool") @@ -158,8 +160,8 @@ dat_wrangled <- dat %>% Quarter = "2") %>% dplyr::mutate(across((c("Year", "Ship", "COMMON_NAME")), as.factor)) %>% dplyr::select(wCPUE, GEAR_TEMPERATURE, BOTTOM_DEPTH, COMMON_NAME, EFFORT, - Year, Ship, Lon, Lat, lat, lon, sx, sy, ctime, - TimeShotHour, timeOfYear, Gear, Quarter, HaulDur, hauljoin) + Year, Ship, Lon, Lat, lat, lon, sx, sy, + ctime, TimeShotHour, timeOfYear, Gear, Quarter, HaulDur, hauljoin) head(dat_wrangled) @@ -190,7 +192,7 @@ pred_grid <- pred_grid %>% head(pred_grid) ``` - + ```{r survey_shapefile, eval=FALSE, echo=FALSE} @@ -219,7 +221,7 @@ Here we want to match covariate data to the prediction grid. dat_cov <- surveyIndex::pred_grid_ebs %>% dplyr::select(-Shape_Area) %>% - dplyr::mutate( + dplyr::mutate( sx = ((lon - mean(lon, na.rm = TRUE))/1000), sy = ((lat - mean(lat, na.rm = TRUE))/1000)) @@ -238,7 +240,7 @@ For this, we are going to create a raster of depth in the Bering sea from the su ```{r covar_depth_load, echo=FALSE} load(file = system.file(paste0("VigA_bottom_depth_raster_", - min(YEARS),"-",max(YEARS), ".rdata"), + min(YEARS),"-",max(YEARS), ".rdata"), package = "surveyIndex" ) ) ``` @@ -269,7 +271,7 @@ extrap_data <- stars::st_extract(x = extrap_data0, # to make future runs of this faster: save(extrap_data0, extrap_data, file = paste0("../inst/VigA_bottom_depth_raster_", - min(YEARS),"-",max(YEARS), ".rdata")) + min(YEARS),"-",max(YEARS), ".rdata")) ``` @@ -278,7 +280,7 @@ save(extrap_data0, extrap_data, plot(extrap_data0, main = "Interpolated Bottom Depths") dat_cov <- cbind.data.frame(dat_cov, - "BOTTOM_DEPTH" = extrap_data$var1.pred) %>% + "BOTTOM_DEPTH" = extrap_data$var1.pred) %>% stats::na.omit() head(dat_cov) @@ -301,12 +303,12 @@ for (i in 1:length(YEARS)) { } extrap_data0 <- coldpool::ebs_bottom_temperature[[tmp]] %>% - as(., Class = "SpatialPointsDataFrame") %>% - sf::st_as_sf() %>% + as(., Class = "SpatialPointsDataFrame") %>% + sf::st_as_sf() %>% sf::st_transform(crs = crs_latlon) %>% stars::st_rasterize() %>% - stars::st_extract(x = ., - at = as.matrix(dat_cov[,c("lon", "lat")])) + stars::st_extract(x = ., + at = as.matrix(dat_cov[,c("lon", "lat")])) names(extrap_data0) <- paste0("GEAR_TEMPERATURE", YEARS) @@ -325,62 +327,62 @@ Now, we need to fill in the data with the zeros! ```{r catch_haul_fill_0} - # Identify vars that will be used -------------------------------------------- - - varsbyyr <- unique( # c("GEAR_TEMPERATURE", "cpi") - gsub(pattern = "[0-9]+", - replacement = "", - x = names(dat_cov)[grepl(names(dat_cov), - pattern = "[0-9]+")])) - - vars <- unique( # c("BOTTOM_DEPTH") - names(dat_cov)[!grepl(names(dat_cov), - pattern = "[0-9]+")]) - vars <- vars[!(vars %in% c("LONG", "LAT", "lon", "lat", "sx", "sy"))] - - ## Fill catch data with zeros --------------------------------------------------------- - - data_hauls <- dat_wrangled %>% - dplyr::select(Year, sx, sy, - dplyr::all_of(varsbyyr), dplyr::all_of(vars), - Ship, hauljoin, - lat, lon, - ctime, TimeShotHour, timeOfYear, Gear, Quarter, - EFFORT, HaulDur) %>% - # dplyr::filter(!is.na(GEAR_TEMPERATURE)) %>% - na.omit() %>% - dplyr::distinct() - - data_catch <- dat_wrangled %>% - dplyr::select(COMMON_NAME, wCPUE, hauljoin) - - dat_catch_haul <- dplyr::left_join(x = data_hauls, - y = data_catch, - by = c("hauljoin")) %>% - dplyr::mutate(wCPUE = ifelse(is.na(wCPUE), 0, wCPUE)) - - head(dat_catch_haul) +# Identify vars that will be used -------------------------------------------- + +varsbyyr <- unique( # c("GEAR_TEMPERATURE", "cpi") + gsub(pattern = "[0-9]+", + replacement = "", + x = names(dat_cov)[grepl(names(dat_cov), + pattern = "[0-9]+")])) + +vars <- unique( # c("BOTTOM_DEPTH") + names(dat_cov)[!grepl(names(dat_cov), + pattern = "[0-9]+")]) +vars <- vars[!(vars %in% c("LONG", "LAT", "lon", "lat", "sx", "sy"))] + +## Fill catch data with zeros --------------------------------------------------------- + +data_hauls <- dat_wrangled %>% + dplyr::select(Year, sx, sy, + dplyr::all_of(varsbyyr), dplyr::all_of(vars), + Ship, hauljoin, + lat, lon, Lat, Lon, + ctime, TimeShotHour, timeOfYear, Gear, Quarter, + EFFORT, HaulDur) %>% + # dplyr::filter(!is.na(GEAR_TEMPERATURE)) %>% + na.omit() %>% + dplyr::distinct() + +data_catch <- dat_wrangled %>% + dplyr::select(COMMON_NAME, wCPUE, hauljoin) + +dat_catch_haul <- dplyr::left_join(x = data_hauls, + y = data_catch, + by = c("hauljoin")) %>% + dplyr::mutate(wCPUE = ifelse(is.na(wCPUE), 0, wCPUE)) + +head(dat_catch_haul) ``` - + ```{r catch_haul_datras } - allpd <- lapply(YEARS, FUN = surveyIndex::get_prediction_grid, x = dat_cov, - vars = vars, varsbyyr = varsbyyr) - names(allpd) <- as.character(YEARS) - - head(allpd[1][[1]]) - +allpd <- lapply(YEARS, FUN = surveyIndex::get_prediction_grid, x = dat_cov, + vars = vars, varsbyyr = varsbyyr) +names(allpd) <- as.character(YEARS) + +head(allpd[1][[1]]) + ``` ### 5a. Covariate Data ```{r cov_datras} - ## split data by species, make into DATRASraw + Nage matrix - ds <- split(dat_catch_haul,dat_catch_haul$COMMON_NAME) - ds <- lapply(ds, surveyIndex::get_datrasraw) - ## OBS, response is added here in "Nage" matrix -- use wCPUE - ds <- lapply(ds,function(x) { x[[2]]$Nage <- matrix(x$wCPUE,ncol=1); colnames(x[[2]]$Nage)<-1; x } ) - +## split data by species, make into DATRASraw + Nage matrix +ds <- split(dat_catch_haul,dat_catch_haul$COMMON_NAME) +ds <- lapply(ds, surveyIndex::get_datrasraw) +## OBS, response is added here in "Nage" matrix -- use wCPUE +ds <- lapply(ds,function(x) { x[[2]]$Nage <- matrix(x$wCPUE,ncol=1); colnames(x[[2]]$Nage)<-1; x } ) + ds ``` @@ -392,10 +394,10 @@ fm <- list( "fm_1_s_t_st" = "Year + s(sx,sy,bs=c('ts'),k=376) + s(sx,sy,bs=c('ts'),k=10,by=Year)", - + # Mdoel with simple covariates "fm_2_cov" = -"s(BOTTOM_DEPTH,bs='ts',k=10) + + "s(BOTTOM_DEPTH,bs='ts',k=10) + s(log(GEAR_TEMPERATURE+3),bs='ts',k=10)" ) ``` @@ -406,16 +408,16 @@ s(log(GEAR_TEMPERATURE+3),bs='ts',k=10)" Here are all of the models we want to try fitting: ```{r model_combos} - comb <- tidyr::crossing( - "SPECIES" = SPECIES, - "fm_name" = gsub(pattern = " ", replacement = "_", x = names(fm))) %>% - dplyr::left_join( - x = ., - y = data.frame("fm" = gsub(pattern = "\n", replacement = "", - x = unlist(fm), fixed = TRUE), - "fm_name" = gsub(pattern = " ", replacement = "_", - x = names(fm))), - by = "fm_name") +comb <- tidyr::crossing( + "SPECIES" = SPECIES, + "fm_name" = gsub(pattern = " ", replacement = "_", x = names(fm))) %>% + dplyr::left_join( + x = ., + y = data.frame("fm" = gsub(pattern = "\n", replacement = "", + x = unlist(fm), fixed = TRUE), + "fm_name" = gsub(pattern = " ", replacement = "_", + x = names(fm))), + by = "fm_name") comb ``` @@ -425,100 +427,243 @@ load(system.file("VigA_model_fits.rdata", package = "surveyIndex") ) ``` ```{r model_fit, eval = FALSE} - - models <- fittimes <- list() - - for(i in 1:nrow(comb)){ - cat("Fitting ",comb$SPECIES[i],"\n", comb$fm_name[i], ": ", comb$fm[i], "\n") - - temp <- paste0(comb$SPECIES[i], " ", comb$fm_name[i]) - - fittimes[[ temp ]] <- - system.time ( models[[ temp ]] <- - surveyIndex::get_surveyidx( - x = ds[[ comb$SPECIES[i] ]], - ages = 1, - myids = NULL, - predD = allpd, - cutOff = 0, - fam = "Tweedie", - modelP = comb$fm[i], - gamma = 1, - control = list(trace = TRUE, - maxit = 20)) ) - } +models <- fittimes <- list() - save(models, fittimes, file = paste0("../inst/VigA_model_fits.Rdata")) +for(i in 1:nrow(comb)){ + cat("Fitting ",comb$SPECIES[i],"\n", comb$fm_name[i], ": ", comb$fm[i], "\n") + + temp <- paste0(comb$SPECIES[i], " ", comb$fm_name[i]) + + fittimes[[ temp ]] <- + system.time ( models[[ temp ]] <- + surveyIndex::get_surveyidx( + x = ds[[ comb$SPECIES[i] ]], + ages = 1, + myids = NULL, + predD = allpd, + cutOff = 0, + fam = "Tweedie", + modelP = comb$fm[i], + gamma = 1, + control = list(trace = TRUE, + maxit = 20)) ) + +} + +save(models, fittimes, file = paste0("../inst/VigA_model_fits.Rdata")) ``` ```{r model_check1} - ## Check basis dimensions splines (spatial resolution) - # sink(paste0(".", dir_out, "gamcheck.txt")) +## Check basis dimensions splines (spatial resolution) +# sink(paste0(".", dir_out, "gamcheck.txt")) par(mfrow = c(2,2)) - lapply(models,function(x) gam.check(x$pModels[[1]])) - # sink() +lapply(models,function(x) gam.check(x$pModels[[1]])) +# sink() ``` ```{r model_check2} - ## Model summaries - # sink(paste0(".", dir_out, "summaries.txt")) - lapply(models,function(x) summary(x$pModels[[1]])) - # sink() +## Model summaries +# sink(paste0(".", dir_out, "summaries.txt")) +lapply(models,function(x) summary(x$pModels[[1]])) +# sink() + +surveyIndex::get_surveyidx_aic(x = models) + - surveyIndex::get_surveyidx_aic(x = models) - - ``` ``` {r mode_aic} - temp <- sapply(models, `[`, "pModels") - mods <- sapply(temp, `[`, 1) +temp <- sapply(models, `[`, "pModels") +mods <- sapply(temp, `[`, 1) - lapply(X = mods, FUN = AIC) +lapply(X = mods, FUN = AIC) ``` ## 8. Indicies of Abundance +```{r design_based_load, eval = FALSE, echo = FALSE} +design_based <- read.csv(file = "../output/cpue_ebs_plusnw.csv") %>% + dplyr::filter(#srvy == SRVY & + YEAR %in% YEARS & + COMMON_NAME %in% SPECIES) +``` + ```{r indicie_abund} + dat <- data.frame() for (i in 1:length(models)){ temp <- models[[i]] - dat <- dplyr::bind_rows(dat, - data.frame(idx = temp$idx[,1], - Year = rownames(temp$idx), - group = names(models)[i])) + dat0 <- data.frame(idx = temp$idx[,1], + lo = temp$lo[,1], + up = temp$up[,1], + Year = rownames(temp$idx), + group = names(models)[i], + formula = paste0("cpue_kgha ~ ", + as.character(temp$pModels[[1]]$formula)[[3]])) + + dat <- dplyr::bind_rows(dat, dat0) } -dat$facet_group <- sapply(X = strsplit(x = dat$group, split = " "), - `[`, 1) +dat$facet_group <- paste0(sapply(X = strsplit(x = dat$group, split = " fm"), `[`, 1)) +# dat$model <- paste0(sapply(X = strsplit(x = dat$group, split = " fm"), `[`, 2)) -dat$idx[dat$Year == 2020] <- NA +dat[dat$Year == 2020, c("idx", "up", "lo")] <- NA ggplot2::ggplot(data = dat, - mapping = aes(x = Year, - y = idx, - group = group, - color = group)) + + mapping = aes(x = Year, + y = idx, + group = formula, + color = formula)) + geom_line(size = 1.5) + - geom_point(size = 2) + + geom_point(size = 2) + + geom_ribbon(aes(ymin = lo, ymax = up, fill = formula), + alpha=0.1, + linetype="dashed", + color="grey") + + ggtitle("Annual Index Model Results") + facet_wrap(vars(facet_group), scales = "free", ncol = 1) + - theme(legend.position = "bottom") + theme(legend.position = "bottom", + legend.direction = "vertical") ``` - -## 8. Predict and plot +## 9. Predict and plot ```{r predict} -# predict.gam(object = models$`red king crab fm_2_cov`, -# newdata = dat_catch_haul %>% -# dplyr::filter(Year == 2021) %>% -# dplyr::select(Year, sx, sy, GEAR_TEMPERATURE, BOTTOM_DEPTH)) +dat_pred <- dat_catch_haul %>% + dplyr::select(Year, sx, sy, Lon, Lat, GEAR_TEMPERATURE, BOTTOM_DEPTH) + +dat <- data.frame() +for (i in 1:length(models)) { + temp <- models[[i]] + dat0 <- data.frame(idx = + predict.gam( + object = temp$pModels[[1]], + newdata = dat_pred), + group = names(models)[i], + formula = paste0("cpue_kgha ~ ", + as.character(temp$pModels[[1]]$formula)[[3]]) + ) + dat00 <- dplyr::bind_cols(dat0, dat_pred) + dat <- dplyr::bind_rows(dat, dat00) + +# dat_r <- raster::rasterFromXYZ(xyz = dat00[,c("lon", "lat", "idx")]) + +} + +dat$facet_group <- paste0(sapply(X = strsplit(x = dat$group, split = " fm"), `[`, 1)) + +for (i in 1:length(unique(dat$facet_group))){ + + ggplot2::ggplot(data = dat %>% + dplyr::filter(facet_group == unique(dat$facet_group)[i]), + mapping = aes(x = Lon, + y = Lat, + group = group, + color = idx)) + + geom_point() + + ggtitle(paste0("Annual Index Model Results for ", unique(dat$facet_group)[i])) + + facet_grid(cols = vars(group), + rows = vars(Year)) + + theme_bw() + +} +``` + +```{r} + +surveyIndex::plot_surveyidx( + x = models, + dat = ds, + myids = NULL, + predD = allpd) ``` +## 10. Simulations + +```{r sim_gam, echo=FALSE, eval=FALSE} + +# sims <- fittimes_sims <- list() +# for(i in 1:nrow(comb)){ +# +# cat("Simulating ",comb$SPECIES[i],"\n", comb$fm_name[i], ": ", comb$fm[i], "\n") +# +# temp <- paste0(comb$SPECIES[i], " ", comb$fm_name[i]) +# +# fittimes[[ temp ]] <- +# system.time ( sims[[ temp ]] <- +# surveyIndex::get_surveyidx_sim( +# model = models[[i]], +# d = ds[[ comb$SPECIES[i] ]]) ) +# } +# +# par(mfrow = c(2, 2)) # Create a 2 x 2 plotting matrix +# for(i in 1:nrow(comb)){ +# plot(sims[[i]]$sim, main = paste0(names(sims)[i], " sims")) +# plot(sims[[i]]$mu[[1]], main = paste0(names(sims)[i], " mu")) +# } +``` + +```{r sim_gam1, eval=FALSE} +REPS <- 4 +ests <- list() + +for(i in 1:nrow(comb)){ + + cat("Simulating ",comb$SPECIES[i],"\n", comb$fm_name[i], ": ", comb$fm[i], "\n") + temp <- paste0(comb$SPECIES[i], " ", comb$fm_name[i]) + +# for(SPECIES in specLevels){ + ests[[ temp ]] <- list() + + ## simulate data + csim <- surveyIndex::get_surveyidx_sim(models[[i]], ds[[comb$SPECIES[i]]]) + sims <-lapply(1:REPS,function(j) surveyIndex::get_surveyidx_sim( + model = models[[i]], + d = ds[[comb$SPECIES[i]]], + sampleFit = FALSE, + condSim = csim) ) + + ## re-estimate + tmp <- ds[[i]] + for(i in 1:REPS) { + tmp[[2]]$Nage <- matrix(sims[[i]][[1]][,1],ncol=1) + colnames(tmp$Nage)<-1 + + ests[[SPECIES]][[i]] <- + surveyIndex::get_surveyidx( + x = tmp, + ages = 1, + myids=NULL, + predD=allpd, + cutOff=0, + fam="Tweedie", + modelP=fm, + gamma=1, + control=list(trace=TRUE,maxit=10)) + # cat(i, " done.\n") + } + +} + +png("simest.png",width=640*pngscal,height=480*pngscal) +par(mfrow=c(2,2)) + +for(i in 1:nrow(comb)){ + +# for(SPECIES in specLevels){ + surveyIndex::plot_simulation_list( + x = ests[[temp]], + base=models[[temp]], + main=temp, + lwd=2) +} +dev.off() + +```