R/data_calculations.R
In ahart1-r-packages/flukeclimatemse: Run Climate Fluke MSE
# Data calculations for the summer flounder MSE and data source locations
# To Do:
# Search functions for "Data_formatting.R" in other R scripts in the package, should now reference this script
# update #Letter labels for manuscript tables/equations
##### OM truestocksize.R ##############################################################################
#Calculate weight-at-age as a fixed input to the stock projection. From Sinatra.CTL file provided by Dr. Gavin Fay in personal communication
#Mean Weights at age at start of year by stock and sex (fem,mal)
tempF1w <- c(0.001, 0.0140, 0.1069, 0.3142, 0.6284, 1.0216, 1.4611, 1.9182)
tempM1w <- c(0.001, 0.0064, 0.0491, 0.1445, 0.2889, 0.4697, 0.6718, 0.8819)
#Mean Weights at age in middle of year by stock and sex (fem,mal)
tempF2w <- c(0.0016, 0.0471, 0.1960, 0.4595, 0.8172, 1.2375, 1.6891, 2.1458)
tempM2w <- c(0.0007, 0.0216, 0.0901, 0.2112, 0.3757, 0.5689, 0.7765, 0.9865)
w_at_age <- rbind(tempF1w, tempF2w, tempM1w, tempM2w)
w_at_age <- colMeans(w_at_age) # Use mean weight at age for both sexes & start/middle of year measurements
print(w_at_age)
# Mortality at age (0-7+)
M_age <- c(0.26, 0.26, 0.26, 0.25, 0.25, 0.25, 0.25, 0.24) # ages 0-7+ from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A90, page 241, https://repository.library.noaa.gov/view/noaa/23031
## Approximate estimates for selectivity-at-age (age1-8) by fleet pulled from figures in the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) pages 378-381 figures A139-A142, then averaged all values, https://repository.library.noaa.gov/view/noaa/23031
# Combined commercial landings and discards selectivity averaged over different time periods
CommercialSelectivityMatrix <- matrix(c(0.01,0.18,1,1,0.9,0.8,0.7,0.9,
0.01,0.1,0.6,1,1,0.9,0.85,0.7,
0.01,0.05,0.25,0.6,1,1,1,0.75,
0.19,1,1,0.3,0.1,0.1,0.1,0.1,
0.1,0.3,0.85,1,0.85,0.7,0.7,0.8,
0.08,0.2,0.6,1,1,0.8,0.75,0.7), nrow=6, ncol=8, byrow = TRUE)
colnames(CommercialSelectivityMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(CommercialSelectivityMatrix) <- c("Land1982", "Land1995", "Land2008", "Disc1982", "Disc1995", "Disc2008")
CommercialSelectivity <- colMeans(CommercialSelectivityMatrix)
print(CommercialSelectivity) # Use as commercial selectivity
# Recreational catch selectivity averaged over different time periods
RecreationalCatchSelMatrix <- matrix(c(0.05,0.5,1,0.7,0.65,0.6,0.65,0.85,
0.01,0.1,0.6,1,0.98,0.75,0.75,0.7,
0.01,0.05,0.1,0.5,1,1,1,0.7), nrow=3, ncol=8, byrow = TRUE)
colnames(RecreationalCatchSelMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(RecreationalCatchSelMatrix) <- c("Land1982", "Land1995", "Land2008")
RecreationalCatchSelectivity <- colMeans(RecreationalCatchSelMatrix)
print(RecreationalCatchSelectivity) # Use as recreational catch selectivity
# Recreational discards selectivity averaged over different time periods
RecreationalDiscMatrix <- matrix(c(0.1,1,0.1,0.08,0.1,0.1,0.1,0.1,
0.1,0.7,1,0.6,0.2,0.25,0.3,0.35,
0.1,0.5,1,1,0.98,0.5,0.4,0.35), nrow=3, ncol=8, byrow = TRUE)
colnames(RecreationalDiscMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(RecreationalDiscMatrix) <- c("Disc1982", "Disc1995", "Disc2008")
RecreationalDiscSelectivity <- colMeans(RecreationalDiscMatrix)
print(RecreationalDiscSelectivity) # Use as recreational discards selectivity
# R_rho = recruitment autocorrelation
RData <- c(81955, 102427, 46954, 78263, 81397, 53988, 12474, 36963, 44019, 47704, 47264, 43928, 58403, 78348, 59520, 52374, 54518, 44100, 60551, 64979, 67860, 50131, 71270, 40634, 48153, 52646, 62460, 73747, 51331, 31296, 35187, 36719, 42271, 29833, 35853, 42415) # Recruitment timeseries (000s) from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A87, page 238, years 1982-2017, https://repository.library.noaa.gov/view/noaa/23031
MeanR <- mean(RData)
R_residuals <- RData-MeanR # calculate mean-observed
plot(R_residuals)
acf(R_residuals, lag.max = 1, plot = FALSE) # Calculate correltation: Use autocorrelation of Lag1 see additional explanation: https://stats.stackexchange.com/questions/81754/understanding-this-acf-output
# Result to use 0.414
##### OM dosurveys.R ##############################################################################################
# Calculate average NEFSC trawl survey CV using survey CV data (reported in percent) from the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019), Tables A42 & A44, page 187,189, https://repository.library.noaa.gov/view/noaa/23031
SpringCVs <- c(33,16,19,23,15,16,20,15,29,22,16,15,23,20,22,17,18,18,15,21,26,15,21,22,15,13,16,11,22,20,18,26,15,12,11,16,14,14,17,22,12,12)
FallCVs <- c(25,13,27,17,18,15,10,19,11,14,17,12,12,14,15,16,20,14,13,18,21,18,19,18,14,29,16,20,17,21,24,15,18,19,17)
AvgCV <- mean(c(SpringCVs, FallCVs)) # average CV reported as a percent
useAvgCV <- round(AvgCV/100, 3) # Use averageCV reported as a proportion
print(useAvgCV)
## Approximate estimates for selectivity-at-age (age1-8) by fleet pulled from figures in the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) pages 378-381 figures A139-A142, then averaged all values, https://repository.library.noaa.gov/view/noaa/23031
# Combined commercial landings and discards selectivity averaged over different time periods
CommercialSelectivityMatrix <- matrix(c(0.01,0.18,1,1,0.9,0.8,0.7,0.9,
0.01,0.1,0.6,1,1,0.9,0.85,0.7,
0.01,0.05,0.25,0.6,1,1,1,0.75,
0.19,1,1,0.3,0.1,0.1,0.1,0.1,
0.1,0.3,0.85,1,0.85,0.7,0.7,0.8,
0.08,0.2,0.6,1,1,0.8,0.75,0.7), nrow=6, ncol=8, byrow = TRUE)
colnames(CommercialSelectivityMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(CommercialSelectivityMatrix) <- c("Land1982", "Land1995", "Land2008", "Disc1982", "Disc1995", "Disc2008")
CommercialSelectivity <- colMeans(CommercialSelectivityMatrix)
# Recreational catch selectivity averaged over different time periods
RecreationalCatchSelMatrix <- matrix(c(0.05,0.5,1,0.7,0.65,0.6,0.65,0.85,
0.01,0.1,0.6,1,0.98,0.75,0.75,0.7,
0.01,0.05,0.1,0.5,1,1,1,0.7), nrow=3, ncol=8, byrow = TRUE)
colnames(RecreationalCatchSelMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(RecreationalCatchSelMatrix) <- c("Land1982", "Land1995", "Land2008")
RecreationalCatchSelectivity <- colMeans(RecreationalCatchSelMatrix)
# Recreational discards selectivity averaged over different time periods
RecreationalDiscMatrix <- matrix(c(0.1,1,0.1,0.08,0.1,0.1,0.1,0.1,
0.1,0.7,1,0.6,0.2,0.25,0.3,0.35,
0.1,0.5,1,1,0.98,0.5,0.4,0.35), nrow=3, ncol=8, byrow = TRUE)
colnames(RecreationalDiscMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(RecreationalDiscMatrix) <- c("Disc1982", "Disc1995", "Disc2008")
RecreationalDiscSelectivity <- colMeans(RecreationalDiscMatrix)
# Combined selectivity
AllSelectivity <- rbind(CommercialSelectivity, RecreationalCatchSelectivity, RecreationalDiscSelectivity)
avgSelectivity <- round(colMeans(AllSelectivity),2) # Use average selectivity estimated by assessment across commercial & recreational surveys, assume age 0 selectivity same as age 1 selectivity & age 7+ equal to age 7 estimated selectivity
print(avgSelectivity)
useSelectivity <- c(0.06, 0.06, 0.42, 0.66, 0.70, 0.70, 0.59, 0.58) # age 0-7+
##### Stockassess.R #################################################################################################
## Fixed data
W_dat <- 0.5 # ??? update to reflect summer flounder weight gain parameter
M_dat <- 0.25 # In Table #B ??? Reflect mean natural mortality across age and sex in in the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) page 65, https://repository.library.noaa.gov/view/noaa/23031
Rho <- 0.14 # 0.995 # update to reflect summer flounder von Bertalanffy k = growth rate coefficient = 0.14 combined across sexes from the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) page 58
## Parameter starting estimates
# Vector of biomass estimated in the first 2 years rather than assuming population is unfished at start of timeseries
AbundData <- read.csv("/Users/ahart2/Research/flukeclimatemse/data/Abund_at_age_A89_NEFSC66AssessmentReport.csv", header=TRUE)
rows <- AbundData$Year
AbundData <- AbundData[,2:9]
for(icol in 1:ncol(AbundData)){
AbundData[,icol] <- as.numeric(gsub(",","",AbundData[,icol]))
}
colnames(AbundData) <- c("Age0", "Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7plus")
rownames(AbundData) <- rows
Bstartvals <- AbundData %>%
rowSums() %>%
head(.,n=5) %>%
mean()
print(Bstartvals) # In Table #B ???, Average biomass over first 5 years (1982-1986) of timeseries used as starting estimate for first 2 years of biomass, from the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A89, page 240, https://repository.library.noaa.gov/view/noaa/23031
Bstartval <- rep(Bstartvals, 2)
# Fstartval is estimated as the unobserved fishing rate for recreational catch, discard and commercial fishing in the year before observed timeseries begins
Fstartval = 0.2
# Vector of fishing mortalities for years 1 to Nyear -1 accounting for total recreational catch & discards and commercial catch & discards
Fvals = rep(-2.0, (Nyear-1))
# Recruitment timeseries (000s) from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A87, page 238, years 1982-2017, https://repository.library.noaa.gov/view/noaa/23031
RData <- c(81955, 102427, 46954, 78263, 81397, 53988, 12474, 36963, 44019, 47704, 47264, 43928, 58403, 78348, 59520, 52374, 54518, 44100, 60551, 64979, 67860, 50131, 71270, 40634, 48153, 52646, 62460, 73747, 51331, 31296, 35187, 36719, 42271, 29833, 35853, 42415)
startingRMean <- mean(RData) # Use as starting parameter estimate when t=1, will update in subsequent years of the simulation as Recruitment data added
print(startingRMean)
# Annual recruitment deviation starting estimate: Exponentiated log values of recruitment deviations from assess_cor_file.csv
Rdevs <- c(exp(4.89E-01), exp(7.12E-01), exp(-6.81E-02), exp(4.43E-01), exp(4.82E-01), exp(7.17E-02), exp(-1.39E+00), exp(-3.07E-01), exp(-1.32E-01), exp(-5.21E-02), exp(-6.13E-02), exp(-1.34E-01), exp(1.51E-01), exp(4.45E-01), exp(1.71E-01), exp(4.30E-02), exp(8.36E-02), exp(-1.28E-01), exp(1.90E-01), exp(2.62E-01), exp(3.06E-01), exp(4.73E-03), exp(3.59E-01), exp(-2.01E-01), exp(-2.63E-02), exp(5.39E-02), exp(2.23E-01), exp(3.88E-01), exp(2.63E-02), exp(-4.70E-01), exp(-3.55E-01), exp(-3.12E-01), exp(-1.79E-01), exp(-5.34E-01), exp(-3.68E-01), exp(-1.83E-01))
meanRdevs <- mean(Rdevs)
sigmaR <- meanRdevs
print(sigmaR) # Use as starting parameter estimate for recruitment standard deviaiton
# ??? Recruitment random effect vector, normally distributed expect bounds between -2 and 2 so bounding at extremes of -5 ot 5 appropriate
R_randEffect = rep(0.001, Nyear)
##### Management model allocatequota.R ##############################################################################################
# Calculate proportionate participation in recreational fishery during reference period 1980-1989, data provided in personal communication by Dr. Jason McNamee
histPercent <- c(5.49, 5.66, 3.75, 17.63, 39.09, 3.14, 2.95, 16.69, 5.6)
names(histPercent) <- c("MA", "RI", "CT", "NY", "NJ", "DE", "MD", "VA", "NC")
histProp <- histPercent/sum(histPercent)
print(histProp)
# 9/30/20 Try to recreate the above, if successful repeat for most recent 10 years, lack data
library(tidyverse)
recCat <- as_tibble(read.csv("/Users/ahart2/Downloads/RecreationalCatch.csv"))
summary(recCat)
recCat %>%
rowwise() %>%
mutate(MAprop = MA/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
RIprop = RI/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
CTprop = CT/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
NYprop = NY/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
NJprop = NJ/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
DEprop = DE/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
MDprop = MD/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
VAprop = VA/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
NCprop = NC/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC)) %>%
filter(X > 1979) %>%
filter (X < 1990)
# ?????? Will need
# Convert weight-length & vice versa:
# in the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) page 59, https://repository.library.noaa.gov/view/noaa/23031
ahart1-r-packages/flukeclimatemse documentation built on Oct. 6, 2020, 2:48 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Data calculations for the summer flounder MSE and data source locations
# To Do:
# Search functions for "Data_formatting.R" in other R scripts in the package, should now reference this script
# update #Letter labels for manuscript tables/equations
##### OM truestocksize.R ##############################################################################
#Calculate weight-at-age as a fixed input to the stock projection. From Sinatra.CTL file provided by Dr. Gavin Fay in personal communication
#Mean Weights at age at start of year by stock and sex (fem,mal)
tempF1w <- c(0.001, 0.0140, 0.1069, 0.3142, 0.6284, 1.0216, 1.4611, 1.9182)
tempM1w <- c(0.001, 0.0064, 0.0491, 0.1445, 0.2889, 0.4697, 0.6718, 0.8819)
#Mean Weights at age in middle of year by stock and sex (fem,mal)
tempF2w <- c(0.0016, 0.0471, 0.1960, 0.4595, 0.8172, 1.2375, 1.6891, 2.1458)
tempM2w <- c(0.0007, 0.0216, 0.0901, 0.2112, 0.3757, 0.5689, 0.7765, 0.9865)
w_at_age <- rbind(tempF1w, tempF2w, tempM1w, tempM2w)
w_at_age <- colMeans(w_at_age) # Use mean weight at age for both sexes & start/middle of year measurements
print(w_at_age)
# Mortality at age (0-7+)
M_age <- c(0.26, 0.26, 0.26, 0.25, 0.25, 0.25, 0.25, 0.24) # ages 0-7+ from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A90, page 241, https://repository.library.noaa.gov/view/noaa/23031
## Approximate estimates for selectivity-at-age (age1-8) by fleet pulled from figures in the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) pages 378-381 figures A139-A142, then averaged all values, https://repository.library.noaa.gov/view/noaa/23031
# Combined commercial landings and discards selectivity averaged over different time periods
CommercialSelectivityMatrix <- matrix(c(0.01,0.18,1,1,0.9,0.8,0.7,0.9,
0.01,0.1,0.6,1,1,0.9,0.85,0.7,
0.01,0.05,0.25,0.6,1,1,1,0.75,
0.19,1,1,0.3,0.1,0.1,0.1,0.1,
0.1,0.3,0.85,1,0.85,0.7,0.7,0.8,
0.08,0.2,0.6,1,1,0.8,0.75,0.7), nrow=6, ncol=8, byrow = TRUE)
colnames(CommercialSelectivityMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(CommercialSelectivityMatrix) <- c("Land1982", "Land1995", "Land2008", "Disc1982", "Disc1995", "Disc2008")
CommercialSelectivity <- colMeans(CommercialSelectivityMatrix)
print(CommercialSelectivity) # Use as commercial selectivity
# Recreational catch selectivity averaged over different time periods
RecreationalCatchSelMatrix <- matrix(c(0.05,0.5,1,0.7,0.65,0.6,0.65,0.85,
0.01,0.1,0.6,1,0.98,0.75,0.75,0.7,
0.01,0.05,0.1,0.5,1,1,1,0.7), nrow=3, ncol=8, byrow = TRUE)
colnames(RecreationalCatchSelMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(RecreationalCatchSelMatrix) <- c("Land1982", "Land1995", "Land2008")
RecreationalCatchSelectivity <- colMeans(RecreationalCatchSelMatrix)
print(RecreationalCatchSelectivity) # Use as recreational catch selectivity
# Recreational discards selectivity averaged over different time periods
RecreationalDiscMatrix <- matrix(c(0.1,1,0.1,0.08,0.1,0.1,0.1,0.1,
0.1,0.7,1,0.6,0.2,0.25,0.3,0.35,
0.1,0.5,1,1,0.98,0.5,0.4,0.35), nrow=3, ncol=8, byrow = TRUE)
colnames(RecreationalDiscMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(RecreationalDiscMatrix) <- c("Disc1982", "Disc1995", "Disc2008")
RecreationalDiscSelectivity <- colMeans(RecreationalDiscMatrix)
print(RecreationalDiscSelectivity) # Use as recreational discards selectivity
# R_rho = recruitment autocorrelation
RData <- c(81955, 102427, 46954, 78263, 81397, 53988, 12474, 36963, 44019, 47704, 47264, 43928, 58403, 78348, 59520, 52374, 54518, 44100, 60551, 64979, 67860, 50131, 71270, 40634, 48153, 52646, 62460, 73747, 51331, 31296, 35187, 36719, 42271, 29833, 35853, 42415) # Recruitment timeseries (000s) from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A87, page 238, years 1982-2017, https://repository.library.noaa.gov/view/noaa/23031
MeanR <- mean(RData)
R_residuals <- RData-MeanR # calculate mean-observed
plot(R_residuals)
acf(R_residuals, lag.max = 1, plot = FALSE) # Calculate correltation: Use autocorrelation of Lag1 see additional explanation: https://stats.stackexchange.com/questions/81754/understanding-this-acf-output
# Result to use 0.414
##### OM dosurveys.R ##############################################################################################
# Calculate average NEFSC trawl survey CV using survey CV data (reported in percent) from the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019), Tables A42 & A44, page 187,189, https://repository.library.noaa.gov/view/noaa/23031
SpringCVs <- c(33,16,19,23,15,16,20,15,29,22,16,15,23,20,22,17,18,18,15,21,26,15,21,22,15,13,16,11,22,20,18,26,15,12,11,16,14,14,17,22,12,12)
FallCVs <- c(25,13,27,17,18,15,10,19,11,14,17,12,12,14,15,16,20,14,13,18,21,18,19,18,14,29,16,20,17,21,24,15,18,19,17)
AvgCV <- mean(c(SpringCVs, FallCVs)) # average CV reported as a percent
useAvgCV <- round(AvgCV/100, 3) # Use averageCV reported as a proportion
print(useAvgCV)
## Approximate estimates for selectivity-at-age (age1-8) by fleet pulled from figures in the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) pages 378-381 figures A139-A142, then averaged all values, https://repository.library.noaa.gov/view/noaa/23031
# Combined commercial landings and discards selectivity averaged over different time periods
CommercialSelectivityMatrix <- matrix(c(0.01,0.18,1,1,0.9,0.8,0.7,0.9,
0.01,0.1,0.6,1,1,0.9,0.85,0.7,
0.01,0.05,0.25,0.6,1,1,1,0.75,
0.19,1,1,0.3,0.1,0.1,0.1,0.1,
0.1,0.3,0.85,1,0.85,0.7,0.7,0.8,
0.08,0.2,0.6,1,1,0.8,0.75,0.7), nrow=6, ncol=8, byrow = TRUE)
colnames(CommercialSelectivityMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(CommercialSelectivityMatrix) <- c("Land1982", "Land1995", "Land2008", "Disc1982", "Disc1995", "Disc2008")
CommercialSelectivity <- colMeans(CommercialSelectivityMatrix)
# Recreational catch selectivity averaged over different time periods
RecreationalCatchSelMatrix <- matrix(c(0.05,0.5,1,0.7,0.65,0.6,0.65,0.85,
0.01,0.1,0.6,1,0.98,0.75,0.75,0.7,
0.01,0.05,0.1,0.5,1,1,1,0.7), nrow=3, ncol=8, byrow = TRUE)
colnames(RecreationalCatchSelMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(RecreationalCatchSelMatrix) <- c("Land1982", "Land1995", "Land2008")
RecreationalCatchSelectivity <- colMeans(RecreationalCatchSelMatrix)
# Recreational discards selectivity averaged over different time periods
RecreationalDiscMatrix <- matrix(c(0.1,1,0.1,0.08,0.1,0.1,0.1,0.1,
0.1,0.7,1,0.6,0.2,0.25,0.3,0.35,
0.1,0.5,1,1,0.98,0.5,0.4,0.35), nrow=3, ncol=8, byrow = TRUE)
colnames(RecreationalDiscMatrix) <- c("Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7", "Age8")
rownames(RecreationalDiscMatrix) <- c("Disc1982", "Disc1995", "Disc2008")
RecreationalDiscSelectivity <- colMeans(RecreationalDiscMatrix)
# Combined selectivity
AllSelectivity <- rbind(CommercialSelectivity, RecreationalCatchSelectivity, RecreationalDiscSelectivity)
avgSelectivity <- round(colMeans(AllSelectivity),2) # Use average selectivity estimated by assessment across commercial & recreational surveys, assume age 0 selectivity same as age 1 selectivity & age 7+ equal to age 7 estimated selectivity
print(avgSelectivity)
useSelectivity <- c(0.06, 0.06, 0.42, 0.66, 0.70, 0.70, 0.59, 0.58) # age 0-7+
##### Stockassess.R #################################################################################################
## Fixed data
W_dat <- 0.5 # ??? update to reflect summer flounder weight gain parameter
M_dat <- 0.25 # In Table #B ??? Reflect mean natural mortality across age and sex in in the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) page 65, https://repository.library.noaa.gov/view/noaa/23031
Rho <- 0.14 # 0.995 # update to reflect summer flounder von Bertalanffy k = growth rate coefficient = 0.14 combined across sexes from the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) page 58
## Parameter starting estimates
# Vector of biomass estimated in the first 2 years rather than assuming population is unfished at start of timeseries
AbundData <- read.csv("/Users/ahart2/Research/flukeclimatemse/data/Abund_at_age_A89_NEFSC66AssessmentReport.csv", header=TRUE)
rows <- AbundData$Year
AbundData <- AbundData[,2:9]
for(icol in 1:ncol(AbundData)){
AbundData[,icol] <- as.numeric(gsub(",","",AbundData[,icol]))
}
colnames(AbundData) <- c("Age0", "Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7plus")
rownames(AbundData) <- rows
Bstartvals <- AbundData %>%
rowSums() %>%
head(.,n=5) %>%
mean()
print(Bstartvals) # In Table #B ???, Average biomass over first 5 years (1982-1986) of timeseries used as starting estimate for first 2 years of biomass, from the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A89, page 240, https://repository.library.noaa.gov/view/noaa/23031
Bstartval <- rep(Bstartvals, 2)
# Fstartval is estimated as the unobserved fishing rate for recreational catch, discard and commercial fishing in the year before observed timeseries begins
Fstartval = 0.2
# Vector of fishing mortalities for years 1 to Nyear -1 accounting for total recreational catch & discards and commercial catch & discards
Fvals = rep(-2.0, (Nyear-1))
# Recruitment timeseries (000s) from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A87, page 238, years 1982-2017, https://repository.library.noaa.gov/view/noaa/23031
RData <- c(81955, 102427, 46954, 78263, 81397, 53988, 12474, 36963, 44019, 47704, 47264, 43928, 58403, 78348, 59520, 52374, 54518, 44100, 60551, 64979, 67860, 50131, 71270, 40634, 48153, 52646, 62460, 73747, 51331, 31296, 35187, 36719, 42271, 29833, 35853, 42415)
startingRMean <- mean(RData) # Use as starting parameter estimate when t=1, will update in subsequent years of the simulation as Recruitment data added
print(startingRMean)
# Annual recruitment deviation starting estimate: Exponentiated log values of recruitment deviations from assess_cor_file.csv
Rdevs <- c(exp(4.89E-01), exp(7.12E-01), exp(-6.81E-02), exp(4.43E-01), exp(4.82E-01), exp(7.17E-02), exp(-1.39E+00), exp(-3.07E-01), exp(-1.32E-01), exp(-5.21E-02), exp(-6.13E-02), exp(-1.34E-01), exp(1.51E-01), exp(4.45E-01), exp(1.71E-01), exp(4.30E-02), exp(8.36E-02), exp(-1.28E-01), exp(1.90E-01), exp(2.62E-01), exp(3.06E-01), exp(4.73E-03), exp(3.59E-01), exp(-2.01E-01), exp(-2.63E-02), exp(5.39E-02), exp(2.23E-01), exp(3.88E-01), exp(2.63E-02), exp(-4.70E-01), exp(-3.55E-01), exp(-3.12E-01), exp(-1.79E-01), exp(-5.34E-01), exp(-3.68E-01), exp(-1.83E-01))
meanRdevs <- mean(Rdevs)
sigmaR <- meanRdevs
print(sigmaR) # Use as starting parameter estimate for recruitment standard deviaiton
# ??? Recruitment random effect vector, normally distributed expect bounds between -2 and 2 so bounding at extremes of -5 ot 5 appropriate
R_randEffect = rep(0.001, Nyear)
##### Management model allocatequota.R ##############################################################################################
# Calculate proportionate participation in recreational fishery during reference period 1980-1989, data provided in personal communication by Dr. Jason McNamee
histPercent <- c(5.49, 5.66, 3.75, 17.63, 39.09, 3.14, 2.95, 16.69, 5.6)
names(histPercent) <- c("MA", "RI", "CT", "NY", "NJ", "DE", "MD", "VA", "NC")
histProp <- histPercent/sum(histPercent)
print(histProp)
# 9/30/20 Try to recreate the above, if successful repeat for most recent 10 years, lack data
library(tidyverse)
recCat <- as_tibble(read.csv("/Users/ahart2/Downloads/RecreationalCatch.csv"))
summary(recCat)
recCat %>%
rowwise() %>%
mutate(MAprop = MA/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
RIprop = RI/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
CTprop = CT/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
NYprop = NY/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
NJprop = NJ/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
DEprop = DE/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
MDprop = MD/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
VAprop = VA/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC),
NCprop = NC/sum(MA, RI, CT, NY, NJ, DE, MD, VA, NC)) %>%
filter(X > 1979) %>%
filter (X < 1990)
# ?????? Will need
# Convert weight-length & vice versa:
# in the 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) page 59, https://repository.library.noaa.gov/view/noaa/23031
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.