Commit a5193b7f authored by bberges's avatar bberges
Browse files

commit update to scripts

parent e58995d3
This diff is collapsed.
# script that derives the HERAS index directly from the StoX projects.
rm(list=ls())
library(Rstox)
library(FLFleet)
library(ggplotFL)
library(mgcv)
library(tidyverse)
library(ggplot2)
#Set up directories
#path <- 'J:/git/HERAS/'
path <- 'J:/git/heras_index_kbwot'
try(setwd(path),silent=TRUE)
mainPath <- file.path(".")
dataPath <- file.path(".","data")
stoxPath <- file.path(".","StoX")
outPath <- file.path(".","output")
functionPath <- file.path(".","functions")
resultsPath <- file.path(".","results")
figurePath <- file.path(".","figures",'kbwot','2_baseline')
load(file.path(dataPath,'miscellaneous','NSH_HAWG2020_sf.Rdata'))
load(file.path(resultsPath,'HERAS_baseline.RData'))
load(file.path(resultsPath,'HERAS_baseline_bootstrap.RData'))
SSB_WBSS_assessment <- read.csv(file.path(dataPath,'WBSS','assessment_2020.csv'))
HERAS_WBSS <- read.csv(file.path(dataPath,'WBSS','WBSS_HERAS.csv'),check.names=FALSE)
wa_WBSS <- read.csv(file.path(dataPath,'WBSS','WBSS_wa.csv'),check.names=FALSE)
mat_WBSS <- read.csv(file.path(dataPath,'WBSS','WBSS_mProp.csv'),check.names=FALSE)
################################################################
# compute df
################################################################
################
# NSAS
################
NSAS_df <- as.data.frame(HERAS_baseline$HERAS.NSAS$HERAS_strata@index[,ac(2017:2020),'abu','all'])
temp <- as.data.frame(HERAS_baseline$HERAS.NSAS$HERAS_all@index[,ac(2017:2020),'abu','all'])
NSAS_df <- rbind(NSAS_df,temp)
NSAS_df$age <- as.factor(NSAS_df$age)
NSAS_df <- NSAS_df[NSAS_df$data != 0,]
# NSAS N matyrity df
NSAS_mat_df <- as.data.frame(HERAS_baseline$HERAS.NSAS$HERAS_strata@index[,ac(2017:2020),'abu',c('IMM','MAT')])
temp <- as.data.frame(HERAS_baseline$HERAS.NSAS$HERAS_all@index[,ac(2017:2020),'abu',c('IMM','MAT')])
NSAS_mat_df <- rbind(NSAS_mat_df,temp)
NSAS_mat_df$age <- as.factor(NSAS_mat_df$age)
NSAS_mat_df <- NSAS_mat_df[NSAS_df$data != 0,]
# NSAS N df
NSAS_N_bootstrap <- as.data.frame(HERAS_baseline_bootstrap$stk.NSAS@stock.n[,ac(2017:2020)])
NSAS_N_bootstrap$age <- as.factor(NSAS_N_bootstrap$age)
NSAS_N <- as.data.frame(HERAS_baseline$stk.NSAS@stock.n[,ac(2017:2020)])
NSAS_N$age <- as.factor(NSAS_N$age)
# NSAS SSB df
SSB_NSAS_assessment <- ssb(NSH.sam)
SSB_NSAS_HERAS <- apply(NSH.tun$HERAS@index*
NSH@stock.wt[dimnames(NSH.tun$HERAS@index)$age,dimnames(NSH.tun$HERAS@index)$year]*
NSH@mat[dimnames(NSH.tun$HERAS@index)$age,dimnames(NSH.tun$HERAS@index)$year],2,sum)
SSB_NSAS_HERAS <- as.data.frame(SSB_NSAS_HERAS)
SSB_NSAS_HERAS_new <- as.data.frame(HERAS_baseline$stk.NSAS@stock[,ac(2017:2020)])
SSB_NSAS_HERAS <- rbind(SSB_NSAS_HERAS,SSB_NSAS_HERAS_new)
temp <- as.data.frame(HERAS_baseline_bootstrap$stk.NSAS@stock[,ac(2017:2020)])
temp$iter <- as.numeric(temp$iter)
SSB_NSAS_HERAS_boot <- temp %>%
group_by(year) %>%
summarise_at( vars(data),
list( Q25=~quantile(., probs = 0.25),
Q50=~quantile(., probs = 0.50),
Q75=~quantile(., probs = 0.75)))
# NSAS % difference df
stockN_plus <- HERAS_baseline$stk.NSAS@stock.n[,ac(2017:2019)]*1e3
stockN_plus[ac(8)] <- apply(stockN_plus[ac(8:9)],2,sum)
A <- apply( NSH.tun$HERAS@index[ac(1:8),ac(2017:2019)]*
NSH@stock.wt[ac(1:8),ac(2017:2019)]*
NSH@mat[ac(1:8),ac(2017:2019)],2,sum)
B <- apply( HERAS_baseline$stk.NSAS@stock.n[ac(1:9),ac(2017:2019)]*
HERAS_baseline$stk.NSAS@stock.wt[ac(1:9),ac(2017:2019)]*
HERAS_baseline$stk.NSAS@mat[ac(1:9),ac(2017:2019)],2,sum)
diff_NSAS <- as.data.frame((NSH.tun$HERAS@index[,ac(2017:2019)]-
stockN_plus[ac(1:8)])/
NSH.tun$HERAS@index[,ac(2017:2019)])
diff_NSAS <- rbind(diff_NSAS,as.data.frame((A-B)/A))
diff_NSAS$age[diff_NSAS$age == 'all'] <- 'SSB'
diff_NSAS$year <- as.factor(diff_NSAS$year)
# NSAS bio df
A <- as.data.frame(HERAS_baseline$stk.NSAS@stock.wt[ac(0:8),ac(2017:2020)]*1e-3)
A$source <- 'new'
A$type <- 'weight'
B <- as.data.frame(NSH@stock.wt[ac(0:8),ac(2017:2019)])
B$source <- 'assessment'
B$type <- 'weight'
C <- as.data.frame(HERAS_baseline$stk.NSAS@mat[ac(0:8),ac(2017:2020)])
C$source <- 'new'
C$type <- 'maturity'
D <- as.data.frame(NSH@mat[ac(0:8),ac(2017:2019)])
D$source <- 'assessment'
D$type <- 'maturity'
NSAS_bio <- rbind(A,B,C,D)
################
# WBSS
################
WBSS_df <- as.data.frame(HERAS_baseline$HERAS.WBSS$HERAS_strata@index[,ac(2017:2020),'abu','all'])
temp <- as.data.frame(HERAS_baseline$HERAS.WBSS$HERAS_all@index[,ac(2017:2020),'abu','all'])
WBSS_df <- rbind(WBSS_df,temp)
WBSS_df$age <- as.factor(WBSS_df$age)
WBSS_df <- WBSS_df[WBSS_df$data != 0,]
# WBSS N matyrity df
WBSS_mat_df <- as.data.frame(HERAS_baseline$HERAS.WBSS$HERAS_strata@index[,ac(2017:2020),'abu',c('IMM','MAT')])
temp <- as.data.frame(HERAS_baseline$HERAS.WBSS$HERAS_all@index[,ac(2017:2020),'abu',c('IMM','MAT')])
WBSS_mat_df <- rbind(WBSS_mat_df,temp)
WBSS_mat_df$age <- as.factor(WBSS_mat_df$age)
WBSS_mat_df <- NSAS_mat_df[WBSS_mat_df$data != 0,]
# WBSS N df
WBSS_N_bootstrap <- as.data.frame(HERAS_baseline_bootstrap$stk.WBSS@stock.n[,ac(2017:2020)])
WBSS_N_bootstrap$age <- as.factor(WBSS_N_bootstrap$age)
WBSS_N <- as.data.frame(HERAS_baseline$stk.WBSS@stock.n[,ac(2017:2020)])
WBSS_N$age <- as.factor(WBSS_N$age)
# WBSS SSB df
HERAS_WBSS <- HERAS_WBSS %>%
gather(key=age,value=abu,colnames(HERAS_WBSS)[2:dim(HERAS_WBSS)[2]])
wa_WBSS <- wa_WBSS %>%
gather(key=age,value=wa,colnames(wa_WBSS)[2:dim(wa_WBSS)[2]])
wa_WBSS_36 <- wa_WBSS[!is.na(match(wa_WBSS$age,unique(HERAS_WBSS$age))),]
mat_WBSS_36 <- mat_WBSS[!is.na(match(mat_WBSS$age,unique(HERAS_WBSS$age))),]
mat_WBSS_36 <- rep(mat_WBSS_36$mat,length(unique(HERAS_WBSS$year)))
HERAS_WBSS$SSB <- HERAS_WBSS$abu*wa_WBSS_36$wa*mat_WBSS_36
HERAS_WBSS_SSB <- HERAS_WBSS %>% group_by(year) %>% summarize(SSB=sum(SSB))
for(idxYear in ac(2017:2020)){
#wa_WBSS_year <- wa_WBSS[wa_WBSS$year == idxYear,]
#wa_WBSS_year_09 <- rbind(wa_WBSS_year,c('2017','9',wa_WBSS_year$wa[dim(wa_WBSS_year)[1]]))
#wa_WBSS_year_09$wa <- an(wa_WBSS_year_09$wa)
#mat_WBSS_09 <- rbind(mat_WBSS,c('9',mat_WBSS$mat[dim(mat_WBSS)[1]]))
#mat_WBSS_09$mat <- an(mat_WBSS_09$mat)
#temp1 <- sweep(HERAS_baseline_bootstrap$stk.WBSS@stock.n[,idxYear]*1e3,2,wa_WBSS_year_09$wa,"*")
#temp2 <- sweep(temp1,2,mat_WBSS_09$mat,"*")
#HERAS_baseline_bootstrap$stk.WBSS@stock[,idxYear] <- apply(temp2,6,'sum')
if(idxYear == '2020'){
wa_WBSS_year <- wa_WBSS[wa_WBSS$year == '2019',]$wa
}else{
wa_WBSS_year <- wa_WBSS[wa_WBSS$year == idxYear,]$wa
}
wa_WBSS_year_09 <- an(c(wa_WBSS_year,wa_WBSS_year[length(wa_WBSS_year)]))
mat_WBSS_09 <- c(mat_WBSS$mat,mat_WBSS$mat[length(mat_WBSS$mat)])
#HERAS_WBSS[HERAS_WBSS$year == '2017',]
temp1 <- sweep(HERAS_baseline_bootstrap$stk.WBSS@stock.n[,idxYear]*1e3,1:6,wa_WBSS_year_09,"*")
temp2 <- sweep(temp1,1:6,mat_WBSS_09,"*")
HERAS_baseline_bootstrap$stk.WBSS@stock[,idxYear] <- quantSums(temp2)
}
temp <- as.data.frame(HERAS_baseline_bootstrap$stk.WBSS@stock[,ac(2017:2020)])
temp$iter <- as.numeric(temp$iter)
SSB_WBSS_HERAS_boot <- temp %>%
group_by(year) %>%
summarise_at( vars(data),
list( Q25=~quantile(., probs = 0.25,na.rm=TRUE),
Q50=~quantile(., probs = 0.50,na.rm=TRUE),
Q75=~quantile(., probs = 0.75,na.rm=TRUE)))
# WBSS % difference df
B <- subset(HERAS_WBSS,year==2017 | year==2018 | year==2019)
B <- arrange(B,year)
for(idxYear in ac(2017:2019)){
wa_WBSS_year <- wa_WBSS[wa_WBSS$year == idxYear & (wa_WBSS$age == 3 | wa_WBSS$age == 4 | wa_WBSS$age == 5 | wa_WBSS$age == 6),]
mat_WBSS_36 <- mat_WBSS[(mat_WBSS$age == 3 | mat_WBSS$age == 4 | mat_WBSS$age == 5 | mat_WBSS$age == 6),]
temp1 <- sweep(HERAS_baseline$stk.WBSS@stock.n[ac(3:6),idxYear]*1e3,1:6,wa_WBSS_year$wa,'*')
temp2 <- sweep(temp1,1:6,mat_WBSS_36$mat,'*')
HERAS_baseline$stk.WBSS@stock[,idxYear] <- quantSums(temp2)
}
HERAS_baseline$stk.WBSS@stock[,ac(2017:2019)] <- (subset(HERAS_WBSS_SSB,year==2017 | year==2018 | year==2019)$SSB -
as.data.frame(HERAS_baseline$stk.WBSS@stock[,ac(2017:2019)])$data)/
subset(HERAS_WBSS_SSB,year==2017 | year==2018 | year==2019)$SSB
diff_WBSS <- as.data.frame(HERAS_baseline$stk.WBSS@stock.n[ac(3:6),ac(2017:2019)]*1e3)
diff_WBSS$data <- (B$abu-diff_WBSS$data)/B$abu
diff_WBSS <- rbind(diff_WBSS,as.data.frame(HERAS_baseline$stk.WBSS@stock[,ac(2017:2019)]))
diff_WBSS$age[diff_WBSS$age == 'all'] <- 'SSB (3-6)'
diff_WBSS$year <- as.factor(diff_WBSS$year)
# NSAS bio df
temp_FLR <- HERAS_baseline$stk.WBSS[ac(0:8),ac(2017:2019)]
for(idxYear in ac(2017:2019)){
wa_WBSS_year <- wa_WBSS[wa_WBSS$year == idxYear,]
mat_WBSS_temp <- mat_WBSS
temp_FLR@mat[,idxYear] <- mat_WBSS_temp$mat
temp_FLR@stock.wt[,idxYear] <- wa_WBSS_year$wa
}
A <- as.data.frame(HERAS_baseline$stk.WBSS@stock.wt[,ac(2017:2020)]*1e-3)
A$source <- 'HERAS'
A$type <- 'weight'
B <- as.data.frame(temp_FLR@stock.wt)
B$source <- 'assessment'
B$type <- 'weight'
C <- as.data.frame(HERAS_baseline$stk.WBSS@mat[,ac(2017:2020)])
C$source <- 'HERAS'
C$type <- 'maturity'
D <- as.data.frame(temp_FLR@mat)
D$source <- 'assessment'
D$type <- 'maturity'
WBSS_bio <- rbind(A,B,C,D)
################################################################
# plots
################################################################
################
# NSAS
################
# plot per strata
scaling_factor <- 1.5
png(file.path(figurePath,'NSAS_N_strata.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(NSAS_df,aes(x=year,y=data))+
theme_bw()+
geom_bar(aes(fill = age),stat="identity")+
ylab('abundance')+
theme(axis.text.x = element_text(angle=90))+
facet_wrap(~area,scale='free_y')
print(p)
dev.off()
# plot per strata MAT/IMM
scaling_factor <- 1.5
png(file.path(figurePath,'NSAS_N_mat_strata.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(NSAS_mat_df,aes(x=year,y=data))+
theme_bw()+
geom_bar(aes(fill = season),stat="identity")+
ylab('abundance')+
theme(axis.text.x = element_text(angle=90))+
labs(fill='')+
facet_wrap(~area,scale='free_y')
print(p)
dev.off()
# assessment plot
scaling_factor <- 1
png(file.path(figurePath,'NSAS_SSB_assessment.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot()+
theme_bw()+
geom_ribbon(data = subset(SSB_NSAS_assessment,year >1985),
aes(x=year,ymin = lbnd, ymax = ubnd,fill='assessment uncertainty'),alpha=0.5)+
geom_line(data = subset(SSB_NSAS_assessment,year >1985),
aes(x=year,y=value,color='assessment trajectory'))+
geom_line(data = subset(SSB_NSAS_HERAS,year>1985),
aes(x=year,y=data,color='HERAS index'))+
scale_color_grey(start=0,end=0.5)+
scale_fill_grey(start=0.4,end=1)+
ylab('SSB')+
labs(colour='',shape='',fill='')
#geom_point( data = subset(SSB_NSAS_HERAS_new,year>1985),
# aes(x=year,y=data,shape='HERAS newly derived'))
print(p)
dev.off()
# compare assessment with baseline
scaling_factor <- 0.8
png(file.path(figurePath,'NSAS_diff_baseline.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(diff_NSAS,aes(x=age,y=data,col=year))+
theme_bw()+
geom_point()+
ylab('% difference')
print(p)
dev.off()
# compare bio
scaling_factor <- 1
png(file.path(figurePath,'NSAS_bio_baseline.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(NSAS_bio,aes(x=age,y=data,col=source))+
theme_bw()+
geom_line()+
facet_grid(type~year,scale='free_y')+
ylab('weight age age/maturity at age')
print(p)
dev.off()
# comparion assessment with bootstraping
scaling_factor <- 0.8
png(file.path(figurePath,'NSAS_SSB_bootstrap.png'),
width = 14*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot()+
theme_bw()+
geom_ribbon(data = subset(SSB_NSAS_assessment,year >2005),
aes(x=year,ymin = lbnd, ymax = ubnd,fill='assessment uncertainty'),alpha=0.5)+
geom_ribbon(data = SSB_NSAS_HERAS_boot,
aes(x=year,ymin = Q25, ymax = Q75,fill='index uncertainty'),alpha=0.5)+
geom_line(data = subset(SSB_NSAS_assessment,year >2005),
aes(x=year,y=value,color='assessment'),size=1.2)+
geom_line(data = subset(SSB_NSAS_HERAS,year>2005),
aes(x=year,y=data,color='index'),size=1.2)+
geom_line(data = SSB_NSAS_HERAS_boot,
aes(x=year,y=Q50,color='median'),size=1.2)+
ylab('SSB')+
labs(colour='',shape='',fill='',shape='',linetype='')+
scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
scale_colour_manual(values=c("#999999", "#56B4E9", "#E69F00"))
print(p)
dev.off()
# comparion baseline and bootstrap runs
scaling_factor <- 1
png(file.path(figurePath,'NSAS_bootstrap_comp.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(NSAS_N_bootstrap,aes(x=age,y=data))+
theme_bw()+
geom_boxplot()+
geom_line(data=NSAS_N,aes(x = age,y=data,group=1,col='baseline'))+
labs(colour='',shape='',fill='',shape='',linetype='')+
ylab('abundance')+
facet_wrap(~year,scale='free_y')
print(p)
dev.off()
#temp <- as.data.frame(HERAS_baseline_bootstrap$HERAS.NSAS@stock[,ac(2017:2020)])
################
# WBSS
################
# plot per strata
scaling_factor <- 1
png(file.path(figurePath,'WBSS_N_strata.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(WBSS_df,aes(x=year,y=data))+
theme_bw()+
geom_bar(aes(fill = age),stat="identity")+
ylab('abundance')+
theme(axis.text.x = element_text(angle=90))+
facet_wrap(~area,scale='free_y')
print(p)
dev.off()
# plot per strata MAT/IMM
scaling_factor <- 1
png(file.path(figurePath,'WBSS_N_mat_strata.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(WBSS_mat_df,aes(x=year,y=data))+
theme_bw()+
geom_bar(aes(fill = season),stat="identity")+
ylab('abundance')+
theme(axis.text.x = element_text(angle=90))+
labs(fill='')+
facet_wrap(~area,scale='free_y')
print(p)
dev.off()
# assessment plot
scaling_factor <- 1
png(file.path(figurePath,'WBSS_SSB_assessment.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot()+
theme_bw()+
geom_ribbon(data = SSB_WBSS_assessment,
aes(x=Year,ymin = Low.1, ymax = High.1,fill='assessment uncertainty'),alpha=0.5)+
geom_line(data = SSB_WBSS_assessment,
aes(x=Year,y=SSB,color='assessment trajectory'))+
geom_line(data = HERAS_WBSS_SSB,
aes(x=year,y=SSB,color='HERAS index (age 3-6)'))+
scale_color_grey(start=0,end=0.5)+
scale_fill_grey(start=0.4,end=1)+
ylab('SSB')+
labs(colour='',shape='',fill='')
#geom_point( data = subset(SSB_NSAS_HERAS_new,year>1985),
# aes(x=year,y=data,shape='HERAS newly derived'))
print(p)
dev.off()
# compare assessment with baseline
scaling_factor <- 0.8
png(file.path(figurePath,'WBSS_diff_baseline.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(diff_WBSS,aes(x=age,y=data,col=year))+
theme_bw()+
geom_point()+
ylab('% difference')
print(p)
dev.off()
# compare bio
scaling_factor <- 1
png(file.path(figurePath,'WBSS_bio_baseline.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(WBSS_bio,aes(x=age,y=data,col=source))+
theme_bw()+
geom_line()+
facet_grid(type~year,scale='free_y')+
scale_x_continuous(breaks = c(0, 2, 4,6,8))+
ylab('weight age age/maturity at age')
print(p)
dev.off()
# comparion assessment with bootstraping
scaling_factor <- 0.8
png(file.path(figurePath,'WBSS_SSB_bootstrap.png'),
width = 14*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot()+
theme_bw()+
geom_ribbon(data = subset(SSB_WBSS_assessment,Year>2005),
aes(x=Year,ymin = Low.1, ymax = High.1,fill='assessment uncertainty'),alpha=0.5)+
geom_ribbon(data = subset(SSB_WBSS_HERAS_boot,year>2005),
aes(x=year,ymin = Q25, ymax = Q75,fill='index uncertainty (all)'),alpha=0.5)+
geom_line(data = subset(SSB_WBSS_assessment,Year>2005),
aes(x=Year,y=SSB,color='assessment'),size=1.2)+
geom_line(data = subset(HERAS_WBSS_SSB,year>2005),
aes(x=year,y=SSB,color='index (3-6)'),size=1.2)+
geom_line(data = subset(SSB_WBSS_HERAS_boot,year>2005),
aes(x=year,y=Q50,color='median (all)'),size=1.2)+
ylab('SSB')+
labs(colour='',shape='',fill='',shape='',linetype='')+
scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
scale_colour_manual(values=c("#999999", "#56B4E9", "#E69F00"))
print(p)
dev.off()
# comparion baseline and bootstrap runs
scaling_factor <- 1
png(file.path(figurePath,'WBSS_bootstrap_comp.png'),
width = 12*scaling_factor, height = 8*scaling_factor, units = "cm", res = 300, pointsize = 10)
p <- ggplot(WBSS_N_bootstrap,aes(x=age,y=data))+
theme_bw()+
geom_boxplot()+
geom_line(data=WBSS_N,aes(x = age,y=data,group=1,col='baseline'))+
labs(colour='',shape='',fill='',shape='',linetype='')+
ylab('abundance')+
facet_wrap(~year,scale='free_y')
print(p)
dev.off()
\ No newline at end of file
# script that derives the HERAS index directly from the StoX projects.
rm(list=ls())
library(Rstox)
library(FLFleet)
library(ggplotFL)
library(mgcv)
library(tidyverse)
library(ggplot2)
#Set up directories
#path <- 'J:/git/HERAS/'
path <- 'J:/git/heras_index_kbwot'
try(setwd(path),silent=TRUE)
mainPath <- file.path(".")
dataPath <- file.path(".","data")
stoxPath <- file.path(".","StoX")
outPath <- file.path(".","output")
functionPath <- file.path(".","functions")
resultsPath <- file.path(".","results")
figurePath <- file.path(".","figures",'kbwot','3_SA_error')
assessmentError <- file.path(".","assessment",'SA_error','assessment')
load(file.path(dataPath,'miscellaneous','NSH_HAWG2020_sf.Rdata'))
NSH_HAWG <- NSH.sam
load(file.path(assessmentError,'NSH_HAWG2020_sf_SA_error_2018.Rdata'))
NSH_error_2018 <- NSH.sam
load(file.path(assessmentError,'NSH_HAWG2020_sf_SA_error_2017_2018.Rdata'))
NSH_error_2017_2018 <- NSH.sam
load(file.path(resultsPath,'HERAS_baseline.RData'))
load(file.path(resultsPath,'HERAS_sensitivity_SA_error.RData'))
SSB_WBSS_assessment <- read.csv(file.path(dataPath,'WBSS','assessment_2020.csv'))
HERAS_WBSS <- read.csv(file.path(dataPath,'WBSS','WBSS_HERAS.csv'),check.names=FALSE)
wa_WBSS <- read.csv(file.path(dataPath,'WBSS','WBSS_wa.csv'),check.names=FALSE)
mat_WBSS <- read.csv(file.path(dataPath,'WBSS','WBSS_mProp.csv'),check.names=FALSE)
################################################################
# compute df
################################################################
################
# NSAS
################
# NSAS % difference df
A <- apply( HERAS_baseline$stk.NSAS@stock.n[,ac(2017:2020)]*
HERAS_baseline$stk.NSAS@stock.wt[,ac(2017:2020)]*
HERAS_baseline$stk.NSAS@mat[,ac(2017:2020)],2,sum)
B <- apply( HERAS_sensitivity_SA_error$stk.NSAS@stock.n[,ac(2017:2020)]*
HERAS_sensitivity_SA_error$stk.NSAS@stock.wt[,ac(2017:2020)]*
HERAS_sensitivity_SA_error$stk.NSAS@mat[,ac(2017:2020)],2,sum)
C <- as.data.frame((HERAS_baseline$HERAS.NSAS$HERAS_strata@index[,ac(2017:2020),'abu','all']-
HERAS_sensitivity_SA_error$HERAS.NSAS$HERAS_strata@index[,ac(2017:2020),'abu','all'])/
HERAS_baseline$HERAS.NSAS$HERAS_strata@index[,ac(2017:2020),'abu','all'])