FLR

FS#45: FLPar & pmin

Laurie Kell — 2011-03-24T08:42:43Z

FLPar doesn´t work with all methods e.g.

this works
pmin(1:2/FLPar(1),0)
[1] 0 0

but this doesn´t

pmin(1/FLPar(1:2),0)
Error in checkSlotAssignment(object, name, value) :
assignment of an object of class "logical" is not valid for slot ".Data" in an object of class "FLPar"; is(value, "array") is not TRUE
Calls: pmin -> > -> > -> @<- -> slot<- -> checkSlotAssignment

FS#44: covar in FLSR

Laurie Kell — 2011-01-25T10:50:20Z

rickerCovA <- function(){
logl <- function(a, b, c, rec, ssb, covar){
loglAR1(log(rec), log(a*(1+c*covar[[1]])*ssb*exp(-b*ssb)))}

initial <- structure(function(rec, ssb, covar) {
# The function to provide initial values
res <-coefficients(lm(c(log(rec/ssb))~c(ssb)))
return(FLPar(a=max(exp(res[1])), b=-max(res[2]), c=0.0))},

# lower and upper limits for optim()
lower=rep(-Inf, 3),
upper=rep( Inf, 3))

model <- rec~a*(1+c*covar[[1]])*ssb*exp(-b*ssb)

return(list(logl=logl, model=model, initial=initial))}

#### Modified so temperature affects larval K
rickerCovB <- function(){
logl <- function(a, b, c, rec, ssb, covar){
loglAR1(log(rec), log(a*(1+c*covar[[1]])*ssb*exp(-b*ssb)))}

initial <- structure(function(rec, ssb, covar) {
# The function to provide initial values
res <-coefficients(lm(c(log(rec/ssb))~c(ssb)))
return(FLPar(a=max(exp(res[1])), b=-max(res[2]), c=0.0))},

# lower and upper limits for optim()
lower=rep(-Inf, 3),
upper=rep( Inf, 3))

model <- rec~a*ssb*exp(-b*(1+c*covar[[1]])*ssb)

return(list(logl=logl, model=model, initial=initial))}

nao <-read.table("http://www.cdc.noaa.gov/data/correlation/nao.data", skip=1, nrow=62, na.strings="-99.90")
dnms <-list(quant="nao", year=1948:2009, unit="unique", season=1:12, area="unique")
nao <-FLQuant(unlist(nao[,-1]), dimnames=dnms, units="nao")

#### include NAO as covar (note that it must be a FLQuants with a single component
#### called “covar” that matches the year span of the data) and adjust the model.

herSR <-as.FLSR(her4)
model(herSR) <-ricker()
herSR <-fmle(herSR)

herCovA <-as.FLSR(her4)
herCovA <-transform(herCovA,ssb=ssb/1000,rec=rec/1000)
model(herCovA)<-rickerCovA()
covar(herCovA)<-FLQuants(covar=seasonMeans(trim(nao, year=dimnames(ssb(herCovA))$year)))
herCovA <-fmle(herCovA,fixed=list(c=0))

FS#43: plot(FLStock)

Laurie Kell — 2011-01-25T10:48:30Z

Issues a warning for a valid plot, i.e.
data(ple4)
plot(ple4)

this is because coercion of data.frame for rec & ssb does

rbind(as.data.frame(ssb(ple4)),as.data.frame(rec(ple4)))

although

rbind(as.data.frame(rec(ple4)),as.data.frame(ssb(ple4)))

doesn´t

FS#42: Profile fails if parameter negative

Laurie Kell — 2010-12-07T09:02:26Z

model(nsher)<-bevholt()
nsher <-fmle(nsher,fixed=list(b=-0.1))
profile(nsher)

FS#41: plotting

Laurie Kell — 2010-10-26T09:44:51Z

Often FLStock::plot or FLStocks::plot is used to summarise projections for different management measures and stock hypotheses. However, in most cases the y-axis ranges arn´t comparable making it difficult to compare results. For example the y-axis don´t start at 0 and if a stock recovers to different levels in scenarios (the whole point of running them!) then the y-axes will differ across scenarios. Since you can´t specify the max y values.

While if you have an FLstock object with iterations the range is not set by the interquartile range but the actual data e.g.

plot(ple4)

stock.n(ple4[[1]])<-propagate(stock.n(ple4[[1]]),500)
stock.n(ple4[[1]])[1,,,,,1]<-stock.n(ple4[[1]])[1,,,,,1]*10
plot(ple4[[1]])

FS#40: FLModel logliklihood and gradient functions are different

Laurie Kell — 2010-10-26T09:38:47Z

In FLModel there is a slot for the gradient, which in FLSR is by default empty. However if you can provide this then you can both speed up fitting and make it more robust.

However, the syntax for the logliklihood and gradient functions are different e.g.

obj<-as.FLSR(ple4,model=bevholt)

bevholt()$logl

gr(obj)<-function(a,b,rec,ssb) {params(obj)["a"]<-a;params(obj)["b"]<-b;FLPar(computeD(obj)[,1:2])}
obj<-fmle(obj)

gr(obj)<-function(par,rec,ssb) {params(obj)["a"]<-par[1];params(obj)["b"]<-par[2];FLPar(computeD(obj)[,1:2])}
obj<-fmle(obj)

The syntax for logl and gr should be the same since If there is a good reason for specifying the parameters separately for logl then the same logic must also apply to gr. However, I think specifying the parameters separately (as currently done in logl) implies that FLPar is redundant. Therefore can we change all functions/methods in FLModel/FLSR to use FLPar rather than the parameters seperately.

This is not a trivial point as if we don´t fix this now we will run into bigger problems in the future, I have already found this when working with ADMb & ADol-C.

Laurie

FS#39: FLPar doesn´t print correctly

Laurie Kell — 2010-10-18T11:40:28Z

FLPar doesn´t print correctly e.g.

> FLPar(array(1:6,c(3,2,1)))
An object of class "FLPar"

params 1 2
a 2 5
b 5 2
c 2 5
units: NA NA NA
> FLPar(array(1:6,c(3,2,1)))@.Data
, , iter = 1

params 1 2
a 1 4
b 2 5
c 3 6

Laurie

FS#38: Broken FLxxx creators

Laurie Kell — 2010-10-01T14:33:59Z

The FLXXX creators are broken

> data(ple4.indices)
> data(ple4)
> pleXSA<-FLXSA(ple4,ple4.indices)
[1] 0
[1] 0
Error in validityMethod(as(object, superClass)) :
Not all 'quant' names are the same. Check using qapply(x, quant)
> qapply(FLAssess(), quant)
$catch.n
[1] "quant"

$stock.n
[1] "quant"

$harvest
[1] "quant"

> qapply(FLQuant(), quant)
Error in function (classes, fdef, mtable) :
unable to find an inherited method for function "qapply", for signature "FLQuant", "standardGeneric"
> FLQuant()
An object of class "FLQuant"
, , unit = unique, season = all, area = unique

year
quant 1
all NA

units: NA

FS#37: rec method

Laurie Kell — 2010-09-30T12:04:18Z

The rec method has a bug, since rec.age is numeric therefore if it is anything other than 1 then you will get the wrong recruitment.

setMethod('rec', signature(object='FLStock'),
# function(object, rec.age=ac(dims(object)$min))
function(object, rec.age=object@range["min"])
{
if(dims(object)$quant != 'age')
stop("rec(FLStock) only defined for age-based objects")
res <- quantSums(stock.n(object)[rec.age,])
dimnames(res) <- list(age=rec.age)
return(res)
}
) # }}}

FS#36: Update to wiki page and R documentation

Alex Hanke — 2010-09-22T19:52:54Z

In the wiki tutorial Getting Data into FLR the bullet :'ICCAT VPA2Box format' could become a section heading with the following info: The readFLStock() function can be used to read in stock and CPUE data used to conduct stock assessment. The 'ICCAT Adapt file' format is supported by specifying type='Adapt'. The 'Adapt file format' consists of a single file of input data containing the biological parameters, catches and catch per unit effort data and it typically carries the *.d01 file extension. Sample Code: path2<-'C:\\Bluefin CPUE\\ICCAT\\Assessment\\Analyses\\West VPA\\Run3d\\bftw7009.d01' BFTStock<-readFLStock(file=path2,type='Adapt') In the R documentation for 'IOfunctions {FLCore} : Input/Output of FLR objects, it would help to see an example like the one above and perhaps use the same text in the details section. Finally, the comments for the 'type' argument need to be corrected so that 'adapt' is capitalized. I was stumped for a while trying to figure out why my code would not work and had to look at the source code for readFLStock before I guessed that my type argument should have been "Adapt" not "adapt". Regards, -Alex-