more dat-file work

2023_runs/README.md

@@ -7,6 +7,7 @@ sep23 is the base selected change model from september which has:
   - spawning biomass weight-at-age set to the RE model estimates from fishery data for the A-season
   - the new AVO time series
   - modest process error in ATS selectivity variability
+  - Pete Hulson's publication on effective sample size for bottom-trawl survey data
 
 
 

2023_runs/dat/README.md

@@ -1,5 +1,24 @@
 # Notes
 
+## 2023 November assessment
+
+[x] Create fixed-effects fishery wt-age datafile "wtage2023.dat" (also done in script "wt.R" in WtAgeRE, but calls on specimen data from survey and sampler data from fishery) 
+[x] Generated random-effect SSB wt-at-age (done in script "wt.R" in WtAgeRE)
+[x] Expand Paul Spencer's input files to work for 2023 (but placeholders only)
+[ ] Put SSB wt-at-age into appropriate input datafile
+[ ] Update fishery catch
+[ ] Update fishery catch-at-age 
+[ ] Update survey cpue so that design-based estimator can be used from scratch
+[ ] Update BTS design based estimates
+[x] Update BTS Covariance matrix ("cov_2023.dat")
+[ ] Add in Hulson's BTS sample size
+[ ] Update AVO time series
+[ ] Update ATS age composition for the 2022 data (previously used BTS ALK)
+
+
+pm_22.dat: Full final dataset with
+
+## September work
 Changed "std_tune.dat" to "obs_err_init.dat" to accommodate adding in datafiles that have tuned stdevs (in directory "tune")
 
 Covariance matrix is derived from Stan's DD algorthm before VAST

2023_runs/dat/comp.R

@@ -0,0 +1,23 @@
+
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+library(ggridges)
+
+df <- data.frame(read.table("agew.dat"))
+names(df) <- 1:15
+df$Year <- (1982:2018)
+df.l <- gather(df,"Age","N",-Year) %>% mutate(Source="DB")
+df2 <- data.frame(read.table("age1.dat",header=TRUE))
+# put the two in one dataframe
+df <- rbind(df2,df.l) 
+
+# Compute proportions
+dfp <- df %>% group_by(Year,Source)  %>% mutate(Age=as.numeric(Age),Proportion=N/sum(N)) %>% filter(Year<=1999)
+
+# Plot
+ggplot(dfp,aes(x=Age,y=as.factor(Year),height = Proportion, fill=Source,color=Source)) + 
+     geom_density_ridges2(stat = "identity",scale = 3, alpha = .1)+
+     xlim(c(1,15)) + ylab("Year") + xlab("Age (years)") + scale_y_discrete(limits=rev(levels(as.factor(dfp$Year)))) + theme_few() 
+     geom_density_ridges(stat = "binline", bins=10,  scale = 0.95)+
+

2023_runs/dat/rescale_cov.R

@@ -0,0 +1,102 @@
+radian
+1
+install.packages("MBESS")
+install.packages("GGally")
+library(MBESS)
+library(GGally)
+#c1 <- as.matrix(read.table("dat/cov_biom.dat",header=F,as.is=T))
+#c1 <- as.matrix(read.table("cov22_orig.dat",header=F,as.is=T))
+#c1 <- as.matrix(read.csv("stan_cov_matrix_2022_db.csv",header=F,as.is=T))
+#c1 <- as.matrix(read.table("stan22.dat",header=F,as.is=T))
+c1 <- as.matrix(read.table("stan_2022_cov.dat",header=F,as.is=T))
+#c1 <- as.matrix(read.table("stan_cov_21.dat",header=F,as.is=T))
+#c1 <- as.matrix(read.table("stan_cov_21.dat",header=F,as.is=T))
+#c1 <- as.matrix(read.csv("varcov.csv",header=F,as.is=T))
+#c1 <- as.matrix(read.table("dat/cov_biom.dat",header=F,as.is=T))
+#c1 <- as.matrix(read.table("dat/rescaled_cov_2018.dat",header=F,as.is=T))
+#c1 <- as.matrix(read.table("dat/rescaled_cov_2018_fitsrv.dat",header=F,as.is=T))
+#c1 <- as.matrix(read.table("dat/thorCov.dat",header=F,as.is=T))
+
+##########################
+#c1 <- as.matrix(read.table("cov_biom19.dat",header=F,as.is=T))
+c1
+d1 <- sqrt(diag(c1))
+#d2 <- as.vector(d1/1e6) #kg
+d2 <- as.vector(d1/1e3)
+#d2 <- as.vector(d1)
+c2 <- cov2cor(c1)
+c3 <- sweep(sweep(c2, 1L, d2, "*"), 2L, d2, "*")
+all.equal(c3,c1)
+?write.table
+q
+c3[1:4,1:4 ]
+
+write.table(c3,file="cov_2022.dat",row.names=F,col.names=F)
+##########################
+
+ set up a model with recruitment set to 1, to generate yield per recruit and spawning biomass per recruit (spr).  That just uses selectivity, 
+ fishing rate, size-at-age, and maturity at age
+- Calculate equilibrium recruitment using 
+
+spr = 9.197e-7
+- based on that, I can get total yield (ypr times rinf) spawning biomss (spr * rinf), etc.
+alpha = exp(-4.43) 
+beta = exp(-6.42)
+ssb = seq(0,1000,100)
+plot(ssb,alpha*ssb*exp(-beta*ssb))
+
+# note here the 1e^ 6 is to reflect that the stock recruit function predicts millions of recruits.  I could have easily folded that into alpha.
+rinf = log(alpha * spr * 1E6) / (beta * spr)  
+rinf
+
+
+c1 <- as.matrix(read.table("dat/thorCov.dat",header=F))
+df <- as.numeric(read.table("dat/vast_biom.dat",header=F,as.is=T))
+names(c1)= 1:38
+names(df)= 1:37
+m <- as.data.frame(mvrnorm(n = 1000, df, c1, tol = 1e-6, empirical = FALSE, EISPACK = FALSE))
+colnames(m) <- 1982:2018
+ggpairs(m,columns=3:7)
+ggpairs(m,columns=33:37)
+pm <- ggpairs(tips, 2:3)
+data(tips, package = "reshape")
+pm <- ggpairs(tips[, 1:3])
+p_(pm)
+pm <- ggpairs(tips, 1:3, columnLabels = c("Total Bill", "Tip", "Sex"))
+p_(pm)
+d1 <- sqrt(diag(c1))
+d2 <- as.vector(d1/1000)
+c2 <- cov2cor(c1)
+c3 <- sweep(sweep(c2, 1L, d2, "*"), 2L, d2, "*")
+all.equal(c3,c1)
+write.table(c3,file="rescaled_cov_2019.dat")
+
+c1 <- as.matrix(read.csv("dat/thorCov_In.csv",header=F,as.is=T))
+c1 <- as.matrix(read.csv("Cov_Index.csv",header=F,as.is=T))
+d1 <- sqrt(diag(c1))
+d2 <- as.vector(d1/1000)
+c2 <- cov2cor(c1)
+c3 <- sweep(sweep(c2, 1L, d2, "*"), 2L, d2, "*")
+all.equal(c3,c1)
+write.table(c3,file="dat/thorCov.dat",row.names=F,col.names=F)
+d2 <- as.vector(d1/1000)
+n <- length(d2)
+d2[n] <- d2[n]*.5
+
+c2 <- cov2cor(c1)
+c3 <- sweep(sweep(c2, 1L, d2, "*"), 2L, d2, "*")
+all.equal(c3,c1)
+write.table(c3,file="dat/rescaled_cov_2018_fitsrv.dat")
+
+c1 <- as.matrix(read.csv("varcov.csv",header=F,as.is=T))
+d1 <- sqrt(diag(c1))
+d2 <- as.vector(d1/1000)
+n <- length(d2)
+d2[n] <- d2[n]*.5
+c2 <- cov2cor(c1)
+c3 <- sweep(sweep(c2, 1L, d2, "*"), 2L, d2, "*")
+all.equal(c3,c1)
+write.table(c3,file="dat/rescaled_cov_2018_fitsrv.dat")
+
+q()
+n

2023_runs/dat/selvar23.dat

@@ -0,0 +1,60 @@
+1964 0 0 0
+1965 0.5 0 0
+1966 0.5 0 0
+1967 0.5 0 0
+1968 0.5 0 0
+1969 0.5 0 0
+1970 0.5 0 0
+1971 0.5 0 0
+1972 0.5 0 0
+1973 0.5 0 0
+1974 0.5 0 0
+1975 0.5 0 0
+1976 0.5 0 0
+1977 0.5 0 0
+1978 0.5 0 0
+1979 0.5 0 0
+1980 0.5 0 0
+1981 0.5 0 0
+1982 0.5 0 0
+1983 0.5 0 0
+1984 0.5 0 0
+1985 0.5 0 0
+1986 0.5 0 0
+1987 0.5 0 0
+1988 0.5 0 0
+1989 0.5 0 0
+1990 0.5 0 0
+1991 0.5 0 0
+1992 0.5 0 0
+1993 0.5 0 0
+1994 0.5 0 0
+1995 0.5 0 0.138
+1996 0.5 0 0.138
+1997 0.5 0 0.138
+1998 0.5 0 0.138
+1999 0.5 0 0.138
+2000 0.5 0 0.138
+2001 0.5 0 0.138
+2002 0.5 0 0.138
+2003 0.5 0 0.138
+2004 0.5 0 0.138
+2005 0.5 0 0.138
+2006 0.5 0 0.138
+2007 0.5 0 0.138
+2008 0.5 0 0.138
+2009 0.5 0 0.138
+2010 0.5 0 0.138
+2011 0.5 0 0.138
+2012 0.5 0 0.138
+2013 0.5 0 0.138
+2014 0.5 0 0.138
+2015 0.5 0 0.138
+2016 0.5 0 0.138
+2017 0.5 0 0.138
+2018 0.5 0 0.138
+2019 1.9 0 0.138
+2020 1.9 0 0.138
+2021 0 0 0.138
+2022 0 0 0.138
+2023 0 0 0.138

2023_runs/dat/ssb23.csv

@@ -0,0 +1,57 @@
+year,3,4,5,6,7,8,9,10,11,12,13,14,15,source
+1970,0.311072,0.427458,0.55009,0.674562,0.797415,0.916081,1.028772,1.134344,1.232157,1.321963,1.403799,1.477905,1.544656,Predicted
+1971,0.307228,0.452313,0.576279,0.701144,0.82365,0.941422,1.052838,1.156889,1.253046,1.341142,1.421276,1.49373,1.558912,Predicted
+1972,0.321568,0.433307,0.585158,0.711117,0.834228,0.952199,1.063499,1.167202,1.262849,1.350331,1.429793,1.501553,1.566041,Predicted
+1973,0.313585,0.464989,0.584425,0.738544,0.862507,0.980459,1.091068,1.193594,1.287737,1.373516,1.451177,1.521113,1.58381,Predicted
+1974,0.329852,0.457189,0.6163,0.738007,0.890127,1.008925,1.119505,1.221329,1.314283,1.398545,1.47449,1.542613,1.603475,Predicted
+1975,0.298991,0.579313,0.720039,0.883094,1.001329,1.144476,1.250469,1.345787,1.430982,1.506773,1.573952,1.633328,1.685689,Predicted
+1976,0.332728,0.441983,0.729979,0.872966,1.034031,1.147123,1.28293,1.380174,1.465961,1.541319,1.607317,1.664999,1.71534,Predicted
+1977,0.294369,0.63373,0.759139,1.051895,1.190692,1.340931,1.438571,1.555963,1.633143,1.698221,1.752966,1.798972,1.837634,Predicted
+1978,0.324975,0.425699,0.772108,0.899593,1.190521,1.324595,1.468092,1.557697,1.666336,1.73448,1.790565,1.836587,1.874294,Predicted
+1979,0.32522,0.483367,0.592591,0.941505,1.066786,1.352017,1.477959,1.611767,1.690814,1.788555,1.845853,1.891417,1.92743,Predicted
+1980,0.309436,0.567516,0.738667,0.851723,1.197264,1.31383,1.586623,1.697742,1.815398,1.877776,1.958924,2.000129,2.030383,Predicted
+1981,0.29829,0.499527,0.767809,0.941965,1.052375,1.391079,1.497887,1.759051,1.857499,1.962077,2.011437,2.079959,2.109153,Predicted
+1982,0.28816,0.424556,0.632569,0.902848,1.075247,1.181115,1.513338,1.612421,1.865168,1.954929,2.05086,2.091833,2.152377,Predicted
+1983,0.338035,0.465312,0.611216,0.82203,1.089842,1.25587,1.352645,1.674029,1.761304,2.001863,2.079492,2.163657,2.193437,Predicted
+1984,0.333252,0.4394,0.572117,0.719624,0.929027,1.193194,1.354018,1.444591,1.759219,1.839519,2.073137,2.144034,2.221794,Predicted
+1985,0.369905,0.47498,0.588734,0.723692,0.869227,1.073532,1.330424,1.482577,1.563703,1.868579,1.939175,2.163378,2.22532,Predicted
+1986,0.304231,0.439473,0.548282,0.663136,0.797126,0.940158,1.140893,1.393528,1.541044,1.617383,1.917496,1.98347,2.203279,Predicted
+1987,0.349399,0.40203,0.542521,0.652877,0.76637,0.896842,1.034853,1.229605,1.475721,1.616508,1.68615,1.979767,2.039562,Predicted
+1988,0.321568,0.430801,0.487801,0.62958,0.738802,0.849367,0.97566,1.108692,1.298017,1.538533,1.673746,1.737981,2.026454,Predicted
+1989,0.275451,0.412535,0.52665,0.585088,0.725601,0.831551,0.937447,1.058175,1.185143,1.368209,1.602495,1.731666,1.790154,Predicted
+1990,0.28947,0.422648,0.567631,0.684074,0.740464,0.875682,0.974076,1.070966,1.181882,1.298723,1.47171,1.696219,1.816089,Predicted
+1991,0.303474,0.455308,0.597385,0.744991,0.859126,0.909551,1.036256,1.124504,1.21034,1.309847,1.415331,1.577302,1.791333,Predicted
+1992,0.340274,0.393124,0.549769,0.693264,0.839622,0.950532,0.996355,1.117575,1.199848,1.279516,1.372883,1.472413,1.62872,Predicted
+1993,0.392955,0.516439,0.578744,0.738174,0.879218,1.019239,1.121106,1.156151,1.265629,1.335782,1.403385,1.485052,1.57345,Predicted
+1994,0.361138,0.534251,0.665318,0.729857,0.887321,1.023282,1.156051,1.249273,1.2749,1.374656,1.435133,1.493352,1.566091,Predicted
+1995,0.316379,0.471433,0.650466,0.783277,0.846281,0.999777,1.130077,1.256098,1.341968,1.360007,1.45221,1.505361,1.556611,Predicted
+1996,0.333901,0.390836,0.549886,0.730096,0.861871,0.922196,1.071871,1.197615,1.318673,1.399421,1.41236,1.499618,1.548064,Predicted
+1997,0.354105,0.444891,0.507783,0.668588,0.847253,0.975036,1.029664,1.172548,1.290894,1.404315,1.477463,1.48303,1.563275,Predicted
+1998,0.309331,0.400946,0.494246,0.557878,0.718031,0.895011,1.020389,1.072152,1.211914,1.327037,1.437251,1.507287,1.509895,Predicted
+1999,0.329153,0.431508,0.529679,0.624911,0.686843,0.842601,1.01331,1.131213,1.174832,1.306188,1.412945,1.515044,1.57736,Predicted
+2000,0.300434,0.408751,0.515378,0.614808,0.708932,0.768001,0.919673,1.085512,1.19811,1.236252,1.362158,1.463627,1.560697,Predicted
+2001,0.337797,0.453002,0.569508,0.678547,0.775854,0.86449,0.915727,1.058065,1.213734,1.315836,1.34353,1.459302,1.551131,Predicted
+2002,0.331896,0.44944,0.570636,0.688907,0.796393,0.889684,0.972589,1.016996,1.151892,1.29988,1.394336,1.414615,1.523333,Predicted
+2003,0.345542,0.441977,0.565429,0.688366,0.805105,0.908631,0.996271,1.072442,1.109511,1.236834,1.377283,1.464428,1.477751,Predicted
+2004,0.333214,0.498176,0.602803,0.728668,0.849481,0.96073,1.056421,1.134723,1.200719,1.227287,1.344158,1.474469,1.551969,Predicted
+2005,0.277321,0.403093,0.571805,0.677536,0.802429,0.920728,1.02839,1.119806,1.19345,1.254639,1.276422,1.388651,1.514547,Predicted
+2006,0.307204,0.443658,0.578357,0.749699,0.853116,0.972025,1.081786,1.179272,1.2596,1.3218,1.371598,1.382333,1.484051,Predicted
+2007,0.306246,0.492451,0.638847,0.776474,0.945239,1.041992,1.151392,1.24982,1.334957,1.402541,1.452055,1.489549,1.488579,Predicted
+2008,0.305784,0.457111,0.651412,0.800194,0.935722,1.09906,1.188069,1.288239,1.376611,1.451368,1.50862,1.548115,1.576076,Predicted
+2009,0.351046,0.532902,0.696418,0.894309,1.03993,1.167289,1.318968,1.394083,1.479114,1.55186,1.611064,1.653231,1.678375,Predicted
+2010,0.269194,0.446664,0.633652,0.79868,0.995241,1.137422,1.259873,1.405702,1.474443,1.552896,1.619094,1.671947,1.708072,Predicted
+2011,0.246138,0.431514,0.617696,0.80725,0.970019,1.160741,1.29459,1.40711,1.54212,1.599693,1.66703,1.722446,1.765044,Predicted
+2012,0.277195,0.394309,0.587637,0.776162,0.963654,1.121093,1.30421,1.428994,1.531637,1.656452,1.703878,1.761374,1.807428,Predicted
+2013,0.344036,0.437337,0.563045,0.758905,0.945201,1.126933,1.276152,1.449471,1.563581,1.655206,1.769055,1.805845,1.853221,Predicted
+2014,0.410493,0.488791,0.58986,0.717857,0.911703,1.092792,1.267094,1.407456,1.571127,1.675277,1.756989,1.861224,1.888867,Predicted
+2015,0.373447,0.448184,0.528505,0.63017,0.757642,0.950132,1.129287,1.301282,1.439133,1.60021,1.701779,1.780988,1.882841,Predicted
+2016,0.29962,0.434905,0.512939,0.594232,0.695042,0.820304,1.009639,1.185034,1.352933,1.486555,1.643423,1.74091,1.816236,Predicted
+2017,0.335214,0.417055,0.558643,0.638534,0.718192,0.814778,0.934012,1.116163,1.283729,1.443549,1.569129,1.718197,1.808264,Predicted
+2018,0.357598,0.430726,0.517694,0.660791,0.739354,0.815576,0.907259,1.02065,1.196434,1.357429,1.510708,1.629944,1.772977,Predicted
+2019,0.433823,0.546216,0.629467,0.719417,0.85989,0.931668,0.998208,1.078352,1.179169,1.341976,1.490055,1.630805,1.738123,Predicted
+2020,0.374715,0.48507,0.600214,0.684275,0.773512,0.912141,0.981288,1.044693,1.121421,1.218713,1.37801,1.522685,1.660197,Predicted
+2021,0.302947,0.426239,0.539358,0.655317,0.738661,0.826045,0.962029,1.028024,1.087995,1.161178,1.254941,1.410816,1.552236,Predicted
+2022,0.340108,0.446412,0.577403,0.692791,0.806753,0.884937,0.964956,1.092163,1.148595,1.198696,1.262054,1.346289,1.493099,Predicted
+2023,0.319566,0.490613,0.604994,0.738365,0.851658,0.960207,1.030665,1.101476,1.218651,1.264729,1.304523,1.357884,1.432609,Predicted
+2024,0.313754,0.449219,0.627224,0.743655,0.87522200000000006,0.983852,1.085745,1.148271,1.21044,1.318695,1.355893,1.387076,1.432245,Predicted
+2025,0.313754,0.443407,0.58583,0.765885,0.880512,1.007415,1.10939,1.203351,1.257234,1.310483,1.409859,1.438446,1.461437,Predicted