Nothing
# FLSR - Stock-recruitment relationships
# FLCore/R/FLSR.R
# Copyright 2003-2012 FLR Team. Distributed under the GPL 2 or later
# Maintainer: Dorleta GarcĂa, AZTI
# $Id: FLSR.R 1778 2012-11-23 08:43:57Z imosqueira $
# FLSR() {{{
setMethod('FLSR', signature(model='ANY'),
function(model, ...)
{
# TODO if no proper rec.age
args <- list(...)
# If both rec and ssb given
if(all(c('rec', 'ssb') %in% names(args)))
{
res <- FLModel(model, ..., class='FLSR')
}
# if rec given, then ssb is dims$min years less
else if ('rec' %in% names(args))
{
drec <- dims(args[['rec']])
ssb <- FLQuant(dimnames=dimnames(window(args[['rec']],
start=drec$minyear-drec$min, end=drec$maxyear-drec$min)))
res <- FLModel(model, ssb=ssb, ..., class='FLSR')
}
# ssb
else if ('ssb' %in% names(args))
{
dssb <- dims(args[['ssb']])
rec <- FLQuant(dimnames=dimnames(window(args[['ssb']],
start=dssb$minyear-1, end=dssb$maxyear-1)))
res <- FLModel(model, rec=rec, ..., class='FLSR')
}
else
res <- FLModel(model, ..., class='FLSR')
# check if years in 'rec' and 'ssb' dimnames match with 'rec' age
if(isTRUE(try(dims(rec(res))$minyear - dims(ssb(res))$minyear != dims(rec(res))$min)))
warning("year dimnames for 'rec' and 'ssb' do not match with recruitment age")
return(res)
}
)
setMethod('FLSR', signature(model='missing'),
function(...)
return(FLSR(formula(NULL), ...))) # }}}
# sr() {{{
sr <- function(sr, ...)
{
# if logl present, run fmle()
try <- try(fmle(sr), silent=TRUE, ...)
if(is(try, 'try-error'))
# else if model present run nls()
try <- try(nls(sr), silent=TRUE, ...)
if(is(try, 'try-error'))
stop('neither nls() nor fmle() could be run on this object')
return(try)
} #}}}
## as.FLSR {{{
setMethod("as.FLSR", signature(object="FLStock"),
function(object, rec.age = dims(stock.n(object))$min, ...)
{
# check rec.age
if(rec.age < dims(stock.n(object))$min)
stop("Supplied recruitment age less than minimum age class")
args <- list(...)
slots <- names(args)[ifelse(length(which(names(args) == "rec.age"))>0,
-which(names(args) == "rec.age"), 1:length(args))]
# calculate ssb and create FLSR object incorprating rec.age
rec <- object@stock.n[as.character(rec.age),]
ssb <- ssb(object)
# now alter stock and recruitment to factor in the recruitement age
if((dim(rec)[2]-1) <= rec.age)
stop("FLStock recruitment data set too short")
rec <- rec[,(1+rec.age):dim(rec)[2]]
units(rec) <- units(slot(object, "stock.n"))
ssb <- ssb[,1:(dim(ssb)[2] - rec.age)]
units(ssb) <- units(slot(object, "stock.wt"))
# create the FLSR object
sr <- FLSR(rec=rec, ssb=ssb, name=object@name,
fitted = FLQuant(dimnames = dimnames(rec), units=units(rec)),
residuals = FLQuant(dimnames = dimnames(rec)),
desc = "'rec' and 'ssb' slots obtained from a 'FLStock' object", ...)
validObject(sr)
return(sr)
}
)
setMethod("as.FLSR", signature(object="FLBiol"),
function(object, rec.age = "missing", ...)
{
validObject(object)
# recruitment delay set using minage
# from the FLStock object
if(missing(rec.age))
rec.age <- dims(n(object))$min
else{
if(rec.age < dims(n(object))$min)
stop("Supplied recruitment age less than minimum age class")
}
if (all(is.na(slot(object, "n"))) || all(is.na(slot(object, "wt"))) ||
all(is.na(slot(object, "fec"))) || all(is.na(slot(object, "spwn"))))
stop("biol must have 'n', 'wt', 'm', 'fec' and 'spwn'")
args <- list(...)
slots <- names(args)[ifelse(length(which(names(args) == "rec.age"))>0,
-which(names(args) == "rec.age"), 1:length(args))]
# calculate ssb and create FLSR object incorporating rec.age
rec <- dimSums(object@n[as.character(rec.age),])
ssb <- apply(slot(object, "n") * exp(- slot(object, "m")*slot(object, "spwn"))*
slot(object, "wt")*slot(object, "fec"), 2:5, sum)
# now alter the stock and recruitment
# vectors to factor in the recruitement age
if((dim(rec)[2]-1) <= rec.age)
stop("FLBiol recruitment data set too short")
rec <- rec[,(1+rec.age):dim(rec)[2],,,]
ssb <- ssb[,1:(dim(ssb)[2] - rec.age),,,]
# create the FLSR object
sr <- FLSR(rec = rec,ssb = ssb, name = object@name,
desc = "'rec' and 'ssb' slots obtained from a 'FLBiol' object", ...)
slot(sr, "fitted") <- FLQuant(dimnames = dimnames(slot(sr, "rec")))
slot(sr, "residuals") <- FLQuant(dimnames = dimnames(slot(sr, "rec")))
units(slot(sr, "rec")) <- units(slot(object, "n"))
units(slot(sr, "ssb")) <- units(slot(object, "wt"))
units(slot(sr, "fitted")) <- units(slot(sr, "rec"))
return(sr)
}
) # }}}
# plot {{{
setMethod("plot", signature(x="FLSR", y="missing"),
function(x, main="Functional form", log.resid=FALSE, cex=0.8)
{
years <- as.numeric(dimnames(residuals(x))$year)
scales <- list(y=list(rot=90), tck=c(1,0))
# initial device settings
trellis.device(new=FALSE)
trellis.par.set(list(layout.heights = list(bottom.padding = -0.3,
axis.xlab.padding = 0.3, xlab = -0.3), layout.widths = list(left.padding = -0.3,
right.padding = -0.3, ylab.axis.padding = -0.3)))
# panel functions
srpanel <- function(x, y, ...) {
panel.xyplot(x, y, col='black', cex=cex)
panel.loess(x,y, col='blue', lty=4)
panel.abline(a=0, b=0, lty=2, col='gray60')
}
respanel <- function(x, y, ...) {
panel.xyplot(x, y, col='black', cex=cex)
panel.lmline(x, y, ..., col='red')
panel.abline(a=0, b=0, lty=2, col='gray60')
}
# get dimensions to condition on (skip quant)
condnames <- names((fitted(x)))[c(3:5)][dim(fitted(x))[c(3:5)]!=1]
cond <- paste(condnames, collapse="+")
if(cond != "") cond <- paste("|", cond)
# 1. SR values with fitted curve
ssb <- FLQuant(seq(0, max(ssb(x) + (ssb(x)*0.15), na.rm=TRUE), length=dim(ssb(x))[2]),
dimnames=dimnames(ssb(x))[1:5])
fitted <- predict(x, ssb=ssb)
print(xyplot(formula(paste("fitted~ssb", cond)), ylab='Recruits', xlab='SSB',
model.frame(FLQuants(rec=propagate(x@rec, dims(x)$iter), ssb=ssb, fitted=fitted)),
col='red', main=main, xlim=c(0, max(ssb, na.rm=TRUE)),
ylim=c(0, max(x@rec, na.rm=TRUE)+(max(x@rec,na.rm=TRUE)/10)),
groups=iter, type='l', scales=scales), split=c(1,1,2,3), more=TRUE)
# Add model line & lowess
trellis.focus("panel", 1, 1)
lpoints(x@ssb, x@rec, col='black', cex=cex)
llines(lowess(x)$ssb, lowess(x)$rec, col='blue', lty=4)
trellis.unfocus()
# residuals or log residuals?
if(log.resid)
resid <- model.frame(FLQuants(resid=log(rec(x) / fitted(x))))
else
resid <- model.frame(FLQuants(resid=residuals(x)))
# 2. Residuals plotted against year
print(xyplot(formula(paste("resid~year", cond)), ylab='Residuals', xlab='',
data=resid, scales=scales, panel=srpanel, groups=iter,
main='Residuals by year'), split=c(2,1,2,3), more=TRUE)
# 3. Residuals at time t vs. residuals at time t+1
print(xyplot(formula(paste("resid1~resid", cond)), ylab='Residuals at t+1',
xlab='Residuals at t', data=cbind(resid, resid1=c(resid$resid[-1], NA)),
panel=respanel, main='AR(1) Residuals', scales=scales), split=c(1,2,2,3), more=TRUE)
# 4. Residuals plotted against SSB
print(xyplot(formula(paste("resid~ssb", cond)), ylab='Residuals', xlab='SSB',
data=cbind(resid, ssb=c(x@ssb)),
panel=srpanel, main='Residuals by SSB', scales=scales), split=c(2,2,2,3), more=TRUE)
# 5. qqplot of residuals
print(qqmath(formula(paste("~resid", cond)), ylab='Residuals',
xlab='Sample Quantiles', data=resid, scales=scales,
panel = function(x, ...) {
panel.qqmath(x, ..., , col='gray40', cex=cex)
panel.qqmathline(x, ..., col='red')
}, main='Normal Q-Q Plot'), split=c(1,3,2,3), more=TRUE)
# 6. Residuals plotted against Recruits
print(xyplot(formula(paste("resid~fitted", cond)), ylab='Residuals',
xlab='Recruits hat', data=cbind(resid, fitted=c(x@fitted)),
panel=srpanel, main='Residuals by Estimated Recruits', scales=scales),
split=c(2,3,2,3), more=FALSE)
invisible()
}
)
# }}}
# lowess {{{
setMethod('lowess', signature(x='FLSR', y='missing', f='ANY', delta='ANY', iter='ANY'),
function(x, f=2/3, iter=3, delta=0.01 * diff(range(ssb(x)[!is.na(ssb(x))])))
{
# output object
rec <- FLQuant(dimnames=dimnames(rec(x))[1:5], iter=dims(x)$iter, units=units(rec(x)))
ssb <- FLQuant(dimnames=dimnames(ssb(x))[1:5], iter=dims(x)$iter, units=units(ssb(x)))
for(i in seq(dims(x)$iter))
{
idx <- array(as.logical(is.na(iter(rec(x), i)) + is.na(iter(ssb(x), i))),
dim=dim(iter(rec(x),i)))
out <- lowess(iter(rec(x),i)@.Data[!idx]~iter(ssb(x),i)@.Data[!idx],
f=f, delta=delta, iter=iter)
iter(rec, i)[!idx][order(ssb(x)[!idx])] <- out$y
iter(ssb, i)[!idx][order(ssb(x)[!idx])] <- out$x
}
return(FLQuants(rec=rec, ssb=ssb))
}
) # }}}
# fmle {{{
setMethod("fmle", signature(object="FLSR", start="ANY"),
function(object, start, ...)
{
res <- callNextMethod()
# AR1 models
if('rho' %in% dimnames(params(object))$params)
{
n <- dim(rec(res))[2]
rho <- c(params(res)['rho',])
residuals(res) <- as.numeric(NA)
residuals(res)[,-1] <- (rec(res)[,-1] - rho*rec(res)[,-n] - fitted(res)[,-1] +
rho*fitted(res)[,-n])
}
# lognormal models
else if(object@logerror)
residuals(res) <- log(rec(res)) - log(fitted(res))
return(res)
}
) # }}}
# ab {{{
setMethod('ab', signature(x='FLSR', model='missing'),
function(x)
{
res <- x
model(res) <- sub('SV', '', SRModelName(model(x)))
params(res) <- ab(params(x), SRModelName(model(x)))
return(res)
}
) # }}}
# sv {{{
setMethod('sv', signature(x='FLSR', model='missing'),
function(x, spr0=params(x)['spr0',])
{
res <- x
model(res) <- SRModelName(model(x))
model(res) <- paste(SRModelName(model(x)), 'SV', sep='')
params(res) <- sv(params(x), SRModelName(model(x)), spr0=spr0)
return(res)
}
) # }}}
# parscale {{{
setMethod('parscale', signature(object='FLSR'),
function(object) {
rec <- rec(object)
ssb <- ssb(object)
res <- switch(SRModelName(model(object)),
bevholt =c(a=mean(rec,na.rm=T), b=mean(ssb,na.rm=T)),
ricker =c(a=mean(rec/ssb,na.rm=T), b=mean(ssb,na.rm=T)),
segreg =c(a=mean(rec/ssb,na.rm=T), b=mean(ssb,na.rm=T)),
shepherd =c(a=mean(rec,na.rm=T), b=mean(ssb,na.rm=T), c=1),
cushing =c(a=mean(rec/ssb,na.rm=T), b=1),
bevholtSV =c(s=1,v=mean(ssb,na.rm=T),spr0=mean(ssb/rec,na.rm=T)),
rickerSV =c(s=1,v=mean(ssb,na.rm=T),spr0=mean(ssb/rec,na.rm=T)),
segregSV =c(s=1,v=mean(ssb,na.rm=T),spr0=mean(ssb/rec,na.rm=T)),
cushingSV =c(s=1,v=mean(ssb,na.rm=T),spr0=mean(ssb/rec,na.rm=T)),
shepherdSV =c(s=1,v=mean(ssb,na.rm=T),c=1,spr0=mean(ssb/rec,na.rm=T)))
if(is.null(res))
stop("SR model not recognized")
else
return(res)
}
) # }}}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.