R/estN.R

In scape: Statistical Catch-at-Age Plotting Environment

estN.int <- function(P, Phat)
{
  numerator <- sapply(rownames(P), function(t) sum(Phat[t,]*(1-Phat[t,])))
  denominator <- sapply(rownames(P), function(t) sum((P[t,]-Phat[t,])^2))
  nhat <- as.numeric(numerator / denominator)

  nhat
}



estN <- function(model, what="CAc", series=NULL, init=NULL, FUN=mean,
                 ceiling=Inf, digits=0)
{
  ## 1  Parse args
  what <- match.arg(what, c("CAc","CAs","CLc","CLs"))
  x <- if(class(model)=="scape") model[[what]] else model  # allow data frame
  if(is.null(x))
    stop("element '", what, "' not found")
  x$Column <- if(substring(what,1,2)=="CA") x$Age else x$Length

  ## 2  Extract series
  if(is.null(series))
    series <- unique(x$Series)
  if(length(series) > 1)
  {
    output <- lapply(series, function(s)
                     estN(model=model, what=what, series=s, init=init, FUN=FUN,
                          ceiling=ceiling, digits=digits))
    names(output) <- series
  }
  else
  {
    ok.series <- x$Series %in% series
    if(!any(ok.series)) stop("please check if the 'series' argument is correct")
    x <- x[!is.na(x$Obs) & ok.series,]

    ## 3  Calculate nhat
    P <- tapply(x$Obs, list(x$Year,x$Column), mean)     # pool sexes
    Phat <- tapply(x$Fit, list(x$Year,x$Column), mean)  # pool sexes
    nhat <- estN.int(P, Phat)

    ## 4  Prepare init
    years <- unique(x$Year)
    if(is.null(init))
      init <- tapply(x$SS, x$Year, unique)
    else if(class(init) == "scape")
      init <- tapply(init[[what]]$SS, init[[what]]$Year, unique)
    else if(identical(init,FALSE))
      init <- nhat
    if(length(init)!=1 && length(init)!=length(years))
      stop("'init' should be NULL, model, vector, FALSE, or 1")

    ## 5  Prepare output
    output <- as.numeric(FUN(nhat) * init / FUN(init))
    output <- pmin(output, ceiling)
    if(!is.null(digits))
      output <- round(output, digits=digits)
    if(length(output) == length(years))
      names(output) <- years
  }

  output
}