setwd("~/ass/2013/WKMSYREF/rmd/cod-iceg/")
SIM <- FALSE
rbx2012 <- read.adcam("~/ass/2012/01/05ass/Adcam", "smx", assYear = 2012, calcSurBio = T)
rbya <- rbx2012$rbya
names(rbya) <- c("year", "age", "catch.n", "catch.wt", "stock.wt", "stock1.wt",
"mat", "stock.n", "z", "harvest", "m", "Chat", "rC", "Uobs", "Uhat", "Ures",
"U2obs", "U2hat", "U2res", "run", "model", "assYear")
rbya$catch.wt <- ifelse(rbya$catch.wt == -1, 0, rbya$catch.wt)
rbya$stock.wt <- ifelse(rbya$stock.wt == -1, 0, rbya$stock.wt)
# rbya$catch.n <- ifelse(rbya$catch.n== -1,0,rbya$catch.n)
rbya$harvest <- ifelse(rbya$harvest == -1, 0, rbya$harvest)
rbya$m <- ifelse(rbya$age %in% c(1, 2), 1e-04, rbya$m)
rbya$mat <- ifelse(rbya$age %in% c(1, 2) | is.na(rbya$mat), 0, rbya$mat)
rbya <- rbya[rbya$year %in% c(1955:2012) & rbya$age %in% c(3:15), ]
rbyaStored <- rbya
rby <- rbx2012$rby
rby <- rby[rby$year %in% c(1955:2012), ]
# convert rbyaa to FLQuants
require(FLCore)
## Loading required package: FLCore
## Loading required package: grid
## Loading required package: lattice
##
## FLCore 2.4 "The Duke of Prawns"
## ---------------------------------------------------------
## * Note that FLR packages are under constant development,
## please report bugs to flr-team@flr-project.org.
## * New documentation can be found in vignetes,
## run vignette(package="FLCore") to get a list.
## * For more information go to http:\\flr-project.org or
## subscribe the mailing list flr-list@flr-project.org.
## ---------------------------------------------------------
## Attaching package: 'FLCore'
## The following object(s) are masked from
## 'file:/home/einarhj/stasi/MfS.RData':
##
## cv
## The following object(s) are masked from 'package:plyr':
##
## desc
## The following object(s) are masked from 'package:ggplot2':
##
## %+%
## The following object(s) are masked from 'package:base':
##
## cbind, rbind
x <- FLQuants()
for (i in 3:17) {
df <- rbya[, c(1, 2, i)]
nome <- names(rbya[i])
names(df)[3] <- "data"
x[[nome]] <- as.FLQuant(df)
}
iCodLong <- FLStock(name = "iCod", desc = "from scratch", catch.n = x$catch.n,
catch.wt = x$catch.wt, landings.n = x$catch.n, landings.wt = x$catch.wt,
stock.n = x$stock.n, stock.wt = x$stock.wt, m = x$m, mat = x$mat, harvest = x$harvest)
catch(iCodLong) <- rby$oY
discards.n(iCodLong) <- 0
discards.wt(iCodLong) <- 0
discards(iCodLong) <- 0
landings(iCodLong) <- rby$oY
harvest.spwn(iCodLong) <- c(0.085, 0.18, 0.248, 0.296, 0.382, 0.437, rep(0.477,
6))
m.spwn(iCodLong) <- rep(0.25, 12)
units(iCodLong)[1:17] <- as.list(c(rep(c("tonnes", "thousands", "kg"), 4), "NA",
"NA", "f", "NA", "NA"))
range(iCodLong) <- c(3, 14, 14, 1955, 2012, 5, 10)
stockFLR <- iCodLong
save(stockFLR, file = "stockFLR.RData")
plot(iCodLong)
## Warning: invalid factor level, NAs generated
rby2 <- data.frame(year = as.data.frame(catch(iCodLong))$year, n3 = as.data.frame(rec(iCodLong))$data,
oY = as.data.frame(catch(iCodLong))$data, ssb = as.data.frame(ssb(iCodLong))$data,
f = as.data.frame(fbar(iCodLong))$data)
d <- melt(rby[, c("year", "n3", "oY", "ssb", "f")], id.vars = "year")
d2 <- melt(rby2[, c("year", "n3", "oY", "ssb", "f")], id.vars = "year")
p <- ggplot(d, aes(year, value)) + geom_line(lwd = 1, col = "blue") + facet_wrap(~variable,
scale = "free_y")
p + geom_line(data = d2, aes(year, value), col = "red", alpha = 0.75)
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_path).
# rm(list=ls())
require(MASS)
## Loading required package: MASS
require(FLCore)
source("../../software/SimulateRefPts2.R")
source("../../software/plotting.R")
source("../../software/myMH.R")
attach("~/stasi/icesSS.RData")
## A function only available for R 2.15.1
paste0 <- function(...) paste(..., sep = "")
stock <- "cod-iceg"
runid = stock
odir <- gsub(" ", "-", runid)
i <- icesSS$assYear %in% 2012 & icesSS$stock %in% stock
dat <- icesSS[i, c("ssb", "r", "year")]
names(dat)[2] <- "rec"
dat$ssb <- dat$ssb/1000
dat$rec <- dat$rec/1000
# res <- dir.create(odir) if (res == FALSE) cat('Make sure you want to
# write over existing files.\n Better to #empty folder before
# proceeding.\n')
save(odir, runid, dat, file = paste("data.rData", sep = "/"))
fit <- fitModels(dat, runid)
## Posterior summary estimates for each model:
##
## mean (SD)
## a:
## 1.286 ( 0.01347 ) acc rate (try for 0.4) : 0.133
##
## b:
## 0.00228 ( 6.115e-08 ) acc rate (try for 0.4) : 0.133
##
## cv:
## 0.381 ( 0.001367 ) acc rate (try for 0.4) : 0.133
##
## Posterior summary estimates for each model:
##
## mean (SD)
## a:
## 3.262e+21 ( 6.51e+43 ) acc rate (try for 0.4) : 0.164
##
## b:
## 1.987e+19 ( 2.471e+39 ) acc rate (try for 0.4) : 0.164
##
## cv:
## 0.3619 ( 0.00132 ) acc rate (try for 0.4) : 0.164
##
## Posterior summary estimates for each model:
##
## mean (SD)
## a:
## 0.006094 ( 8.903e-08 ) acc rate (try for 0.4) : 0.184
##
## b:
## 1.201e+09 ( 2.812e+18 ) acc rate (try for 0.4) : 0.184
##
## cv:
## 0.3607 ( 0.001186 ) acc rate (try for 0.4) : 0.184
save(fit, file = paste("fit.rData", sep = "/"))
# attach('stockFLR.RData')
stk <- stockFLR
if (SIM) {
sim <- EqSim(fit, stk, Fscan = seq(0, 1, length = 40))
save(sim, file = paste("sim.rData", sep = "/"))
}
if (!SIM) attach("sim.rData")
# png(paste0(odir, '/data.png'), 800, 800)
with(dat, {
ssb.scl <- 10
rec.scl <- 100
plot(ssb/ssb.scl, rec/rec.scl, type = "n", las = 1, xlab = sprintf("SSB (%is)",
ssb.scl), ylab = sprintf("Recruits (%is)", rec.scl), main = paste(runid,
"stock-recruit data"))
lines(ssb/ssb.scl, rec/rec.scl, type = "b", cex = 3)
text(x = ssb/ssb.scl, y = rec/rec.scl, label = substring(year, 3, 4), cex = 0.7)
})
# dev.off()
# plot the fits ----------------------------
# png(paste0(odir, '/HS-fit.png'), 800, 800)
plot(fit$HS)
# dev.off()
# png(paste0(odir, '/RK-fit.png'), 800, 800)
plot(fit$RK)
# dev.off()
# png(paste0(odir, '/BH-fit.png'), 800, 800)
plot(fit$BH)
# dev.off()
par(mfrow = c(2, 2))
LLmodplot(fit)
SRplot(fit)
# dev.off()
summary(fit$HS)
## a b cv llik
## Min. :1.54e+07 Min. :1.60e-08 Min. :0.280 Min. :-29.8
## 1st Qu.:8.69e+07 1st Qu.:9.00e-08 1st Qu.:0.336 1st Qu.:-23.2
## Median :4.29e+08 Median :3.74e-07 Median :0.359 Median :-22.6
## Mean :1.20e+09 Mean :1.35e-06 Mean :0.361 Mean :-22.9
## 3rd Qu.:1.80e+09 3rd Qu.:1.90e-06 3rd Qu.:0.382 3rd Qu.:-22.1
## Max. :9.43e+09 Max. :1.12e-05 Max. :0.505 Max. :-21.8
summary(fit$RK)
## a b cv llik
## Min. :0.98 Min. :0.00149 Min. :0.289 Min. :-34.9
## 1st Qu.:1.21 1st Qu.:0.00210 1st Qu.:0.354 1st Qu.:-26.3
## Median :1.27 Median :0.00227 Median :0.377 Median :-25.5
## Mean :1.29 Mean :0.00228 Mean :0.381 Mean :-25.8
## 3rd Qu.:1.36 3rd Qu.:0.00245 3rd Qu.:0.404 3rd Qu.:-24.9
## Max. :1.76 Max. :0.00319 Max. :0.562 Max. :-24.2
summary(fit$BH)
## a b cv llik
## Min. :136 Min. :0.00e+00 Min. :0.262 Min. :-30.1
## 1st Qu.:160 1st Qu.:0.00e+00 1st Qu.:0.337 1st Qu.:-23.3
## Median :165 Median :0.00e+00 Median :0.359 Median :-22.6
## Mean :165 Mean :6.95e-15 Mean :0.362 Mean :-22.9
## 3rd Qu.:171 3rd Qu.:2.50e-16 3rd Qu.:0.382 3rd Qu.:-22.1
## Max. :201 Max. :1.63e-13 Max. :0.497 Max. :-21.8
table(fit$SR[1, ])/nrow(fit$SR)
## , , b = 0.00160192729167472, sigma = 0.399245693558867, llik = -27.6599059950524
##
## a
## mod 1.07020752855216
## RK 2e-04
# png(paste0(odir, '/refpts.png'), 800, 800) rfpts <- Eqplot(sim, stk,
# fit, Blim = 0.4 * max(fit $ dat $ ssb))
rfpts <- Eqplot(sim, stk, fit, Blim = 120)
# dev.off()
write.table(as.data.frame(t(c(round(rfpts[1:2]), round(rfpts, 2)[-c(1:2)]))),
file = paste0(odir, "/refpts.csv"), sep = ",", row.names = FALSE)
## Warning: cannot open file 'cod-iceg/refpts.csv': No such file or directory
## Error: cannot open the connection