-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathWHAMinputs.R
126 lines (109 loc) · 5.15 KB
/
WHAMinputs.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
# originally used for bluefish RTA 2022
# modified for herring RTA 2024
# create csv for WHAM input
# long input best?
# going wide...
library(dplyr)
library(ggplot2)
library(tidyr)
WHAMinputs <- function(infile, strata=NULL, outfile) {
splitoutput <- read.csv(infile)
#defaults to standard zoop model strata
if(is.null(strata)) {
stratlook <- data.frame(Stratum = c("Stratum_1",
"Stratum_2",
"Stratum_3",
"Stratum_4",
"Stratum_5",
"Stratum_6",
"Stratum_7"),
Region = c("AllEPU",
"her_sp",
"her_fa",
"MAB",
"GB",
"GOM",
"SS"))
} else {
stratlook <- strata
}
# for herring model, do herring_spring, herring_fall, and herring_larval areas
zoopindex <- splitoutput |>
dplyr::left_join(stratlook) |>
dplyr::select(Time, Region, Estimate, SE=Std..Error.for.Estimate) |>
tidyr::pivot_longer(c(Estimate, SE), names_to = "Var") |>
dplyr::group_by(Var) |>
tidyr::pivot_wider(names_from = Region, values_from = value) |>
dplyr::select(Time, dplyr::starts_with("her_")) |>
tidyr::pivot_longer(!c(Time, Var), names_to = "Region", values_to = "value") |>
tidyr::pivot_wider(names_from = "Var", values_from = "value") |>
tidyr::pivot_wider(names_from = "Region", values_from = c("Estimate", "SE"),
names_glue = "{Region}_{.value}", names_vary = "slowest")
readr::write_csv(zoopindex, outfile)
}
stratlook2 <- data.frame(Stratum = c("Stratum_1",
"Stratum_2",
"Stratum_3",
"Stratum_4",
"Stratum_5",
"Stratum_6",
"Stratum_7",
"Stratum_8",
"Stratum_9"),
Region = c("AllEPU",
"her_sp",
"her_fa",
"her_larv",
"no_larv",
"MAB",
"GB",
"GOM",
"SS"))
# make data files
# fall small copepod doy models did not converge, use without covariates
WHAMinputs(infile = "pyindex/smallcopeALL_fall_500_biascorrect/Index.csv",
outfile = "WHAMfits/fallsmallcopeALLindex.csv")
# spring large copepod models with doy covariate coverged and doy improved fit
WHAMinputs(infile = "pyindex/lgcopeALL_spring_500_biascorrect_doy/Index.csv",
outfile = "WHAMfits/springlargecopeindex.csv")
# didn't try covariates in the small copeopods larval area model, outa time
WHAMinputs(infile = "pyindex/smallcopeALL_sepfeb_yrshift_500_larvarea_biascorrect/Index.csv",
strata = stratlook2,
outfile = "WHAMfits/sepfebsmallcopeALLlarvareaindex.csv")
# tried covariates in the zoopvolume index and none converged
WHAMinputs(infile = "pyindex/zoopvol_spring_500_biascorrect/Index.csv",
outfile = "WHAMfits/springzoopvolumeindex.csv")
# bias corrected fall results
#splitoutput <- read.csv("pyindex/allagg_fall_500_lennosst_ALLsplit_biascorrect/Index.csv")
# # code for visualizing summed SE approximation vs calculated SE
# in2off <- splitoutput %>%
# left_join(stratlook) %>%
# dplyr::select(Time, Region, Estimate, Std..Error.for.Estimate) %>%
# tidyr::pivot_longer(c(Estimate, Std..Error.for.Estimate), names_to = "Var") %>%
# dplyr::group_by(Var) %>%
# tidyr::pivot_wider(names_from = Region, values_from = value) %>%
# dplyr::mutate(AlbInshore = MABGBalbinshore,
# NewOffshore = bfoff,
# BigOld = bfin,
# BigNew = bfall,
# CompBigNew = bfin + bfoff,
# AlbOld = AlbInshore + BigOld,
# AlbNew = AlbInshore + BigNew,
# StateWaters = MABGBstate,
# FedWaters = MABGBfed) %>%
# dplyr::select(Time, AlbInshore, NewOffshore, BigOld, BigNew, CompBigNew, AlbOld, AlbNew, StateWaters, FedWaters) %>%
# tidyr::pivot_longer(!c(Time, Var), names_to = "Region", values_to = "value") %>%
# tidyr::pivot_wider(names_from = "Var", values_from = "value")
#
# ggplot(in2off, aes(x=Time, y=Estimate, colour = Region)) +
# geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
# geom_point()+
# geom_line()
#
# compare <- in2off %>%
# filter(Region %in% c("CompBigNew", "BigNew"))
#
# ggplot(compare, aes(x=Time, y=Estimate, colour = Region)) +
# geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
# geom_point()+
# geom_line()