Nothing
R/FLXSA.R
In FLXSA: eXtended Survivor Analysis for FLR
Defines functions
validFLXSA.control
validFLXSA
FLXSA.control
is.FLXSA
is.FLXSA.control
Documented in
FLXSA.control
is.FLXSA
is.FLXSA.control
# FLXSA.R -
# FLXSA/R/FLXSA.R
# Copyright 2003-2011 FLR Team. Distributed under the GPL 2 or later
# Maintainer: Iago Mosqueira, JRC
# $Id: $
# {{{ FLXSA.control class
validFLXSA.control <- function(object){
if (object@tol <= 0)
return("value of tol must be > 0")
if (object@maxit <= 0)
return("value of maxit must be > 0")
if (object@min.nse < 0)
return("value of min.nse must be > 0")
if (object@fse < 0)
return("value of fse must be > 0")
if (object@rage < -1)
return("value of rage must be >= -1")
if (object@qage < 0)
return("value of qage must be >= 0")
if (object@shk.yrs < 0)
return("value of shk.yrs must be >= 0")
if (object@shk.ages < 0)
return("value of shk.ages must be >= 0")
if (object@window <= 5)
return("value of window must be >= 5")
# Everything is fine
return(TRUE)
}
setClass("FLXSA.control",
representation(
tol ="numeric",
maxit ="integer",
min.nse ="numeric",
fse ="numeric",
rage ="integer",
qage ="integer",
shk.n ="logical",
shk.f ="logical",
shk.yrs ="integer",
shk.ages ="integer",
window ="integer",
tsrange ="integer",
tspower ="integer",
vpa ="logical"),
prototype=prototype(
tol =as.double(10e-10),
maxit =as.integer(30),
min.nse =as.double(0.3),
fse =as.double(0.5),
rage =as.integer(-1),
qage =as.integer(10),
shk.n =TRUE,
shk.f =TRUE,
shk.yrs =as.integer(5),
shk.ages =as.integer(5),
window =as.integer(100),
tsrange =as.integer(20),
tspower =as.integer(3),
vpa =FALSE),
validity=validFLXSA.control
)
# }}}
# FLXSA class {{{
validFLXSA <- function(object){
# All FLQuant objects must have same dimensions
return(TRUE)
Dim <- dim(object@stock.n)
if (!all(dim(object@harvest) == Dim))
return("n and f arrays must have same dimensions")
# Everything is fine
return(TRUE)
}
setClass("FLXSA",
contains='FLAssess',
representation(
survivors="FLQuant",
se.int ="FLQuant",
se.ext ="FLQuant",
n.fshk ="FLQuant",
n.nshk ="FLQuant",
var.fshk ="FLQuant",
var.nshk ="FLQuant",
q.hat ="FLQuants",
q2.hat ="FLQuants",
diagnostics="data.frame",
control ="FLXSA.control"),
prototype=prototype(
survivors=FLQuant(quant='age'),
se.int =FLQuant(quant='age'),
se.ext =FLQuant(quant='age'),
n.fshk =FLQuant(quant='age'),
n.nshk =FLQuant(quant='age'),
var.fshk =FLQuant(quant='age'),
var.nshk =FLQuant(quant='age'),
catch.n =FLQuant(quant='age'),
stock.n =FLQuant(quant='age'),
harvest =FLQuant(quant='age'),
q.hat =FLQuants(),
q2.hat =FLQuants(),
diagnostics=new("data.frame"),
control =new("FLXSA.control")),
validity=validFLXSA
) # }}}
# FLXSA() {{{
setGeneric("FLXSA", function(stock, indices, ...)
standardGeneric("FLXSA"))
setMethod("FLXSA", signature(stock="FLStock", indices="FLIndex"),
function(stock, indices, control=FLXSA.control(), diag.flag=TRUE) {
FLXSA(stock=stock, indices=FLIndices(indices), control=control, diag.flag=diag.flag)
}
)
setMethod("FLXSA", signature(stock="FLStock", indices="FLIndices"),
function(stock, indices, control=FLXSA.control(), diag.flag=TRUE) {
Call <- match.call()
# check FLIndices input
for (i in 1:length(indices)) {
# startf & endf present in @range
if (is.na(indices[[i]]@range["startf"]) || is.na(indices[[i]]@range["endf"]))
stop(paste("Must supply startf & endf for range in FLIndex",i))
# @range has all elements
if (!all(c("min","max","plusgroup","minyear","maxyear","startf","endf") %in%
names(indices[[i]]@range)))
stop("Range must have names 'min','max','plusgroup','minyear','maxyear',
'startf','endf'")
# adjust ranges to dims()
indices[[i]]@range["min"] <- max(indices[[i]]@range["min"], dims(indices[[i]])$min,
stock@range["min"])
indices[[i]]@range["max"] <- min(indices[[i]]@range["max"], dims(indices[[i]])$max,
stock@range["max"])
indices[[i]]@range["minyear"] <- max(indices[[i]]@range["minyear"],
dims(indices[[i]])$minyear, stock@range["minyear"])
indices[[i]]@range["maxyear"] <- min(indices[[i]]@range["maxyear"],
dims(indices[[i]])$maxyear, stock@range["maxyear"])
# trim according to range
age <- indices[[i]]@range["min"]:indices[[i]]@range["max"]
year<- indices[[i]]@range["minyear"]:indices[[i]]@range["maxyear"]
indices[[i]] <- trim(indices[[i]], age=age, year=year)
}
# Double check validity
if (!validObject(stock))
stop("stock is not valid!")
if (!validObject(indices))
stop("FLIndices is not valid!")
if (!validObject(control))
stop("control is not valid!")
# adjust range for minage
if ("minage" %in% names(stock@range))
minage <- stock@range["minage"]
else if ("min" %in% names(stock@range))
minage <- stock@range["min"]
else if ("minquant" %in% names(stock@range))
minage <- stock@range["minquant"]
else
stop("'minage' not found in range")
# adjust range for maxage
if ("maxage" %in% names(stock@range))
maxage <- stock@range["maxage"]
else if ("max" %in% names(stock@range))
maxage <- stock@range["max"]
else if ("maxquant" %in% names(stock@range))
maxage <- stock@range["maxquant"]
else
stop("'maxage' not found in range")
# adjust plsugroup
if ("plusgroup" %in% names(stock@range))
stock@range["plusgroup"] <- maxage
if (maxage<minage | stock@range["maxyear"] < stock@range["minyear"])
stop("Error in range")
if (is.na(stock@range["plusgroup"]))
stop("Plus Group must be specified")
# trim stock
stock <- trim(stock, year=stock@range["minyear"]:stock@range["maxyear"])
if (!is.na(stock@range["plusgroup"]) &
stock@range["plusgroup"] < dims(stock@catch.n)$max)
stock <- setPlusGroup(stock, stock@range["plusgroup"])
stock@m <- stock@m[as.character(minage:maxage),,,,]
# check catch.n is available
if (all(is.na(stock@catch.n)))
stop("catch.n is not available")
stock@stock.n[] <- new('FLQuant')
stock@harvest[] <- new('FLQuant')
# fqs
fqs <- function(assess) {
assess@index <- new("FLQuants", lapply(assess@index,FLQuant))
assess@index.hat <- new("FLQuants", lapply(assess@index.hat,FLQuant))
assess@index.var <- new("FLQuants", lapply(assess@index.var,FLQuant))
assess@index.res <- new("FLQuants", lapply(assess@index.res,FLQuant))
assess@q.hat <- new("FLQuants", lapply(assess@q.hat,FLQuant))
assess@q2.hat <- new("FLQuants", lapply(assess@q2.hat,FLQuant))
if (validObject(assess))
return(assess)
else
stop("not valid")
}
#
iters.stock <-max(unlist(qapply(stock, function(x) dims(x)$iter)))
iters.indices<-max(unlist(lapply(indices@.Data,
function(x) max(unlist(qapply(x, function(x2) dims(x2)$iter))))))
if ((iters.stock>1 && iters.indices>1) && missing(diag.flag))
diag.flag <- FALSE
if ((iters.stock>1 && iters.indices>1) && diag.flag)
return("Multiple iters only allowed if diag.flag=FALSE")
if(!diag.flag) {
res<-.Call("FLXSA", iter(stock, 1), lapply(indices, iter, 1), control, FALSE)
iters <- max(iters.stock,iters.indices)
if (iters>1) {
res@stock.n<-propagate(FLQuant(res@stock.n@.Data),iters)
res@harvest<-propagate(FLQuant(res@harvest@.Data),iters)
for (i in as.character(2:iters)) {
res. <- .Call("FLXSA", iter(stock,i), lapply(indices, iter,i), control, FALSE)
iter(res@stock.n,i)<-FLQuant(res.@stock.n@.Data)
iter(res@harvest,i)<-FLQuant(res.@harvest@.Data)
}
}
res@harvest@units <- "f"
return(res)
}
res <-.Call("FLXSA", stock, indices, control, diag.flag)
res <-fqs(.Call("FLXSA", stock, indices, control, diag.flag))
if (class(res) != "FLXSA")
return(res)
res@call <- as.character(Call)
# put wts amd nhats into a data frame
df <- as.data.frame(res@wts)
df1 <- (df[4])
df1[df1 >= 1, 1] <- paste("index", df1[df1 >= 1, 1])
df1[df1 == -1, 1] <- "fshk"
df1[df1 == -2, 1] <- "nshk"
df <- cbind(df[-4], df1)
names(df) <- c("w", "nhat", "yrcls", "age", "year", "source")
for(i in 1:length(indices)) {
v <- paste("index", i)
if (length(slot(indices[[i]],"name"))>0)
df$source[df$source==v] <- slot(indices[[i]],"name")
}
wts.df <-df[,c(4,5,1,6)]
names(wts.df)[3] <-"data"
nhat.df<-df[,c(4,5,2,6)]
names(nhat.df)[3]<-"data"
index.names<-sort(unique(df[,"source"]))
index.names<-index.names[substr(index.names,1,5)=="index"]
fill.flq<-function(obj) {
dms <-dims(obj)
dmns<-dimnames(obj)
dmns[[1]]<-dms[[2]]:dms[[3]]
dmns[[2]]<-dms[[5]]:dms[[6]]
res<-as.FLQuant(NA,dimnames=dmns)
res[dimnames(obj)[[1]],dimnames(obj)[[2]],,,]<-obj[dimnames(obj)[[1]],
dimnames(obj)[[2]],,,]
return(res)
}
res2 <- new("FLXSA")
res2@index.var<-new('FLQuants')
res2@index.hat<-new('FLQuants')
res2@index <-new('FLQuants')
j=0
for (i in index.names) {
j=j+1
res2@index.var[[j]]<-1.0/fill.flq(as.FLQuant(wts.df[wts.df[, "source"]==i,-4]))
res2@index.hat[[j]]<- fill.flq(as.FLQuant(nhat.df[nhat.df[,"source"]==i,-4]))
res2@index.var[[j]]<-1.0/as.FLQuant(wts.df[wts.df[, "source"]==i,-4])
res2@index.hat[[j]]<-as.FLQuant(nhat.df[nhat.df[,"source"]==i,-4])
dmns <- dimnames(res2@index.hat[[j]])
index <- trim(indices[[j]]@index,age=dmns$age,year=dmns$year)
res2@index[ is.na(index)]<-NA
res2@index.hat[is.na(index)]<-NA
res2@index.var[is.na(index)]<-NA
res2@index.res[is.na(index)]<-NA
}
# F shrinkage
fshk <- df[df[,"source"]=="fshk",]
if (length(fshk[,1])>0) {
y.range <- range(fshk[,"year"])
a.range<-range(fshk[,"age"])
max.yr <-fshk[fshk[,"age"] ==a.range[2],]
max.age<-fshk[fshk[,"year"]==y.range[2],]
res2@n.fshk <-as.FLQuant(NA,dimnames=list(age=a.range[1]:a.range[2],
year=y.range[1]:y.range[2]))
res2@n.fshk[as.character(max.age[,"age"]),as.character(y.range[2])]<-max.age[,
"nhat"]
res2@n.fshk[as.character(a.range[2]),as.character(max.yr[,"year"])]<-max.yr[,"nhat"]
res2@var.nshk <-as.FLQuant(NA,dimnames=list(age=a.range[1]:a.range[2],
year=y.range[1]:y.range[2]))
res2@var.nshk[as.character(max.age[,"age"]),as.character(y.range[2])]<-
1/max.age[,"w"]
res2@var.nshk[ as.character(a.range[2]),as.character(max.yr[,"year"])]<-
1/max.yr[,"w"]
}
# N shrinkage
nshk <-df[df[,"source"]=="nshk",]
if (length(nshk[,1])>0) {
y.range<-range(nshk[,"year"])
a.range<-range(nshk[,"age"])
max.yr <-nshk[nshk[,"age"] ==a.range[2],]
max.age<-nshk[nshk[,"year"]==y.range[2],]
res2@n.nshk <-as.FLQuant(NA,dimnames=list(age=a.range[1]:a.range[2],
year=y.range[1]:y.range[2]))
res2@n.nshk[,as.character(y.range[2])] <- max.age[,"nhat"]
res2@n.nshk[ as.character(a.range[2])] <- max.yr[,"nhat"]
res2@var.nshk <- as.FLQuant(NA,dimnames=list(age=a.range[1]:a.range[2],
year=y.range[1]:y.range[2]))
res2@var.nshk[,as.character(y.range[2])]<-1/max.age[,"w"]
res2@var.nshk[ as.character(a.range[2])]<-1/max.yr[,"w"]
}
res2@diagnostics<- df
res2@index.hat <- res@index.hat
res2@stock.n <- FLQuant(res@stock.n@.Data)
res2@harvest <- FLQuant(res@harvest@.Data)
res2@survivors <- FLQuant(res@survivors@.Data)
res2@se.int <- FLQuant(res@se.int@.Data)
res2@se.ext <- FLQuant(res@se.ext@.Data)
res2@n.fshk <- FLQuant(res@n.fshk@.Data)
res2@n.nshk <- FLQuant(res@n.nshk@.Data)
res2@var.fshk <- FLQuant(res@var.fshk@.Data)
res2@var.nshk <- FLQuant(res@var.nshk@.Data)
res2@harvest@units <- "f"
res2@q.hat <- res@q.hat
res2@q2.hat <- res@q2.hat
res2@index.res <- res@index.res
res2@index <- res@index
res2@index.hat <-FLQuants(res@index.hat)
res2@index.res <-FLQuants(res@index.res)
res2@index.var <-FLQuants(res@index.var)
res2@q.hat <-FLQuants(res@q.hat)
res2@q2.hat <-FLQuants(res@q2.hat)
for (i in 1:length(indices)) {
res2@index.range[[i]]<-indices[[i]]@range
res2@index.name[i] <-indices[[i]]@name
res2@index[[i]] <-res@index[[i]]
res2@index.hat[[i]] <-FLQuant(res@index.hat[[i]],
dimnames=dimnames(res@index.hat[[i]]))
res2@index.res[[i]] <-FLQuant(res@index.res[[i]],
dimnames=dimnames(res@index.res[[i]]))
res2@index.var[[i]] <-FLQuant(res@index.var[[i]],
dimnames=dimnames(res@index.var[[i]]))
res2@q.hat[[i]] <-FLQuant(res@q.hat[[i]],
dimnames=dimnames(res@q.hat[[i]]))
res2@q2.hat[[i]] <-FLQuant(res@q2.hat[[i]],
dimnames=dimnames(res@q2.hat[[i]]))
dmns<-dimnames(indices[[i]]@index)
na. <-is.na(indices[[i]]@index[
dmns$age[dmns$age %in% dimnames(res2@index[[i]])$age],
dmns$year[dmns$year %in% dimnames(res2@index[[i]])$year]])
res2@index[[i]][na.] <- NA
res2@index.hat[[i]][na.]<- NA
res2@index.res[[i]][na.]<- NA
res2@index.var[[i]][na.]<- NA
dimnames(res2@index[[i]] )$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@index.hat[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@index.res[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@index.var[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@q.hat[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@q2.hat[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@index[[i]] )$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@index.hat[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@index.res[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@index.var[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@q.hat[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@q2.hat[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@index[[i]] )$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@index.hat[[i]])$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@index.res[[i]])$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@index.var[[i]])$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@q.hat[[i]])$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@q2.hat[[i]])$area <-dimnames(indices[[i]]@index)$area
}
res2@control <- res@control
res2@call <- res@call
res2@desc <- res@desc
res2@range <- stock@range
units(res2@harvest)<-"f"
res2@control@shk.n <-as.logical(res2@control@shk.n)
res2@control@shk.f <-as.logical(res2@control@shk.f)
res2@control@vpa <-as.logical(res2@control@vpa)
return(res2)
}
)
# }}}
# FLXSA.control{{{
FLXSA.control <- function(x=NULL, tol=10e-10, maxit=30, min.nse=0.3, fse=0.5, rage=0, qage=10, shk.n=TRUE,
shk.f=TRUE, shk.yrs=5, shk.ages=5, window=100, tsrange=20, tspower=3, vpa=FALSE){
if (is.null(x)){
res <- new("FLXSA.control", tol=tol, maxit=as.integer(maxit), min.nse=min.nse, fse=fse,
rage=as.integer(rage), qage=as.integer(qage), shk.n=as.logical(shk.n)[1],
shk.f=as.logical(shk.f)[1], shk.yrs=as.integer(shk.yrs), shk.ages=as.integer(shk.ages),
window=as.integer(window), tsrange=as.integer(tsrange), tspower=as.integer(tspower),
vpa=as.logical(vpa)[1])
} else { # We reuse an FLXSA.control object embedded in an FLXSA object
if (!is.null(x) & !(is(x, "FLXSA") | is(x, "FLXSA.control")))
stop("FLXSA must be an 'FLXSA' or an 'FLXSA.control' object!")
if (is(x, "FLXSA"))
res <- x@control
if (is(x, "FLXSA.control"))
res <- x
if (is.null(x))
res <- new("FLXSA.control")
# ... and possibly redefine some of its parameters
if (!missing(tol))
res@tol <- tol
if (!missing(maxit))
res@maxit <- as.integer(maxit)
if (!missing(min.nse))
res@min.nse <- min.nse
if (!missing(fse))
res@fse <- fse
if (!missing(rage))
res@rage <- as.integer(rage)
if (!missing(qage))
res@qage <- as.integer(qage)
if (!missing(shk.n))
res@shk.n <- as.logical(shk.n)[1]
if (!missing(shk.f))
res@shk.f <- as.logical(shk.f)[1]
if (!missing(shk.yrs))
res@shk.yrs <- as.integer(shk.yrs)
if (!missing(shk.ages))
res@shk.ages <- as.integer(shk.ages)
if (!missing(window))
res@window <- as.integer(window)
if (!missing(tsrange))
res@tsrange <- as.integer(tsrange)
if (!missing(tspower))
res@tspower <- as.integer(tspower)
if (!missing(vpa))
res@vpa <- as.logical(vpa)[1]
# Verify that this object is valid
test <- validObject(res)
if (!test)
stop("Invalid object:", test)
}
return(res)
}
# }}}
# assess {{{
setMethod("assess", signature(control="FLXSA.control"),
function(control, stock, indices, ...){
if (!is(stock, "FLStock")) stop("stock not of type FLStock")
if ( is(indices, "FLIndex")) indices<-FLIndices(list(indices))
if (!is(indices, "FLIndices")) stop("indices not of type FLIndices")
print("FLXSA")
FLXSA(stock,indices,control)
}
)
# }}}
# is.FLXSA {{{
is.FLXSA <- function(x)
return(inherits(x, "FLXSA"))
# }}}
# Test if an object is of FLXSA.control class
is.FLXSA.control <- function(x)
return(inherits(x, "FLXSA.control"))
## show (a replacement for print of S3 classes)
setMethod("show", signature(object="FLXSA.control"),
function(object){
n.<-slotNames(object)
for (i in 1:length(n.))
cat(n.[i],"\t\t",slot(object,n.[i]),"\n")
}
)
# DiagsXSA::FLXSA
# Author : MP, KH & JJP, EJ
# Modified to become summary fof FLXSA
# Updated : RDS 11/11/10
# ---------------------------------------------------------------------------------
setGeneric("diagnostics", function(object, ...){
standardGeneric("diagnostics")
}
)
setMethod("diagnostics", signature(object="FLXSA"),
tempfun <- function(object, sections=rep(T, 8), ...){
#browser()
#print(1)
indices<-new("FLIndices")
for (i in 1:length(object@index))
{
indices[[i]] <- FLIndex(index=object@index[[i]])
indices[[i]]@name <- object@index.name[i]
}
control <-object@control
titledat <- paste("FLR XSA Diagnostics ", as.character(Sys.time()),
"\n\nCPUE data from ", object@call[3], "\n\n",
"Catch data for ", dims(object@stock.n)$year," years ",
dims(object@stock.n)$minyear," to ",dims(object@stock.n)$maxyear,
". Ages ",dims(object@stock.n)$min," to ",dims(object@stock.n)$max,".\n\n", sep="")
# print(titledat)
# general tuning series information
idx.info <- NULL
for (i in 1:length(object@index)) {
idx.info <- rbind(idx.info,c(indices[[i]]@name, (dims(object@index[[i]]))$min,
(dims(object@index[[i]]))$max,(dims(object@index[[i]]))$minyear,(dims(object@index[[i]]))$maxyear,
indices[[i]]@range["startf"], indices[[i]]@range["endf"]))
}
dimnames(idx.info) <- list(NULL,c("fleet","first age","last age","first year","last year","alpha","beta"))
# print(as.data.frame(idx.info))
set1 <- paste("\n\n","Time series weights :\n\n")
set2 <- paste(ifelse(control@tsrange==0|control@tspower==0, " Tapered time weighting not applied\n\n",
paste(" Tapered time weighting applied\n", " Power = ",control@tspower,"over ",control@tsrange,
"years\n\n", sep=" ")))
set3 <- "Catchability analysis :\n\n"
set4 <- paste(ifelse(as.numeric(control@rage) < dims(object@stock.n)$min, " Catchability independent of size for all ages\n\n",
paste(" Catchability independent of size for ages > ",control@rage,"\n\n",sep=" ")))
set5 <- paste(ifelse(as.numeric(control@qage) < dims(object@stock.n)$min, " Catchability independent of age for all ages\n\n",
paste(" Catchability independent of age for ages > ",control@qage,"\n\n",sep=" ")))
set6 <- "Terminal population estimation :\n\n"
set7 <- paste(ifelse(control@shk.f, paste(" Survivor estimates shrunk towards the mean F\n",
" of the final ",control@shk.yrs,"years or the ",control@shk.ages,"oldest ages.\n\n",
" S.E. of the mean to which the estimates are shrunk = ", control@fse,"\n",sep=" "),
" Final estimates not shrunk towards mean F\n"))
set8 <- ifelse(as.numeric(control@min.nse)==0, "\n", paste("\n", " Minimum standard error for population\n",
" estimates derived from each fleet = ",control@min.nse,"\n\n", sep=" "))
set9 <- " prior weighting not applied\n\n"
# cat(set1, set2, set3, set4, set5, set6, set7, set8, set9)
regwtstitle <- "Regression weights\n"
### Calculation of time series weighting
yr.range <- max(dims(object@harvest)$minyear,(dims(object@harvest)$maxyear-9)):dims(object@harvest)$maxyear
regWt <- FLQuant(dimnames=list(age = 'all', year = yr.range))
for(y in yr.range)
regWt[,as.character(y)] <- (1-((max(yr.range)-y)/control@tsrange)^control@tspower)^control@tspower
regwts <- matrix(round(regWt,3),dims(regWt)$age,dimnames=list(age="all",year=yr.range))
# cat(regwtstitle)
# print(regwts)
FMtitle <- "\n\n Fishing mortalities\n"
FM <- matrix(round(trim(object@harvest,year=yr.range),3), dims(object@harvest)$age,
dimnames=list(age=dims(object@harvest)$min:dims(object@harvest)$max, year=yr.range))
# cat(FMtitle)
# print(FM)
PNtitle <- "\n\n XSA population number (Thousand)\n"
PN <- (t(matrix(round(trim(object@stock.n,year=yr.range),0), dims(object@stock.n)$age,
dimnames=list(age=dims(object@stock.n)$min:dims(object@stock.n)$max, year=yr.range))))
# cat(PNtitle)
# print(PN)
nextyear <- dims(object@survivors)$maxyear
survtitle <- paste("\n\n Estimated population abundance at 1st Jan ",nextyear,"\n")
survivors <- t(matrix(round(object@survivors[,as.character(nextyear)]),
dimnames=list(age=dims(object@survivors)$min:dims(object@survivors)$max, year=nextyear)))
# cat(survtitle)
# print(survivors)
## tuning info
fleetname <- list()
logQs <- list()
mlqbtitle <- list()
mlqb <- list()
for (f in 1:length(object@index)) {
fleetname[[f]] <- paste("\n\n Fleet: ",indices[[f]]@name,"\n\n","Log catchability residuals.\n\n")
logQs[[f]] <- matrix(round(object@index.res[[f]],3), nrow=dims(object@index.res[[f]])$age,
dimnames=list(age=dimnames(object@index.res[[f]])$age, year=dimnames(object@index.res[[f]])$year))
# print(fleetname[[f]])
# print(logQs[[f]])
if (control@rage < dims(object@index[[f]])$max){
mlqbtitle[[f]] <- paste("\n\n Mean log catchability and standard error of ages with catchability \n",
"independent of year class strength and constant w.r.t. time \n\n")
q.tab <- rbind(Mean_Logq=round(log(object@q.hat[[f]]),4), S.E_Logq=round(sd(matrix(object@index.res[[f]],
dim(object@index.res[[f]])[2],dim(object@index.res[[f]])[1],byrow=T),na.rm=T),4))
colnames(q.tab) <- dimnames(object@q.hat[[f]])$age
if (dims(object@index[[f]])$min <= control@rage ) {
mlqb[[f]] <- q.tab[,as.character((control@rage+1):max(as.numeric(dimnames(object@index[[f]])$age)))]
} else {mlqb[[f]] <- q.tab}
}
# print reg stats powermodel, note that maximum printed age is min(rage, max(age in tun series))
if (dims(object@index[[f]])$min <= control@rage ) {
mlqbtitle[[f]] <- paste("\n Regression statistics \n", "Ages with q dependent on year class strength \n")
mlqb[[f]] <- paste(cbind((matrix(object@q2.hat[[f]][as.character(dims(object@index[[f]])$min:(min(control@rage,dims(object@index[[f]])$max)))],dimnames=list(age=paste("Age ",dims(object@index[[f]])$min:(min(control@rage,dims(object@index[[f]])$max)),sep=""),"slope"))),
(matrix(object@q.hat[[f]][as.character(dims(object@index[[f]])$min:(min(control@rage,dims(object@index[[f]])$max)))],dimnames=list(age=dims(object@index[[f]])$min:(min(control@rage,dims(object@index[[f]])$max)),"intercept")))))
}
# cat(mlqbtitle[[f]])
# print(mlqb[[f]])
}
tysurvtitle <- "\n\n Terminal year survivor and F summaries: \n "
header <- list()
tysurvtext <- list()
#browser()
# cat(tysurvtitle)
agevec <- sort(unique(object@diagnostics$age))
for ( age in 1:length(agevec)){
# cat("age: ", age, "\n")
header[[age]] <- paste("\n ,Age ",agevec[age], " Year class =", max(object@diagnostics$year) - agevec[age] ," \n\n","source \n", sep="")
weights <- object@diagnostics[(object@diagnostics$age==agevec[age]) & (object@diagnostics$year== max(object@diagnostics$year)),]
# calc surivors and scaled wts
weights$survivors <- round(exp(weights$nhat))
weights$scaledWts <- round(weights$w / sum(weights$w) ,3)
row.names(weights) <- weights$source
tysurvtext[[age]] <- weights[ ,c("scaledWts","survivors","yrcls") ]
# cat(header[[age]])
# print(tysurvtext[[age]])
}
## send the text to the screen
if(sections[1]){
cat(titledat)
print(as.data.frame(idx.info))
}
if(sections[2]){
cat(set1, set2, set3, set4, set5, set6, set7, set8, set9)
}
if(sections[3]){
cat(regwtstitle)
print(regwts)
}
if(sections[4]){
cat(FMtitle)
print(FM)
}
if(sections[5]){
cat(PNtitle)
print(PN)
}
if(sections[6]){
cat(survtitle)
print(survivors)
}
if(sections[7]){
for(f in 1:length(object@index)) {
cat(fleetname[[f]])
print(logQs[[f]])
cat(mlqbtitle[[f]])
print(mlqb[[f]])
}
}
if(sections[8]){
cat(tysurvtitle)
for ( age in 1:length(agevec)){
cat(header[[age]])
print(tysurvtext[[age]])
}
}
invisible()
}
)
Try the FLXSA package in your browser
Any scripts or data that you put into this service are public.
FLXSA documentation built on May 2, 2019, 6:06 p.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.
Defines functions validFLXSA.control validFLXSA FLXSA.control is.FLXSA is.FLXSA.control
Documented in FLXSA.control is.FLXSA is.FLXSA.control
# FLXSA.R -
# FLXSA/R/FLXSA.R
# Copyright 2003-2011 FLR Team. Distributed under the GPL 2 or later
# Maintainer: Iago Mosqueira, JRC
# $Id: $
# {{{ FLXSA.control class
validFLXSA.control <- function(object){
if (object@tol <= 0)
return("value of tol must be > 0")
if (object@maxit <= 0)
return("value of maxit must be > 0")
if (object@min.nse < 0)
return("value of min.nse must be > 0")
if (object@fse < 0)
return("value of fse must be > 0")
if (object@rage < -1)
return("value of rage must be >= -1")
if (object@qage < 0)
return("value of qage must be >= 0")
if (object@shk.yrs < 0)
return("value of shk.yrs must be >= 0")
if (object@shk.ages < 0)
return("value of shk.ages must be >= 0")
if (object@window <= 5)
return("value of window must be >= 5")
# Everything is fine
return(TRUE)
}
setClass("FLXSA.control",
representation(
tol ="numeric",
maxit ="integer",
min.nse ="numeric",
fse ="numeric",
rage ="integer",
qage ="integer",
shk.n ="logical",
shk.f ="logical",
shk.yrs ="integer",
shk.ages ="integer",
window ="integer",
tsrange ="integer",
tspower ="integer",
vpa ="logical"),
prototype=prototype(
tol =as.double(10e-10),
maxit =as.integer(30),
min.nse =as.double(0.3),
fse =as.double(0.5),
rage =as.integer(-1),
qage =as.integer(10),
shk.n =TRUE,
shk.f =TRUE,
shk.yrs =as.integer(5),
shk.ages =as.integer(5),
window =as.integer(100),
tsrange =as.integer(20),
tspower =as.integer(3),
vpa =FALSE),
validity=validFLXSA.control
)
# }}}
# FLXSA class {{{
validFLXSA <- function(object){
# All FLQuant objects must have same dimensions
return(TRUE)
Dim <- dim(object@stock.n)
if (!all(dim(object@harvest) == Dim))
return("n and f arrays must have same dimensions")
# Everything is fine
return(TRUE)
}
setClass("FLXSA",
contains='FLAssess',
representation(
survivors="FLQuant",
se.int ="FLQuant",
se.ext ="FLQuant",
n.fshk ="FLQuant",
n.nshk ="FLQuant",
var.fshk ="FLQuant",
var.nshk ="FLQuant",
q.hat ="FLQuants",
q2.hat ="FLQuants",
diagnostics="data.frame",
control ="FLXSA.control"),
prototype=prototype(
survivors=FLQuant(quant='age'),
se.int =FLQuant(quant='age'),
se.ext =FLQuant(quant='age'),
n.fshk =FLQuant(quant='age'),
n.nshk =FLQuant(quant='age'),
var.fshk =FLQuant(quant='age'),
var.nshk =FLQuant(quant='age'),
catch.n =FLQuant(quant='age'),
stock.n =FLQuant(quant='age'),
harvest =FLQuant(quant='age'),
q.hat =FLQuants(),
q2.hat =FLQuants(),
diagnostics=new("data.frame"),
control =new("FLXSA.control")),
validity=validFLXSA
) # }}}
# FLXSA() {{{
setGeneric("FLXSA", function(stock, indices, ...)
standardGeneric("FLXSA"))
setMethod("FLXSA", signature(stock="FLStock", indices="FLIndex"),
function(stock, indices, control=FLXSA.control(), diag.flag=TRUE) {
FLXSA(stock=stock, indices=FLIndices(indices), control=control, diag.flag=diag.flag)
}
)
setMethod("FLXSA", signature(stock="FLStock", indices="FLIndices"),
function(stock, indices, control=FLXSA.control(), diag.flag=TRUE) {
Call <- match.call()
# check FLIndices input
for (i in 1:length(indices)) {
# startf & endf present in @range
if (is.na(indices[[i]]@range["startf"]) || is.na(indices[[i]]@range["endf"]))
stop(paste("Must supply startf & endf for range in FLIndex",i))
# @range has all elements
if (!all(c("min","max","plusgroup","minyear","maxyear","startf","endf") %in%
names(indices[[i]]@range)))
stop("Range must have names 'min','max','plusgroup','minyear','maxyear',
'startf','endf'")
# adjust ranges to dims()
indices[[i]]@range["min"] <- max(indices[[i]]@range["min"], dims(indices[[i]])$min,
stock@range["min"])
indices[[i]]@range["max"] <- min(indices[[i]]@range["max"], dims(indices[[i]])$max,
stock@range["max"])
indices[[i]]@range["minyear"] <- max(indices[[i]]@range["minyear"],
dims(indices[[i]])$minyear, stock@range["minyear"])
indices[[i]]@range["maxyear"] <- min(indices[[i]]@range["maxyear"],
dims(indices[[i]])$maxyear, stock@range["maxyear"])
# trim according to range
age <- indices[[i]]@range["min"]:indices[[i]]@range["max"]
year<- indices[[i]]@range["minyear"]:indices[[i]]@range["maxyear"]
indices[[i]] <- trim(indices[[i]], age=age, year=year)
}
# Double check validity
if (!validObject(stock))
stop("stock is not valid!")
if (!validObject(indices))
stop("FLIndices is not valid!")
if (!validObject(control))
stop("control is not valid!")
# adjust range for minage
if ("minage" %in% names(stock@range))
minage <- stock@range["minage"]
else if ("min" %in% names(stock@range))
minage <- stock@range["min"]
else if ("minquant" %in% names(stock@range))
minage <- stock@range["minquant"]
else
stop("'minage' not found in range")
# adjust range for maxage
if ("maxage" %in% names(stock@range))
maxage <- stock@range["maxage"]
else if ("max" %in% names(stock@range))
maxage <- stock@range["max"]
else if ("maxquant" %in% names(stock@range))
maxage <- stock@range["maxquant"]
else
stop("'maxage' not found in range")
# adjust plsugroup
if ("plusgroup" %in% names(stock@range))
stock@range["plusgroup"] <- maxage
if (maxage<minage | stock@range["maxyear"] < stock@range["minyear"])
stop("Error in range")
if (is.na(stock@range["plusgroup"]))
stop("Plus Group must be specified")
# trim stock
stock <- trim(stock, year=stock@range["minyear"]:stock@range["maxyear"])
if (!is.na(stock@range["plusgroup"]) &
stock@range["plusgroup"] < dims(stock@catch.n)$max)
stock <- setPlusGroup(stock, stock@range["plusgroup"])
stock@m <- stock@m[as.character(minage:maxage),,,,]
# check catch.n is available
if (all(is.na(stock@catch.n)))
stop("catch.n is not available")
stock@stock.n[] <- new('FLQuant')
stock@harvest[] <- new('FLQuant')
# fqs
fqs <- function(assess) {
assess@index <- new("FLQuants", lapply(assess@index,FLQuant))
assess@index.hat <- new("FLQuants", lapply(assess@index.hat,FLQuant))
assess@index.var <- new("FLQuants", lapply(assess@index.var,FLQuant))
assess@index.res <- new("FLQuants", lapply(assess@index.res,FLQuant))
assess@q.hat <- new("FLQuants", lapply(assess@q.hat,FLQuant))
assess@q2.hat <- new("FLQuants", lapply(assess@q2.hat,FLQuant))
if (validObject(assess))
return(assess)
else
stop("not valid")
}
#
iters.stock <-max(unlist(qapply(stock, function(x) dims(x)$iter)))
iters.indices<-max(unlist(lapply(indices@.Data,
function(x) max(unlist(qapply(x, function(x2) dims(x2)$iter))))))
if ((iters.stock>1 && iters.indices>1) && missing(diag.flag))
diag.flag <- FALSE
if ((iters.stock>1 && iters.indices>1) && diag.flag)
return("Multiple iters only allowed if diag.flag=FALSE")
if(!diag.flag) {
res<-.Call("FLXSA", iter(stock, 1), lapply(indices, iter, 1), control, FALSE)
iters <- max(iters.stock,iters.indices)
if (iters>1) {
res@stock.n<-propagate(FLQuant(res@stock.n@.Data),iters)
res@harvest<-propagate(FLQuant(res@harvest@.Data),iters)
for (i in as.character(2:iters)) {
res. <- .Call("FLXSA", iter(stock,i), lapply(indices, iter,i), control, FALSE)
iter(res@stock.n,i)<-FLQuant(res.@stock.n@.Data)
iter(res@harvest,i)<-FLQuant(res.@harvest@.Data)
}
}
res@harvest@units <- "f"
return(res)
}
res <-.Call("FLXSA", stock, indices, control, diag.flag)
res <-fqs(.Call("FLXSA", stock, indices, control, diag.flag))
if (class(res) != "FLXSA")
return(res)
res@call <- as.character(Call)
# put wts amd nhats into a data frame
df <- as.data.frame(res@wts)
df1 <- (df[4])
df1[df1 >= 1, 1] <- paste("index", df1[df1 >= 1, 1])
df1[df1 == -1, 1] <- "fshk"
df1[df1 == -2, 1] <- "nshk"
df <- cbind(df[-4], df1)
names(df) <- c("w", "nhat", "yrcls", "age", "year", "source")
for(i in 1:length(indices)) {
v <- paste("index", i)
if (length(slot(indices[[i]],"name"))>0)
df$source[df$source==v] <- slot(indices[[i]],"name")
}
wts.df <-df[,c(4,5,1,6)]
names(wts.df)[3] <-"data"
nhat.df<-df[,c(4,5,2,6)]
names(nhat.df)[3]<-"data"
index.names<-sort(unique(df[,"source"]))
index.names<-index.names[substr(index.names,1,5)=="index"]
fill.flq<-function(obj) {
dms <-dims(obj)
dmns<-dimnames(obj)
dmns[[1]]<-dms[[2]]:dms[[3]]
dmns[[2]]<-dms[[5]]:dms[[6]]
res<-as.FLQuant(NA,dimnames=dmns)
res[dimnames(obj)[[1]],dimnames(obj)[[2]],,,]<-obj[dimnames(obj)[[1]],
dimnames(obj)[[2]],,,]
return(res)
}
res2 <- new("FLXSA")
res2@index.var<-new('FLQuants')
res2@index.hat<-new('FLQuants')
res2@index <-new('FLQuants')
j=0
for (i in index.names) {
j=j+1
res2@index.var[[j]]<-1.0/fill.flq(as.FLQuant(wts.df[wts.df[, "source"]==i,-4]))
res2@index.hat[[j]]<- fill.flq(as.FLQuant(nhat.df[nhat.df[,"source"]==i,-4]))
res2@index.var[[j]]<-1.0/as.FLQuant(wts.df[wts.df[, "source"]==i,-4])
res2@index.hat[[j]]<-as.FLQuant(nhat.df[nhat.df[,"source"]==i,-4])
dmns <- dimnames(res2@index.hat[[j]])
index <- trim(indices[[j]]@index,age=dmns$age,year=dmns$year)
res2@index[ is.na(index)]<-NA
res2@index.hat[is.na(index)]<-NA
res2@index.var[is.na(index)]<-NA
res2@index.res[is.na(index)]<-NA
}
# F shrinkage
fshk <- df[df[,"source"]=="fshk",]
if (length(fshk[,1])>0) {
y.range <- range(fshk[,"year"])
a.range<-range(fshk[,"age"])
max.yr <-fshk[fshk[,"age"] ==a.range[2],]
max.age<-fshk[fshk[,"year"]==y.range[2],]
res2@n.fshk <-as.FLQuant(NA,dimnames=list(age=a.range[1]:a.range[2],
year=y.range[1]:y.range[2]))
res2@n.fshk[as.character(max.age[,"age"]),as.character(y.range[2])]<-max.age[,
"nhat"]
res2@n.fshk[as.character(a.range[2]),as.character(max.yr[,"year"])]<-max.yr[,"nhat"]
res2@var.nshk <-as.FLQuant(NA,dimnames=list(age=a.range[1]:a.range[2],
year=y.range[1]:y.range[2]))
res2@var.nshk[as.character(max.age[,"age"]),as.character(y.range[2])]<-
1/max.age[,"w"]
res2@var.nshk[ as.character(a.range[2]),as.character(max.yr[,"year"])]<-
1/max.yr[,"w"]
}
# N shrinkage
nshk <-df[df[,"source"]=="nshk",]
if (length(nshk[,1])>0) {
y.range<-range(nshk[,"year"])
a.range<-range(nshk[,"age"])
max.yr <-nshk[nshk[,"age"] ==a.range[2],]
max.age<-nshk[nshk[,"year"]==y.range[2],]
res2@n.nshk <-as.FLQuant(NA,dimnames=list(age=a.range[1]:a.range[2],
year=y.range[1]:y.range[2]))
res2@n.nshk[,as.character(y.range[2])] <- max.age[,"nhat"]
res2@n.nshk[ as.character(a.range[2])] <- max.yr[,"nhat"]
res2@var.nshk <- as.FLQuant(NA,dimnames=list(age=a.range[1]:a.range[2],
year=y.range[1]:y.range[2]))
res2@var.nshk[,as.character(y.range[2])]<-1/max.age[,"w"]
res2@var.nshk[ as.character(a.range[2])]<-1/max.yr[,"w"]
}
res2@diagnostics<- df
res2@index.hat <- res@index.hat
res2@stock.n <- FLQuant(res@stock.n@.Data)
res2@harvest <- FLQuant(res@harvest@.Data)
res2@survivors <- FLQuant(res@survivors@.Data)
res2@se.int <- FLQuant(res@se.int@.Data)
res2@se.ext <- FLQuant(res@se.ext@.Data)
res2@n.fshk <- FLQuant(res@n.fshk@.Data)
res2@n.nshk <- FLQuant(res@n.nshk@.Data)
res2@var.fshk <- FLQuant(res@var.fshk@.Data)
res2@var.nshk <- FLQuant(res@var.nshk@.Data)
res2@harvest@units <- "f"
res2@q.hat <- res@q.hat
res2@q2.hat <- res@q2.hat
res2@index.res <- res@index.res
res2@index <- res@index
res2@index.hat <-FLQuants(res@index.hat)
res2@index.res <-FLQuants(res@index.res)
res2@index.var <-FLQuants(res@index.var)
res2@q.hat <-FLQuants(res@q.hat)
res2@q2.hat <-FLQuants(res@q2.hat)
for (i in 1:length(indices)) {
res2@index.range[[i]]<-indices[[i]]@range
res2@index.name[i] <-indices[[i]]@name
res2@index[[i]] <-res@index[[i]]
res2@index.hat[[i]] <-FLQuant(res@index.hat[[i]],
dimnames=dimnames(res@index.hat[[i]]))
res2@index.res[[i]] <-FLQuant(res@index.res[[i]],
dimnames=dimnames(res@index.res[[i]]))
res2@index.var[[i]] <-FLQuant(res@index.var[[i]],
dimnames=dimnames(res@index.var[[i]]))
res2@q.hat[[i]] <-FLQuant(res@q.hat[[i]],
dimnames=dimnames(res@q.hat[[i]]))
res2@q2.hat[[i]] <-FLQuant(res@q2.hat[[i]],
dimnames=dimnames(res@q2.hat[[i]]))
dmns<-dimnames(indices[[i]]@index)
na. <-is.na(indices[[i]]@index[
dmns$age[dmns$age %in% dimnames(res2@index[[i]])$age],
dmns$year[dmns$year %in% dimnames(res2@index[[i]])$year]])
res2@index[[i]][na.] <- NA
res2@index.hat[[i]][na.]<- NA
res2@index.res[[i]][na.]<- NA
res2@index.var[[i]][na.]<- NA
dimnames(res2@index[[i]] )$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@index.hat[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@index.res[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@index.var[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@q.hat[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@q2.hat[[i]])$unit <-dimnames(indices[[i]]@index)$unit
dimnames(res2@index[[i]] )$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@index.hat[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@index.res[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@index.var[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@q.hat[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@q2.hat[[i]])$season <-dimnames(indices[[i]]@index)$season
dimnames(res2@index[[i]] )$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@index.hat[[i]])$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@index.res[[i]])$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@index.var[[i]])$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@q.hat[[i]])$area <-dimnames(indices[[i]]@index)$area
dimnames(res2@q2.hat[[i]])$area <-dimnames(indices[[i]]@index)$area
}
res2@control <- res@control
res2@call <- res@call
res2@desc <- res@desc
res2@range <- stock@range
units(res2@harvest)<-"f"
res2@control@shk.n <-as.logical(res2@control@shk.n)
res2@control@shk.f <-as.logical(res2@control@shk.f)
res2@control@vpa <-as.logical(res2@control@vpa)
return(res2)
}
)
# }}}
# FLXSA.control{{{
FLXSA.control <- function(x=NULL, tol=10e-10, maxit=30, min.nse=0.3, fse=0.5, rage=0, qage=10, shk.n=TRUE,
shk.f=TRUE, shk.yrs=5, shk.ages=5, window=100, tsrange=20, tspower=3, vpa=FALSE){
if (is.null(x)){
res <- new("FLXSA.control", tol=tol, maxit=as.integer(maxit), min.nse=min.nse, fse=fse,
rage=as.integer(rage), qage=as.integer(qage), shk.n=as.logical(shk.n)[1],
shk.f=as.logical(shk.f)[1], shk.yrs=as.integer(shk.yrs), shk.ages=as.integer(shk.ages),
window=as.integer(window), tsrange=as.integer(tsrange), tspower=as.integer(tspower),
vpa=as.logical(vpa)[1])
} else { # We reuse an FLXSA.control object embedded in an FLXSA object
if (!is.null(x) & !(is(x, "FLXSA") | is(x, "FLXSA.control")))
stop("FLXSA must be an 'FLXSA' or an 'FLXSA.control' object!")
if (is(x, "FLXSA"))
res <- x@control
if (is(x, "FLXSA.control"))
res <- x
if (is.null(x))
res <- new("FLXSA.control")
# ... and possibly redefine some of its parameters
if (!missing(tol))
res@tol <- tol
if (!missing(maxit))
res@maxit <- as.integer(maxit)
if (!missing(min.nse))
res@min.nse <- min.nse
if (!missing(fse))
res@fse <- fse
if (!missing(rage))
res@rage <- as.integer(rage)
if (!missing(qage))
res@qage <- as.integer(qage)
if (!missing(shk.n))
res@shk.n <- as.logical(shk.n)[1]
if (!missing(shk.f))
res@shk.f <- as.logical(shk.f)[1]
if (!missing(shk.yrs))
res@shk.yrs <- as.integer(shk.yrs)
if (!missing(shk.ages))
res@shk.ages <- as.integer(shk.ages)
if (!missing(window))
res@window <- as.integer(window)
if (!missing(tsrange))
res@tsrange <- as.integer(tsrange)
if (!missing(tspower))
res@tspower <- as.integer(tspower)
if (!missing(vpa))
res@vpa <- as.logical(vpa)[1]
# Verify that this object is valid
test <- validObject(res)
if (!test)
stop("Invalid object:", test)
}
return(res)
}
# }}}
# assess {{{
setMethod("assess", signature(control="FLXSA.control"),
function(control, stock, indices, ...){
if (!is(stock, "FLStock")) stop("stock not of type FLStock")
if ( is(indices, "FLIndex")) indices<-FLIndices(list(indices))
if (!is(indices, "FLIndices")) stop("indices not of type FLIndices")
print("FLXSA")
FLXSA(stock,indices,control)
}
)
# }}}
# is.FLXSA {{{
is.FLXSA <- function(x)
return(inherits(x, "FLXSA"))
# }}}
# Test if an object is of FLXSA.control class
is.FLXSA.control <- function(x)
return(inherits(x, "FLXSA.control"))
## show (a replacement for print of S3 classes)
setMethod("show", signature(object="FLXSA.control"),
function(object){
n.<-slotNames(object)
for (i in 1:length(n.))
cat(n.[i],"\t\t",slot(object,n.[i]),"\n")
}
)
# DiagsXSA::FLXSA
# Author : MP, KH & JJP, EJ
# Modified to become summary fof FLXSA
# Updated : RDS 11/11/10
# ---------------------------------------------------------------------------------
setGeneric("diagnostics", function(object, ...){
standardGeneric("diagnostics")
}
)
setMethod("diagnostics", signature(object="FLXSA"),
tempfun <- function(object, sections=rep(T, 8), ...){
#browser()
#print(1)
indices<-new("FLIndices")
for (i in 1:length(object@index))
{
indices[[i]] <- FLIndex(index=object@index[[i]])
indices[[i]]@name <- object@index.name[i]
}
control <-object@control
titledat <- paste("FLR XSA Diagnostics ", as.character(Sys.time()),
"\n\nCPUE data from ", object@call[3], "\n\n",
"Catch data for ", dims(object@stock.n)$year," years ",
dims(object@stock.n)$minyear," to ",dims(object@stock.n)$maxyear,
". Ages ",dims(object@stock.n)$min," to ",dims(object@stock.n)$max,".\n\n", sep="")
# print(titledat)
# general tuning series information
idx.info <- NULL
for (i in 1:length(object@index)) {
idx.info <- rbind(idx.info,c(indices[[i]]@name, (dims(object@index[[i]]))$min,
(dims(object@index[[i]]))$max,(dims(object@index[[i]]))$minyear,(dims(object@index[[i]]))$maxyear,
indices[[i]]@range["startf"], indices[[i]]@range["endf"]))
}
dimnames(idx.info) <- list(NULL,c("fleet","first age","last age","first year","last year","alpha","beta"))
# print(as.data.frame(idx.info))
set1 <- paste("\n\n","Time series weights :\n\n")
set2 <- paste(ifelse(control@tsrange==0|control@tspower==0, " Tapered time weighting not applied\n\n",
paste(" Tapered time weighting applied\n", " Power = ",control@tspower,"over ",control@tsrange,
"years\n\n", sep=" ")))
set3 <- "Catchability analysis :\n\n"
set4 <- paste(ifelse(as.numeric(control@rage) < dims(object@stock.n)$min, " Catchability independent of size for all ages\n\n",
paste(" Catchability independent of size for ages > ",control@rage,"\n\n",sep=" ")))
set5 <- paste(ifelse(as.numeric(control@qage) < dims(object@stock.n)$min, " Catchability independent of age for all ages\n\n",
paste(" Catchability independent of age for ages > ",control@qage,"\n\n",sep=" ")))
set6 <- "Terminal population estimation :\n\n"
set7 <- paste(ifelse(control@shk.f, paste(" Survivor estimates shrunk towards the mean F\n",
" of the final ",control@shk.yrs,"years or the ",control@shk.ages,"oldest ages.\n\n",
" S.E. of the mean to which the estimates are shrunk = ", control@fse,"\n",sep=" "),
" Final estimates not shrunk towards mean F\n"))
set8 <- ifelse(as.numeric(control@min.nse)==0, "\n", paste("\n", " Minimum standard error for population\n",
" estimates derived from each fleet = ",control@min.nse,"\n\n", sep=" "))
set9 <- " prior weighting not applied\n\n"
# cat(set1, set2, set3, set4, set5, set6, set7, set8, set9)
regwtstitle <- "Regression weights\n"
### Calculation of time series weighting
yr.range <- max(dims(object@harvest)$minyear,(dims(object@harvest)$maxyear-9)):dims(object@harvest)$maxyear
regWt <- FLQuant(dimnames=list(age = 'all', year = yr.range))
for(y in yr.range)
regWt[,as.character(y)] <- (1-((max(yr.range)-y)/control@tsrange)^control@tspower)^control@tspower
regwts <- matrix(round(regWt,3),dims(regWt)$age,dimnames=list(age="all",year=yr.range))
# cat(regwtstitle)
# print(regwts)
FMtitle <- "\n\n Fishing mortalities\n"
FM <- matrix(round(trim(object@harvest,year=yr.range),3), dims(object@harvest)$age,
dimnames=list(age=dims(object@harvest)$min:dims(object@harvest)$max, year=yr.range))
# cat(FMtitle)
# print(FM)
PNtitle <- "\n\n XSA population number (Thousand)\n"
PN <- (t(matrix(round(trim(object@stock.n,year=yr.range),0), dims(object@stock.n)$age,
dimnames=list(age=dims(object@stock.n)$min:dims(object@stock.n)$max, year=yr.range))))
# cat(PNtitle)
# print(PN)
nextyear <- dims(object@survivors)$maxyear
survtitle <- paste("\n\n Estimated population abundance at 1st Jan ",nextyear,"\n")
survivors <- t(matrix(round(object@survivors[,as.character(nextyear)]),
dimnames=list(age=dims(object@survivors)$min:dims(object@survivors)$max, year=nextyear)))
# cat(survtitle)
# print(survivors)
## tuning info
fleetname <- list()
logQs <- list()
mlqbtitle <- list()
mlqb <- list()
for (f in 1:length(object@index)) {
fleetname[[f]] <- paste("\n\n Fleet: ",indices[[f]]@name,"\n\n","Log catchability residuals.\n\n")
logQs[[f]] <- matrix(round(object@index.res[[f]],3), nrow=dims(object@index.res[[f]])$age,
dimnames=list(age=dimnames(object@index.res[[f]])$age, year=dimnames(object@index.res[[f]])$year))
# print(fleetname[[f]])
# print(logQs[[f]])
if (control@rage < dims(object@index[[f]])$max){
mlqbtitle[[f]] <- paste("\n\n Mean log catchability and standard error of ages with catchability \n",
"independent of year class strength and constant w.r.t. time \n\n")
q.tab <- rbind(Mean_Logq=round(log(object@q.hat[[f]]),4), S.E_Logq=round(sd(matrix(object@index.res[[f]],
dim(object@index.res[[f]])[2],dim(object@index.res[[f]])[1],byrow=T),na.rm=T),4))
colnames(q.tab) <- dimnames(object@q.hat[[f]])$age
if (dims(object@index[[f]])$min <= control@rage ) {
mlqb[[f]] <- q.tab[,as.character((control@rage+1):max(as.numeric(dimnames(object@index[[f]])$age)))]
} else {mlqb[[f]] <- q.tab}
}
# print reg stats powermodel, note that maximum printed age is min(rage, max(age in tun series))
if (dims(object@index[[f]])$min <= control@rage ) {
mlqbtitle[[f]] <- paste("\n Regression statistics \n", "Ages with q dependent on year class strength \n")
mlqb[[f]] <- paste(cbind((matrix(object@q2.hat[[f]][as.character(dims(object@index[[f]])$min:(min(control@rage,dims(object@index[[f]])$max)))],dimnames=list(age=paste("Age ",dims(object@index[[f]])$min:(min(control@rage,dims(object@index[[f]])$max)),sep=""),"slope"))),
(matrix(object@q.hat[[f]][as.character(dims(object@index[[f]])$min:(min(control@rage,dims(object@index[[f]])$max)))],dimnames=list(age=dims(object@index[[f]])$min:(min(control@rage,dims(object@index[[f]])$max)),"intercept")))))
}
# cat(mlqbtitle[[f]])
# print(mlqb[[f]])
}
tysurvtitle <- "\n\n Terminal year survivor and F summaries: \n "
header <- list()
tysurvtext <- list()
#browser()
# cat(tysurvtitle)
agevec <- sort(unique(object@diagnostics$age))
for ( age in 1:length(agevec)){
# cat("age: ", age, "\n")
header[[age]] <- paste("\n ,Age ",agevec[age], " Year class =", max(object@diagnostics$year) - agevec[age] ," \n\n","source \n", sep="")
weights <- object@diagnostics[(object@diagnostics$age==agevec[age]) & (object@diagnostics$year== max(object@diagnostics$year)),]
# calc surivors and scaled wts
weights$survivors <- round(exp(weights$nhat))
weights$scaledWts <- round(weights$w / sum(weights$w) ,3)
row.names(weights) <- weights$source
tysurvtext[[age]] <- weights[ ,c("scaledWts","survivors","yrcls") ]
# cat(header[[age]])
# print(tysurvtext[[age]])
}
## send the text to the screen
if(sections[1]){
cat(titledat)
print(as.data.frame(idx.info))
}
if(sections[2]){
cat(set1, set2, set3, set4, set5, set6, set7, set8, set9)
}
if(sections[3]){
cat(regwtstitle)
print(regwts)
}
if(sections[4]){
cat(FMtitle)
print(FM)
}
if(sections[5]){
cat(PNtitle)
print(PN)
}
if(sections[6]){
cat(survtitle)
print(survivors)
}
if(sections[7]){
for(f in 1:length(object@index)) {
cat(fleetname[[f]])
print(logQs[[f]])
cat(mlqbtitle[[f]])
print(mlqb[[f]])
}
}
if(sections[8]){
cat(tysurvtitle)
for ( age in 1:length(agevec)){
cat(header[[age]])
print(tysurvtext[[age]])
}
}
invisible()
}
)
Try the FLXSA package in your browser
Any scripts or data that you put into this service are public.
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.