Benthos Indices: Results

Introduction

Models for macrobenthos and megabenthos were run for spring and fall using all spatial and spatio-temporal random effects, and three catchability covariates: mean predator length at the location, number of predator species at the location, and bottom temperature at the location.

Methods

See workflow document for details.

Results

Counts of stations total and with prey

In the table below, Total stations is all stations where benthivore predators were sampled. The fields from Total stomachs through Variance n preds reflect the number of predators sampled at all of these stations. The last four columns reflect the number and percent of stations with positive macrobenthos and megabenthos prey, respectively.

Note that we did not use data from 1973-1979 in the models. As of October 2024, we are including data through 2023 from both NEFSC and NEAMAP.

macrobenagg_stn <- readRDS(here::here("fhdata/macrobenagg_stn_all_modBT.rds"))
megabenagg_stn <- readRDS(here::here("fhdata/megabenagg_stn_all_modBT.rds"))

nstationsyr <- macrobenagg_stn |>
  dplyr::group_by(year, season_ng) |>
  dplyr::summarise(totstations = n(),
                   totstomachs = sum(nstomtot),
                   minnpred = min(nbenthivoresp),
                   meannpred = mean(nbenthivoresp),
                   maxnpred = max(nbenthivoresp),
                   varnpred = var(nbenthivoresp))

macropreystns <- macrobenagg_stn |>
  dplyr::filter(meanmacrobenpywt > 0) |>
  dplyr::group_by(year, season_ng) |>
  dplyr::summarise(macrobenstations = n())

megapreystns <- megabenagg_stn |>
  dplyr::filter(meanmegabenpywt > 0) |>
  dplyr::group_by(year, season_ng) |>
  dplyr::summarise(megabenstations = n())

stationsprey <- dplyr::left_join(nstationsyr, macropreystns) |>
  dplyr::left_join(megapreystns) |>
  dplyr::arrange(desc(season_ng), .by_group = TRUE) |>
  dplyr::mutate(year = as.character(year),
                percmacro = round(macrobenstations/totstations*100),
                percmega = round(megabenstations/totstations*100)) |>
  na.omit()

flextable::flextable(stationsprey) |>
  flextable::colformat_char(j = "year") |>
  flextable::fontsize(size = 8, part = "all") |>
  flextable::set_header_labels(
    year = "Year", season_ng = "Season", totstations = "Total stations", 
    totstomachs = "Total stomachs", minnpred = "Minimum n preds", meannpred = "Mean n preds", 
    maxnpred = "Maximum n preds", varnpred = "Variance n preds",
    macrobenstations =  "Macrobenthos stations",
    megabenstations = "Megabenthos stations",
    percmacro =  "Percent macrobenthos", percmega = "Percent megabenthos")

Year	Season	Total stations	Total stomachs	Minimum n preds	Mean n preds	Maximum n preds	Variance n preds	Macrobenthos stations	Megabenthos stations	Percent macrobenthos	Percent megabenthos
1973	SPRING	148	1,437	1	1.655405	5	0.7443924	110	35	74	24
1973	FALL	183	2,493	1	2.142077	6	1.4522308	151	61	83	33
1974	SPRING	165	1,972	1	2.018182	5	1.3106430	126	70	76	42
1974	FALL	186	1,886	1	1.838710	6	1.1630340	149	57	80	31
1975	SPRING	34	328	1	1.500000	5	0.8030303	27	17	79	50
1975	FALL	162	1,812	1	1.864198	5	1.0000767	135	51	83	31
1976	SPRING	182	2,183	1	2.076923	5	1.6404590	152	49	84	27
1976	FALL	163	2,105	1	2.024540	6	1.3203817	132	66	81	40
1977	SPRING	217	2,591	1	2.552995	7	2.5446322	165	68	76	31
1977	FALL	232	2,395	1	2.043103	8	2.0500821	162	93	70	40
1978	SPRING	271	2,887	1	2.081181	11	2.1785705	187	75	69	28
1978	FALL	495	5,372	1	2.115152	8	2.1871181	334	185	67	37
1979	SPRING	294	2,377	1	1.880952	7	1.1837315	199	79	68	27
1979	FALL	545	5,999	1	2.137615	10	2.1666892	351	214	64	39
1980	SPRING	269	2,466	1	1.836431	6	1.1373245	171	75	64	28
1980	FALL	367	3,844	1	2.013624	10	1.6254969	233	113	63	31
1981	SPRING	262	3,246	1	1.793893	5	0.8539089	144	113	55	43
1981	FALL	230	2,801	1	1.700000	6	0.9532751	96	72	42	31
1982	SPRING	372	4,988	1	2.190860	10	2.3812639	234	104	63	28
1982	FALL	231	2,210	1	1.865801	6	1.1775645	87	57	38	25
1983	SPRING	297	2,455	1	1.666667	7	0.9662162	110	65	37	22
1983	FALL	206	1,736	1	1.728155	5	0.9111058	64	41	31	20
1984	SPRING	309	2,312	1	1.394822	5	0.4929601	135	74	44	24
1984	FALL	255	2,604	1	1.956863	7	1.5375019	106	69	42	27
1985	SPRING	314	6,213	1	3.098726	9	3.7506054	201	118	64	38
1985	FALL	250	4,377	1	2.980000	9	3.3690763	104	80	42	32
1986	SPRING	311	5,871	1	3.395498	9	3.2979152	191	121	61	39
1986	FALL	249	4,112	1	3.116466	12	4.2807358	132	119	53	48
1987	SPRING	304	5,711	1	3.483553	10	3.3264613	234	145	77	48
1987	FALL	269	3,909	1	2.762082	8	2.5551240	148	109	55	41
1988	SPRING	250	4,036	1	3.104000	8	3.0654458	156	102	62	41
1988	FALL	262	4,152	1	3.122137	12	3.7321517	142	100	54	38
1989	SPRING	286	7,273	1	4.625874	11	3.6735738	247	175	86	61
1989	FALL	296	7,547	1	4.405405	12	4.7028859	222	151	75	51
1990	SPRING	292	5,727	1	3.376712	9	2.8129266	221	159	76	54
1990	FALL	310	7,045	1	4.190323	12	4.6400355	204	151	66	49
1991	SPRING	305	6,743	1	4.144262	9	2.9856989	250	189	82	62
1991	FALL	386	7,662	1	4.145078	10	4.2853913	271	165	70	43
1992	SPRING	419	9,469	1	4.527446	11	4.5752018	325	204	78	49
1992	FALL	429	8,653	1	4.391608	13	4.3603032	317	214	74	50
1993	SPRING	438	10,931	1	5.271689	11	4.5553065	367	286	84	65
1993	FALL	424	10,001	1	4.938679	13	5.7409117	317	220	75	52
1994	SPRING	426	10,259	1	4.730047	11	4.4846009	323	260	76	61
1994	FALL	369	8,425	1	4.384824	11	5.2265082	275	219	75	59
1995	SPRING	460	12,275	1	5.313043	12	4.3898077	374	300	81	65
1995	FALL	410	11,162	1	4.739024	13	6.2324587	298	213	73	52
1996	SPRING	460	12,497	1	5.191304	12	5.2312968	381	290	83	63
1996	FALL	338	7,348	1	4.639053	12	5.7981142	229	191	68	57
1997	SPRING	437	9,848	1	4.574371	11	4.4743875	324	225	74	51
1997	FALL	312	5,692	1	4.208333	11	5.4387728	188	178	60	57
1998	SPRING	497	13,900	1	5.895372	13	5.6422568	409	308	82	62
1998	FALL	365	8,929	1	4.835616	12	4.8245522	272	233	75	64
1999	SPRING	506	13,686	1	7.638340	19	8.4768638	425	330	84	65
1999	FALL	353	7,083	1	6.861190	15	8.0062452	288	219	82	62
2000	SPRING	483	11,630	1	8.020704	17	8.9497779	415	321	86	66
2000	FALL	330	5,901	1	6.645455	17	12.6915446	238	187	72	57
2001	SPRING	479	10,616	1	7.893528	18	10.8861296	399	277	83	58
2001	FALL	317	6,162	1	7.400631	19	12.7535439	246	169	78	53
2002	SPRING	472	11,730	1	8.879237	18	10.5565107	396	267	84	57
2002	FALL	317	6,183	1	7.189274	17	12.6792517	243	171	77	54
2003	SPRING	416	9,231	1	7.932692	17	11.5569045	336	230	81	55
2003	FALL	312	6,797	1	8.163462	16	7.5005256	261	190	84	61
2004	SPRING	458	9,605	1	7.788210	18	9.8346966	367	234	80	51
2004	FALL	299	5,225	1	6.685619	15	11.3504972	246	148	82	49
2005	SPRING	422	7,986	1	7.412322	17	13.1550022	340	209	81	50
2005	FALL	320	6,191	1	7.537500	17	12.0048589	278	176	87	55
2006	SPRING	460	9,453	1	8.919565	21	15.5556077	392	253	85	55
2006	FALL	354	7,030	1	8.375706	20	13.5666603	297	197	84	56
2007	SPRING	474	9,895	1	8.702532	18	13.5963149	407	295	86	62
2007	FALL	485	8,807	1	6.806186	16	10.3673170	448	331	92	68
2008	SPRING	480	9,186	1	7.031250	17	9.3246999	439	335	91	70
2008	FALL	487	9,065	1	7.205339	20	12.5791484	432	324	89	67
2009	SPRING	541	14,315	1	9.090573	21	18.2195591	510	384	94	71
2009	FALL	496	11,680	1	8.266129	21	16.3694363	441	340	89	69
2010	SPRING	521	13,137	1	9.180422	24	22.0058467	475	378	91	73
2010	FALL	483	9,918	1	7.511387	20	16.2628369	421	330	87	68
2011	SPRING	485	10,773	1	8.049485	18	16.2413479	437	325	90	67
2011	FALL	468	10,517	1	8.198718	20	20.0996129	432	349	92	75
2012	SPRING	499	13,536	1	10.118236	22	22.8233415	474	350	95	70
2012	FALL	498	11,708	1	8.598394	27	23.5969391	454	362	91	73
2013	SPRING	525	13,754	1	10.118095	19	13.4325918	510	368	97	70
2013	FALL	514	10,647	1	8.272374	25	18.3857070	475	339	92	66
2014	SPRING	438	10,624	1	9.230594	20	16.9192450	426	303	97	69
2014	FALL	497	10,138	1	8.056338	24	19.3556906	431	351	87	71
2015	SPRING	516	11,298	1	8.713178	21	16.8612629	482	332	93	64
2015	FALL	512	10,682	1	8.083984	22	21.8109367	445	354	87	69
2016	SPRING	485	12,366	1	9.668041	22	19.4288319	451	352	93	73
2016	FALL	513	10,480	1	8.083821	23	20.6394447	438	341	85	66
2017	SPRING	397	10,066	1	10.168766	20	16.0042745	369	279	93	70
2017	FALL	273	4,406	1	5.230769	10	3.9061086	204	144	75	53
2018	SPRING	397	8,467	1	8.256927	19	14.9186703	377	299	95	75
2018	FALL	444	8,062	1	7.443694	22	19.6559647	383	334	86	75
2019	SPRING	482	9,705	1	8.095436	17	14.2403533	442	347	92	72
2019	FALL	488	9,653	1	7.991803	22	19.3264214	401	339	82	69
2020	SPRING	122	1,749	1	7.483607	17	11.3096464	97	82	80	67
2020	FALL	142	2,632	1	5.260563	10	3.9103486	126	129	89	91
2021	SPRING	457	9,757	1	8.638950	21	16.8320761	429	344	94	75
2021	FALL	493	8,683	1	7.251521	24	16.0341694	411	348	83	71
2022	SPRING	495	9,954	1	8.224242	20	16.6277512	451	358	91	72
2022	FALL	449	8,693	1	7.772829	23	19.1848851	380	343	85	76
2023	SPRING	216	3,637	1	6.379630	20	16.8970715	188	145	87	67
2023	FALL	470	8,636	1	7.748936	20	20.5722270	347	282	74	60

Spatial distribution of stations

# from each output folder in pyindex, 
outdir <- here::here("pyindex")
moddirs <- list.dirs(outdir)
moddirs <- moddirs[-1]

# only latest
moddirs23 <- moddirs[-grep("1980-2022", moddirs)]

Spatial coverage for stomachs in fall, both groups (same set of predators, same stations)

mapfall <- paste0(moddirs23[1], "/Data_by_year.png")

knitr::include_graphics(mapfall) 

Spatial coverage for stomachs in spring, both groups (same set of predators, same stations)

mapspring <- paste0(moddirs23[2], "/Data_by_year.png")

knitr::include_graphics(mapspring) 

Model selection

See workflow document for details.

Model results

stratlook <- data.frame(Stratum = c("Stratum_1",
                                    "Stratum_2",
                                    "Stratum_3",
                                    "Stratum_4",
                                    "Stratum_5"),
                        Region  = c("AllEPU", 
                                    "MAB",
                                    "GB",
                                    "GOM",
                                    "SS"))

# function to apply extracting info
getmodindex <- function(d.name){
  # read settings
  modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
  modname <- modpath[length(modpath)]
  
  index <- read.csv(file.path(d.name, "Index.csv"))
  # return model indices as a dataframe
  out <- data.frame(modname = modname,
                    index
  )
  
  return(out)
}

modcompareindex <- purrr::map_dfr(moddirs23, getmodindex)

splitoutput <- modcompareindex %>%
  dplyr::mutate(Season = modname |> map(str_split, pattern = "_") |> map_chr(c(1,2))) %>%
  dplyr::mutate(Data = modname |> map(str_split, pattern = "_") |> map_chr(c(1,1))) %>%
  dplyr::left_join(stratlook) %>%
  dplyr::filter(Region %in% c("GOM", "GB", "MAB","SS", "AllEPU")) 

benthosmax <- max(splitoutput$Estimate)


benthostsmean <- splitoutput %>%
  dplyr::group_by(modname, Region) %>%
  dplyr::mutate(fmean = mean(Estimate)) 

Do we have ecodata trends in any time series? Some. Driven by early 80s spike?

filterEPUs <- c("MAB", "GB", "GOM", "SS", "AllEPU")

currentMonth <- lubridate::month(Sys.Date())
    currentYear <- lubridate::year(Sys.Date())
    if (currentMonth > 4) {
      endShade <- currentYear
    } else {
      endShade <- currentYear - 1
    }
    shadedRegion <- c(endShade-9,endShade)

    shade.alpha <- 0.3
    shade.fill <- "lightgrey"
    lwd <- 1
    pcex <- 2
    trend.alpha <- 0.5
    trend.size <- 2
    hline.size <- 1
    line.size <- 2
    hline.alpha <- 0.35
    hline.lty <- "dashed"
    label.size <- 5
    hjust.label <- 1.5
    letter_size <- 4
    errorbar.width <- 0.25
    x.shade.min <- shadedRegion[1]
    x.shade.max <- shadedRegion[2]

    setup <- list(
      shade.alpha = shade.alpha,
      shade.fill =shade.fill,
      lwd = lwd,
      pcex = pcex,
      trend.alpha = trend.alpha,
      trend.size = trend.size,
      line.size = line.size,
      hline.size = hline.size,
      hline.alpha = hline.alpha,
      hline.lty = hline.lty,
      errorbar.width = errorbar.width,
      label.size = label.size,
      hjust.label = hjust.label,
      letter_size = letter_size,
      x.shade.min = x.shade.min,
      x.shade.max = x.shade.max
    )
    

    fix<- splitoutput |>
      dplyr::filter(Data %in% c("macrobenthos"),
                    Region %in% filterEPUs) |>
      dplyr::group_by(Region, Season) |>
      dplyr::summarise(max = max(Estimate))
    
    p <- splitoutput |>
      dplyr::filter(Data %in% c("macrobenthos"),
                    Region %in% filterEPUs) |>
      dplyr::group_by(Region, Season) |>
      dplyr::left_join(fix) |>
      dplyr::mutate(#Value = Value/resca,
        Mean = as.numeric(Estimate),
        SE = Std..Error.for.Estimate,
        Mean = Mean/max,
        SE = SE/max,
        Upper = Mean + SE,
        Lower = Mean - SE) |>
      ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, group = Season))+
      ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
                        xmin = setup$x.shade.min , xmax = setup$x.shade.max,
                        ymin = -Inf, ymax = Inf) +
      ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Season), alpha = 0.5)+
      ggplot2::geom_point()+
      ggplot2::geom_line()+
      ggplot2::ggtitle("Macrobenthos")+
      ggplot2::ylab(expression("Relative benthos biomass"))+
      ggplot2::xlab(ggplot2::element_blank())+
      ggplot2::facet_wrap(Region~Season, ncol = 2)+
      ecodata::geom_gls()+
      ecodata::theme_ts()+
      ecodata::theme_facet()+
      ecodata::theme_title()
    
    p

    fix<- splitoutput |>
      dplyr::filter(Data %in% c("megabenthos"),
                    Region %in% filterEPUs) |>
      dplyr::group_by(Region, Season) |>
      dplyr::summarise(max = max(Estimate))
    
    p <- splitoutput |>
      dplyr::filter(Data %in% c("megabenthos"),
                    Region %in% filterEPUs) |>
      dplyr::group_by(Region, Season) |>
      dplyr::left_join(fix) |>
      dplyr::mutate(#Value = Value/resca,
        Mean = as.numeric(Estimate),
        SE = Std..Error.for.Estimate,
        Mean = Mean/max,
        SE = SE/max,
        Upper = Mean + SE,
        Lower = Mean - SE) |>
      ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, group = Season))+
      ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
                        xmin = setup$x.shade.min , xmax = setup$x.shade.max,
                        ymin = -Inf, ymax = Inf) +
      ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Season), alpha = 0.5)+
      ggplot2::geom_point()+
      ggplot2::geom_line()+
      ggplot2::ggtitle("Megabenthos")+
      ggplot2::ylab(expression("Relative benthos biomass"))+
      ggplot2::xlab(ggplot2::element_blank())+
      ggplot2::facet_wrap(Region~Season, ncol = 2)+
      ecodata::geom_gls()+
      ecodata::theme_ts()+
      ecodata::theme_facet()+
      ecodata::theme_title()
    
    p

Trend comparisons by season, any difference macro and megabenthos?

ggplot(benthostsmean |> filter(Season =="fall"), 
       aes(x=Time, y=((Estimate-fmean)/fmean), colour=Data)) +
  geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
                    ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative benthos biomass scaled to time series mean")  +
  ggtitle("Fall models: trend comparison") +
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

ggplot(benthostsmean |> filter(Season =="spring"), 
       aes(x=Time, y=((Estimate-fmean)/fmean), colour=Data)) +
  geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean, 
                    ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  #scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
  ylab("Relative benthos biomass scaled to time series mean")  +
  ggtitle("Spring models: trend comparison") +
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

Scale comparisons by season, any difference macro and megabenthos? Yes, way more macrobenthos, as expected.

ggplot(splitoutput |> filter(Season =="fall"), 
       aes(x=Time, y=Estimate, colour=Data)) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=Data), linetype = 0, alpha = 0.15)+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/benthosmax, digits = 2))+
  ylab("Relative benthos biomass scaled to maximum")  +
  ggtitle("Fall models: scale comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

ggplot(splitoutput |> filter(Season =="spring"), 
       aes(x=Time, y=Estimate, colour=Data)) +
  #geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=Data), linetype = 0, alpha = 0.15)+
  geom_point(aes(shape=Data))+
  geom_line()+
  #acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
  facet_wrap(~Region, scales = "free_y", ncol = 1) +
  scale_y_continuous(labels=function(x)round(x/benthosmax, digits = 2))+
  ylab("Relative benthos biomass scaled to maximum")  +
  ggtitle("Spring models: scale comparison")+
  theme(#legend.position = c(1, 0),
        #legend.justification = c(1, 0)
        legend.position = "bottom",
        legend.text = element_text(size=rel(0.5)),
        legend.key.size = unit(0.5, 'lines'),
        legend.title = element_text(size=rel(0.5)))

Center of gravity exploration

Has benthos shifted along the coast?

library(FishStatsUtils)

getcogVAST <- function(d.name){
  
  fit <- VAST::reload_model(readRDS(paste0(d.name,"/fit.rds")))
  
  dir.create(paste0(d.name,"/test"))
  
  cogout <- FishStatsUtils::plot_range_index(Sdreport = fit$parameter_estimates$SD, 
                                             Report = fit$Report, 
                                             TmbData = fit$data_list,
                                             year_labels = as.numeric(fit$year_labels),
                                             years_to_plot = fit$years_to_plot,
                                             Znames = colnames(fit$data_list$Z_xm),
                                             PlotDir = paste0(d.name,"/test")) #already have plots, will delete this directory
  
  saveRDS(cogout, paste0(d.name,"/cogout.rds"))
  
  unlink(paste0(d.name,"/test"), recursive=TRUE) #removes directory with unneeded plots
  
}

purrr::map(moddirs23, getcogVAST)

Center of gravity plots with ecodata trends

plotcog <- function(d.name){
  
  cogout <- readRDS(paste0(d.name,"/cogout.rds"))
  modpath <- unlist(str_split(d.name, pattern = "/"))
  modname <- modpath[length(modpath)]
  
  cogdat <- as.data.frame(cogout$COG_Table) |>
    dplyr::mutate(direction = ifelse(m==1, "Eastward", "Northward"))
  
  ggplot2::ggplot(cogdat, ggplot2::aes(x = Year, y = COG_hat)) + 
    ggplot2::geom_point() + 
    ecodata::geom_gls() + 
    ggplot2::labs(y = "Center of gravity, km") +
    ggplot2::facet_wrap(~direction, scales = "free_y") + 
    ecodata::theme_facet()+
    ggplot2::ggtitle(modname)
}

purrr::map(moddirs23, plotcog)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Density maps

for(d.name in moddirs23){
  
  modpath <- unlist(str_split(d.name, pattern = "/"))
  modname <- modpath[length(modpath)]
  
  cat(modname, "\n")
  cat(paste0("![](",d.name, "/ln_density-predicted.png)"))
  cat("\n\n")
  
}

macrobenthos_fall_500_cov

macrobenthos_spring_500_cov

megabenthos_fall_500_cov

megabenthos_spring_500_cov

What prey are in the raw input data?

Macrobenthos composition over time in stomachs

NEFSC data only

macrobenall_stn <- readRDS(here("fhdata/macrobenall_stn.rds"))

 MAB <- data.frame(stratum = c(1010:1080, 1100:1120, 1600:1750, 3010:3450, 3470, 3500, 3510),
                      EPU = "MAB")
 
 GB <- data.frame(stratum = c(1090, 1130:1210, 1230, 1250, 3460, 3480, 3490, 3520:3550),
                  EPU = "GB")
 
 GOM <- data.frame(stratum = c(1220, 1240, 1260:1290, 1360:1400, 3560:3830),
                   EPU = "GOM")
  
EPUlook <- dplyr::bind_rows(MAB, GB, GOM)

# object is called prey
load(here("data/prey.RData")) #original mappings from June 2023

addtomacro <- prey |>
  dplyr::filter(RPATH == "Megabenthos") |> 
  dplyr::filter(PYNAM %in% c("ALBUNEA PARETII", "EMERITA TALPOIDA", "RANILIA MURICATA", "HIPPIDAE"))

megaben <- prey |>
  dplyr::filter(RPATH == "Megabenthos") |>
  dplyr::filter(!PYNAM %in% c("ALBUNEA PARETII", "EMERITA TALPOIDA", "RANILIA MURICATA", "HIPPIDAE"))

macroben <- prey |>
  dplyr::filter(RPATH == "Macrobenthos") |> 
  dplyr::filter(!PYNAM %in% c("PECTINIDAE", "PECTINIDAE SHELL", "PECTINIDAE VISCERA")) |>
  dplyr::bind_rows(addtomacro)

macroben.analcoll <- macroben |>
  dplyr::select(PYNAM, PYCOMNAM, AnalSci, AnalCom, CollSci, Collcom)

collcolors <- data.frame(prey = unique(macroben.analcoll$Collcom),
                        #preycol = colors()[seq(1, length(colors()), 15)])
                        preycol = viridis::turbo(length(unique(macroben.analcoll$Collcom))))


macrobensppwt <- macrobenall_stn %>%
  dplyr::left_join(EPUlook) %>%
  dplyr::filter(macrobenthos=="macroben",
                !is.na(EPU),
                !is.na(season_ng)) %>%
  dplyr::group_by(year, season_ng, EPU, macrobenpynam) %>%
  dplyr::summarise(pysum = sum(macrobenpywt)) %>%
  dplyr::mutate(pyprop = pysum / sum(pysum)) %>%
  dplyr::ungroup() %>%
  dplyr::left_join(macroben.analcoll, by =  dplyr::join_by(macrobenpynam == PYNAM)) |>
  dplyr::left_join(collcolors, by =  dplyr::join_by(Collcom == prey))

preycols <- macrobensppwt$preycol
names(preycols) <- as.factor(macrobensppwt$macrobenpynam)

preyaggcols <- macrobensppwt$preycol
names(preyaggcols) <- as.factor(macrobensppwt$Collcom)


top20spp <- macrobensppwt |>
  dplyr::group_by(season_ng, EPU, Collcom) |>
  dplyr::summarise(meanprop = mean(pyprop, na.rm = TRUE),
                   maxprop = max(pyprop, na.rm = TRUE)) |>
  dplyr::ungroup() |>
  dplyr::arrange(EPU, season_ng, desc(meanprop))

plotDietCompBar <- function(dat, metric, catname, plotcol, title=NULL){   #defines the name of the function and the requied inputs
  
dat %>% 
  filter(!is.na(.data[[metric]])) %>%          #take out the NA values (literally "not is NA")
  ggplot(aes(x=year, y=.data[[metric]], fill=.data[[catname]])) +     #year on x axis, %diet comp on y
  geom_bar(width = 1, stat = "identity", position="fill") +    #stacked bar chart
  scale_fill_manual(values=plotcol) +         #custom colors
  ecodata::theme_facet() +                    #simplify
  facet_grid(fct_relevel(EPU,  "GOM", "GB", "MAB")~
               fct_relevel(season_ng, "SPRING", "FALL")) + #separate by ordered epu and season
  labs(y = "% composition",              #add sensible labels
       title = title) 
    
}   

addSmallLegend <- function(myPlot, pointSize = 2, textSize = 6, spaceLegend = 0.5) {
    myPlot +
        guides(shape = guide_legend(override.aes = list(size = pointSize)),
               color = guide_legend(override.aes = list(size = pointSize)),
               fill = guide_legend(ncol = 1)) +
        theme(legend.title = element_text(size = textSize), 
              legend.text  = element_text(size = textSize),
              legend.key.size = unit(spaceLegend, "lines"))
}

pall <- plotDietCompBar(dat=macrobensppwt, 
                          metric="pysum", 
                        catname="macrobenpynam",
                        plotcol=preycols,
                          title="Diet composition")  

pcat <- plotDietCompBar(dat=macrobensppwt, 
                          metric="pysum", 
                        catname="Collcom",
                        plotcol=preyaggcols,
                          title="Diet composition")

addSmallLegend(pall)

addSmallLegend(pcat)

Collcom trends by EPU

macrobensppwt |>
  dplyr::filter(EPU=="GOM") |>
  ggplot2::ggplot(ggplot2::aes(x=year, y=pyprop, colour = season_ng)) +
  ggplot2::geom_line() +
  ecodata::theme_facet() +
  ggplot2::facet_wrap(~Collcom) +
  ggplot2::ggtitle("GOM Macro") +
  ggplot2::theme(legend.position = "bottom")

macrobensppwt |>
  dplyr::filter(EPU=="GB") |>
  ggplot2::ggplot(ggplot2::aes(x=year, y=pyprop, colour = season_ng)) +
  ggplot2::geom_line() +
  ecodata::theme_facet() +
  ggplot2::facet_wrap(~Collcom) +
  ggplot2::ggtitle("GB Macro")+
  ggplot2::theme(legend.position = "bottom")

macrobensppwt |>
  dplyr::filter(EPU=="MAB") |>
  ggplot2::ggplot(ggplot2::aes(x=year, y=pyprop, colour = season_ng)) +
  ggplot2::geom_line() +
  ecodata::theme_facet() +
  ggplot2::facet_wrap(~Collcom) +
  ggplot2::ggtitle("MAB Macro")+
  ggplot2::theme(legend.position = "bottom")

Megabenthos composition over time in stomachs

NEFSC data only

megabenall_stn <- readRDS(here("fhdata/megabenall_stn.rds"))


megaben.analcoll <- megaben |>
  dplyr::select(PYNAM, PYCOMNAM, AnalSci, AnalCom, CollSci, Collcom)

collcolors <- data.frame(prey = unique(megaben.analcoll$Collcom),
                        #preycol = colors()[seq(1, length(colors()), 15)])
                        preycol = viridis::turbo(length(unique(megaben.analcoll$Collcom))))


megabensppwt <- megabenall_stn %>%
  dplyr::left_join(EPUlook) %>%
  dplyr::filter(megabenthos=="megaben",
                !is.na(EPU),
                !is.na(season_ng)) %>%
  dplyr::group_by(year, season_ng, EPU, megabenpynam) %>%
  dplyr::summarise(pysum = sum(megabenpywt)) %>%
  dplyr::mutate(pyprop = pysum / sum(pysum)) %>%
  dplyr::ungroup() %>%
  dplyr::left_join(megaben.analcoll, by =  dplyr::join_by(megabenpynam == PYNAM)) |>
  dplyr::left_join(collcolors, by =  dplyr::join_by(Collcom == prey))

preycols <- megabensppwt$preycol
names(preycols) <- as.factor(megabensppwt$megabenpynam)

preyaggcols <- megabensppwt$preycol
names(preyaggcols) <- as.factor(megabensppwt$Collcom)


top20spp <- megabensppwt |>
  dplyr::group_by(season_ng, EPU, Collcom) |>
  dplyr::summarise(meanprop = mean(pyprop, na.rm = TRUE),
                   maxprop = max(pyprop, na.rm = TRUE)) |>
  dplyr::ungroup() |>
  dplyr::arrange(EPU, season_ng, desc(meanprop))


pall <- plotDietCompBar(dat=megabensppwt, 
                          metric="pysum", 
                        catname="megabenpynam",
                        plotcol=preycols,
                          title="Diet composition")  

pcat <- plotDietCompBar(dat=megabensppwt, 
                          metric="pysum", 
                        catname="Collcom",
                        plotcol=preyaggcols,
                          title="Diet composition")

addSmallLegend(pall)

addSmallLegend(pcat)

Collcom trends by EPU

megabensppwt |>
  dplyr::filter(EPU=="GOM") |>
  ggplot2::ggplot(ggplot2::aes(x=year, y=pyprop, colour = season_ng)) +
  ggplot2::geom_line() +
  ecodata::theme_facet() +
  ggplot2::facet_wrap(~Collcom) +
  ggplot2::ggtitle("GOM Mega") +
  ggplot2::theme(legend.position = "bottom")

megabensppwt |>
  dplyr::filter(EPU=="GB") |>
  ggplot2::ggplot(ggplot2::aes(x=year, y=pyprop, colour = season_ng)) +
  ggplot2::geom_line() +
  ecodata::theme_facet() +
  ggplot2::facet_wrap(~Collcom) +
  ggplot2::ggtitle("GB Mega")+
  ggplot2::theme(legend.position = "bottom")

megabensppwt |>
  dplyr::filter(EPU=="MAB") |>
  ggplot2::ggplot(ggplot2::aes(x=year, y=pyprop, colour = season_ng)) +
  ggplot2::geom_line() +
  ecodata::theme_facet() +
  ggplot2::facet_wrap(~Collcom) +
  ggplot2::ggtitle("MAB Mega")+
  ggplot2::theme(legend.position = "bottom")

Have the benthivores used as samplers changed over time?

Any other sources to compare with?

NEFSC benthic database available here: https://www.obis.org/dataset/e7c86904-aac7-4a17-a895-99a54c430d80

Full dataset in local data-raw/NEFSCbenthicdatabase folder