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