In einarhjorleifsson/namr: What the package does (one line)

#devtools::install_github("einarhjorleifsson/namr", build_vignettes = TRUE)
library(namr)
library(ggplot2)
library(tidyr)
library(dplyr)

repfile <- system.file("ass/hm/2015namhors.rep", package = "namr")
rbx <- read_rbx_hmac(repfile)
datfile <- system.file("ass/hm/2015namhorsdata.dat", package = "namr")
ibx <- read_ibx_hmac(datfile)
rbx <- read_hmac(repfile)

rbyf <- rbx$rbyf
rbyf %>% 
  select(year, fleet, landings) %>% 
  ggplot() +
  theme_bw() +
  geom_bar(aes(year, landings, fill = fleet), stat = "identity") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = seq(1960, 2020, by=10)) +
  labs(x = NULL, y = NULL, title = "Landings by different fleet and nations")

rbyaf <- rbx$rbyaf
rbyaf %>% 
  filter(!is.na(oC)) %>% 
  ggplot(aes(year, age)) +
  theme_bw() +
  geom_point() +
  facet_grid(fleet ~ .) +
  ylim(0,7) +
  scale_x_continuous(breaks = seq(1975, 2015, 5)) +
  labs(x = NULL, y = "Age", title = "The catch at age observation space")

The above figure shows the catch at age that is input into the model. There are two features that raises questions:

The maximum age is 5 years old for any of the fleet. Yet the model has age groups up to and including age 7. It also report fishing mortality in those older age groups (i.e. age 6 and 7).
Minor: Within the age groups 5 and below there are in some years zero values (interpreted in my R-scrpt as being NA) in some years. Wander if that is an "error".

rbyf %>%
  filter(fleet %in% c("Bulgaria", "Poland",
                      "Rumania", "USSR")) %>% 
  group_by(fleet) %>% 
  mutate(cpue = oU/mean(oU[year %in% c(1980:1986)])) %>% 
  filter(year %in% 1972:1987) %>% 
  ggplot(aes(year, cpue, colour = fleet)) +
  geom_hline(yintercept = 1) +
  theme_bw() +
  geom_point() +
  geom_smooth() +
  scale_colour_brewer(palette = "Set1") +
  labs(x = NULL, y = NULL, colour = "Fleet", title = "Catch per unit effort")

rbyf %>% 
  select(year, fleet, oU, pU) %>% 
  filter(fleet != "Pelagic",
         year <= 2015) %>% 
  ggplot(aes(year)) +
  theme_bw() +
  geom_point(aes(y = oU)) +
  geom_line(aes(y = pU)) +
  facet_wrap(~ fleet, scale = "free_y") +
  labs(x = NULL, y = NULL, title = "Observed and predicted catch per unit effort") +
  expand_limits(y = 0)

rbyaf <- rbx$rbyaf
rbyaf %>% 
  filter(fleet != "Rumania",
         year %in% c(1970:2015),
         age < 6) %>% 
  ggplot() +
  theme_bw() +
  geom_hline(yintercept = 0, col = "grey") +
  geom_point(aes(year, rC)) +
  facet_grid(age ~ fleet) +
  labs(x = NULL, y = NULL, title = "Catch at age residuals")

Questions/issues:

Why no residuals for age 6 and 7???

rbya <- rbx$rbya
rbya <- left_join(rbya,ibx$ibya)
names(rbya) <- c("year","age","n","f","sW","cW")
rbya <- left_join(rbya,ibx$iba)
rbya$bio <- rbya$n * rbya$sW
rbya$m <- 0.45
rbya$pC <- rbya$f/(rbya$f+rbya$m) * (1 - exp(-(rbya$f + rbya$m))) * rbya$n
rbya$pY <- rbya$pC * rbya$cW
ggplot(rbya,aes(year,bio,fill = factor(age))) + 
  geom_bar(stat = "identity") +
  scale_fill_brewer(palette = "Set1") +
  labs(x = NULL,y = NULL, title="Horse mackerel: Total biomass") +
  theme(legend.position=c(0.7,0.7))
ggplot(rbya,aes(year, pY, fill = factor(age))) +
  geom_bar(stat = "identity") +
  scale_fill_brewer(palette = "Set1") +
  labs(x =NULL, y = NULL, title="Horse mackerel: Catch")+theme(legend.position=c(0.8,0.74))

rbya %>% 
  filter(year <= 2013) %>% 
  group_by(year) %>% 
  summarise(pY = sum(pY)) %>% 
  ggplot() + 
  geom_line(aes(year,pY)) +
  geom_point(data=rby[rby$year <= 2013,],aes(year,yield)) +
  labs(x="",y="",title="Obverved (points) vs. predicted (line) yield")

rbya %>% 
  filter(year <= 2013) %>% 
  ggplot(aes(year, f, group = age)) + 
  geom_line() +  
  geom_text(aes(label=age),size=3) +
  labs(x="",y="",title="Fishing mortality at age")

rby %>% 
  filter(year <= 2014) %>% 
  ggplot(aes(x=year)) + 
  geom_line(aes(y = pU), col = "red", lwd = 3)+
  geom_point(aes(y = oU), col = "blue") +
  labs(x = NULL, y = NULL, title = "Observed and predicted survey indices")

rbx$rbyf %>% 
  filter(year < 2014, fleet != "Pelagic") %>% 
  ggplot(aes(x = year)) + 
  geom_line(aes(y = pCPUE), col = "red", lwd = 3) +
  geom_point(aes(y = oCPUE), col = "blue") +
  facet_wrap(~fleet, scale = "free_y") +
  scale_y_continuous(expand = c(0,0)) +
  expand_limits(y = 0) +
  labs(x = NULL, y = NULL, title = "Observed and predicted commercial indices")

rbx$rbyf %>% 
  filter(year < 2014, fleet == "Midwater") %>% 
  ggplot(aes(x = year)) +
  geom_line(aes(y = pCPUE), col = "red", lwd = 3) +
  geom_point(aes(y = oCPUE), col = "blue") +
  scale_y_continuous(expand = c(0,0)) +
  expand_limits(y = 0) +
  labs(x = NULL, y = NULL, title = "Observed and predicted midwater fleet indices")

rba  <- rbx$rba
rba1 <-  melt(rba, id.vars = "age", variable.name = "fleet")
rba1$age <- as.numeric(rba1$age)
ggplot(rba1, aes(age, value, col=fleet)) + 
  geom_line(size=2) +
  labs(x = NULL, y = NULL, title = "Estimated selection pattern") +
  scale_y_continuous(expand = c(0,0)) + 
  theme(legend.position = c(0.8,0.42)) + 
  scale_color_brewer(palette = "Set1")

rbya %>% 
  filter(year %in% seq(1961,2011,by=10)) %>% 
  ggplot(aes(age, n, col = factor(year))) +
  geom_line(size = 2) +
  scale_y_continuous(expand = c(0,0)) + 
  theme(legend.position = c(0.8,0.42)) +
  scale_color_brewer(palette="Set1") +
  labs(x = "Age", y = NULL, title = "Population numbers in selected years")

rbya %>% 
  mutate(bio = n * cW) %>% 
  filter(year %in% seq(1961,2011,by=10)) %>% 
  ggplot(aes(age, bio, col = factor(year))) + 
  geom_line(size = 2) + 
  scale_y_continuous(expand = c(0,0)) + 
  theme(legend.position = c(0.8,0.42)) + 
  scale_color_brewer(palette="Set1") +
  labs(x = "Age", y = NULL, title = "Population biomass in selected years")

dat <- rby[,c("year",paste("yield",1:6,sep=""))]
names(dat)[2:7] <- FLEETS[1:6]
dat1 <- melt(dat,"year", variable = "fleet", value.name = "yield")
dat1$year <- as.numeric(dat1$year)
ggplot(dat1, aes(year, yield, fill=factor(fleet))) + 
  geom_bar(stat = "identity") + 
  scale_fill_brewer(palette="Set1") +
  scale_y_continuous(expand=c(0,0)) + 
  theme(legend.position=c(0.8,0.78)) +
  labs(x = NULL, y = NULL, title = "Catch by fleet")

EINAR: Code below not run, needs to be rewritten

library(plyr)
rbya <- rbx$rbya
rbya <- ddply(rbya,.(year),transform,oC=n/sum(n))
rbya$fleet <- "stock"
rbya <- rbya[,c("year","age","fleet","oC")]
rbyaf <- rbx$rbyaf[!is.na(rbx$rbyaf$oC),c("year","age","fleet","oC")]
rbya <- rbya[rbya$year %in% rbyaf$year,]
rbyaf <- rbind(rbyaf,rbya)
rbyaf$age <- as.numeric(rbyaf$age)


dat  <- rbyaf[rbyaf$fleet %in% c("Poland","USSR"),]
dat$age <- as.numeric(dat$age)
dat <- rbyaf[rbyaf$year %in% dat$year,]
dat <- dat[dat$fleet %in% c("Poland","USSR","stock"),]
ggplot(dat,aes(age,oC,col=factor(fleet)))+geom_line(size=2)+
facet_wrap(~year)+scale_y_continuous(expand=c(0,0))+expand_limits(y=0)+scale_color_brewer(palette="Set1")

dat  <- rbx$rbyaf[ !is.na(rbx$rbyaf$oC) & rbx$rbyaf < 1995,]
ggplot(dat,aes(age,oC,col=factor(fleet)))+geom_line(size=2)+
facet_wrap(~year) +
  scale_y_continuous(expand=c(0,0)) +
  expand_limits(y=0) + 
  scale_colour_brewer(palette="Set1")

dust bin

rby <- rbx$rby
i <- rby$year <= 2015
p <- ggplot(rby[i,],aes(x=year)) + labs(x="",y="")
p + 
  geom_line(aes(y=ssb),col="red",lwd=3) +
  geom_line(aes(y=tBio),col="blue",lwd=2) +
  geom_line(aes(y=landings),col="green",lwd=2)