R/tv.workup.fu.kriging.R

In ices-tools-dev/NephAssess: SG Nephrops assessment functions

# TV workup function for FU where kriging is used
# This function generates estimates based on the standard method
# with an input object (survey.data) containing standard estimates 
# and kriged estimates. The kriged estimates will be appended afterward.


tv.workup.fu.kriging <-
  function(wk.dir, strata.object, survey.data, f.u)
  {
    
    old.dir <- setwd(wk.dir)
    on.exit(setwd(old.dir))
    
    if(missing(survey.data) || missing(f.u) || missing(strata.object)) stop("strata.object, survey.data and f.u must be specified")
    
    f.u <- check.fu(f.u)
    
    inputs <- 
      switch(f.u,
             "clyde" = list( 
               areas = c(MUD = 716.8, "MUDDY SAND" = 664.6, "SANDY MUD" = 698.6)),
             
             "south minch" = list(
               areas = c(MUD = 303.1, "MUDDY SAND" = 2027.6, "SANDY MUD" = 2740.6)),
             
             "moray firth" = list(
               areas = c("MUDDY SAND E" = 646, "MUDDY SAND C" = 698.44, "MUDDY SAND W" = 681.68, "SANDY MUD" = 169)),
             
             "firth forth" = list( 
               areas = c("MUDDY SAND E" = 395.019, "MUDDY SAND C" = 210.596, "MUDDY SAND W" = 138.90, "SANDY MUD" = 170.48)),
             
             "north minch" = list( 
               areas = c(U = 655.6, V = 425, W = 563, X = 131)),
             
             "jura" = list( 
               areas = c(MUD = 90, "MUDDY SAND" = 142, "SANDY MUD" = 150)),
             
             "noup" = list( 
               areas = c("MUDDY SAND" = 409)),
             
             "fladen" = list( 
               areas = c(C = 15634 , MC = 4303.5 , MF = 4966.6, "F" = 3248.1))
      )
    
    if(f.u == "fladen")
    { 
      names(strata.object) [dim(strata.object)[2]-1] <- "strata_type_old"
      names(strata.object) [dim(strata.object)[2]] <- "strata_type"
    }
    
    
    bias <- 
      switch(f.u,
             "clyde" = 1.19,
             "south minch" = 1.32,
             "moray firth" = 1.21,
             "firth forth" = 1.18,
             "north minch" = 1.33,
             "jura" = 1.19,
             "noup"=1.35,
             "fladen" = 1.35)
    
    strata.object <- check.strata(strata.object)
    
    total.area <- sum(inputs$areas)
    
    strata.object <- strata.object[order(strata.object$strata_type),]
    
    stratum.names <- unique(strata.object$strata_type)
    if (!all(stratum.names %in% names(inputs$areas))) stop("stratum names not consistent in inputs")
    
    strata.list <- vector("list", length(stratum.names))
    names(strata.list) <- stratum.names
    
    for (i in stratum.names)
    {
      tmp <- list()
      tmp$strata_type <- strata.object[strata.object$strata_type == i, ]
      av.dens <- tmp$strata_type$average.density
      tmp$area <- unname(inputs$areas [i])
      tmp$no.stations <- length(av.dens)
      tmp$sample.mean <- mean(av.dens)
      tmp$sample.variance <- var(av.dens)
      tmp <- variance.calcs (tmp)
      #tmp $ biomass <- tmp $ total.n * inputs $ mean.wt
      
      strata.list[[i]] <- tmp
    }
    
    ## if any strata have only 1 station fill in the variance by assuming constant CV
    which.fill <- sapply(strata.list, "[[", "no.stations") == 1
    fill.cv <- mean(sapply(strata.list, function(x) x$sample.mean/sqrt(x$sample.variance)), na.rm = TRUE)
    strata.list[which.fill] <- lapply(strata.list[which.fill], function(x){x$sample.variance <- (x$sample.mean/fill.cv)^2; x})
    strata.list[which.fill] <- lapply(strata.list[which.fill], variance.calcs)
    
    
    ## Proportions of the mean
    
    sum.var.mean <- sum(sapply(strata.list, "[[", "mean.var"))
    strata.list <- lapply(strata.list, function(x){x$prop <- x$mean.var/sum.var.mean; x})
    
    
    #################################################
    #####  Final table   - TV Survey results  #######
    #################################################
    
    abundance <- sum(sapply( strata.list, "[[", "total.n"))
    mean.density <- abundance/total.area
    confidence.interval <-  2 * sqrt(sum.var.mean)
    
    
    recent.year <- cbind(rev(survey.data$year)[1] + 1, abundance, mean.density, confidence.interval, NA, NA, NA, NA)
    colnames(recent.year) <- c("year", 
                               "abundance", 
                               "mean.density", 
                               "confidence.interval",
                               "krig.mean",
                               "krig.cv",
                               "krig.abundance",
                               "krig.confidence.interval")
    # Survey results tables
    final.table <- round(rbind(survey.data, recent.year), 3)
    final.table.bias.corrected <- cbind(year = final.table$year, 
                                        abundance = round(final.table$abundance/bias), 
                                        mean.density = round(final.table$mean.density/bias, 3), 
                                        confidence.interval = round(final.table$confidence.interval/bias),
                                        krig.mean = round(final.table$krig.mean/bias, 3),
                                        krig.abundance = round(final.table$krig.abundance/bias),
                                        krig.confidence.interval = round(final.table$krig.confidence.interval/bias)) 
    
    
    # if no fishstats folder create one
    if(!file.exists("fishstats")){dir.create("fishstats")}
    
    write.table(final.table, paste(getwd(), "/", "fishstats/", f.u, "_TV_results.csv", sep = ""), row.names = FALSE, sep = ",") 
    write.table(final.table.bias.corrected, paste(getwd(), "/", "fishstats/", f.u, "_TV_results_bias_corrected.csv", sep = ""), row.names = FALSE, sep = ",") 
    
    
    ####################
    ######  PLOT #######
    ####################
    
    # Not corrected for bias
    
    png(paste(getwd(), "/", "fishstats/", f.u, "_TV_results.png", sep = ""), width=2200, height=1500, pointsize=50)
    par(las = 1, bty = "l", ann = FALSE, pch = 16, mfrow = c(1,1))
    
    with( final.table,
          {
            plot(year, abundance, cex = 0.8, type = "o", lwd = 3,
                 ylim = c(0, max(abundance + confidence.interval, na.rm = TRUE)))
            arrows(year, abundance - confidence.interval, year, abundance + confidence.interval, angle = 90, code = 3, length = 0.2,lwd=3)
          })
    title(main = f.u, ylab = "Nephrops abundance (millions)", xlab = "Year")
    
    dev.off()
    
    
    # Corrected for bias
    
    png(paste(getwd(), "/", "fishstats/", f.u, "_TV_results_bias_corrected.png", sep = ""), width=2200, height=1500, pointsize=50)
    par(las = 1, bty = "l", ann = FALSE, pch = 16, mfrow = c(1,1))
    
    with( final.table,
          {
            plot(year, abundance/bias, cex = 0.8, type = "o", lwd = 3,
                 ylim = c(0, max(abundance/bias + confidence.interval/bias, na.rm = TRUE)))
            arrows(year, abundance/bias - confidence.interval/bias, year, abundance/bias + confidence.interval/bias, angle = 90, code = 3, length = 0.2,lwd=3)
          })
    title(main = f.u, ylab = "Nephrops abundance (millions)", xlab = "Year")
    
    dev.off()
    
    #########################################################
    ##### Final table - Results by stratum last 3 years #####
    #########################################################
    
    results.by.stratum <-
      data.frame(
        stratum.names,
        sapply(strata.list, "[[", "area" ),
        sapply(strata.list, "[[", "no.stations"),
        round( sapply(strata.list, "[[", "sample.mean"), 3),
        round( sapply(strata.list, "[[", "sample.variance"), 3),
        round( sapply(strata.list, "[[", "total.n"), 1),
        round( sapply(strata.list, "[[", "mean.var"), 0),
        round( sapply(strata.list, "[[", "prop"), 3),
        stringsAsFactors = FALSE
      )
    
    # calculate totals for each column and bind to previous table
    sums <- c( "Total", lapply( results.by.stratum[-1], sum ) )
    sums[4:5] <- NA
    names(sums) <- names(results.by.stratum)
    results.by.stratum <- rbind(results.by.stratum, sums)
    
    names(results.by.stratum) <- c("Stratum", "Area (km2)", "Number of stations", "Mean burrow density (no./m2)", "Observed variance", "Abundance (millions)",
                                   "Stratum variance", "Proportion of total variance")
    # convert to character...
    results.by.stratum[] <- lapply( results.by.stratum, as.character )
    results.by.stratum [ is.na(results.by.stratum) ] <- ""
    
    #Correct for bias
    results.by.stratum.bias.corrected<- results.by.stratum
    results.by.stratum.bias.corrected[,"Abundance (millions)"]<- round(as.numeric(results.by.stratum.bias.corrected[,"Abundance (millions)"])/bias, 1)
    results.by.stratum.bias.corrected[,"Mean burrow density (no./m2)"]<- round(as.numeric(results.by.stratum.bias.corrected[,"Mean burrow density (no./m2)"])/bias, 3)
    results.by.stratum.bias.corrected[,"Observed variance"]<- round(as.numeric(results.by.stratum.bias.corrected[,"Observed variance"])/bias^2, 3)
    results.by.stratum.bias.corrected[,"Stratum variance"]<- round(as.numeric(results.by.stratum.bias.corrected[,"Stratum variance"])/bias^2, 0)
    results.by.stratum.bias.corrected [ is.na(results.by.stratum.bias.corrected) ] <- ""
    
    # Write table to a csv
    write.table(results.by.stratum, paste(getwd(), "/", "fishstats/", f.u, "_TV results by stratum.csv", sep = ""), sep = ",", row.names=FALSE)
    write.table(results.by.stratum.bias.corrected, paste(getwd(), "/", "fishstats/", f.u, "_TV results by stratum_bias_corrected.csv", sep = ""), sep = ",", row.names=FALSE)
    
    strata_list <- lapply(strata.list, function(x) {class(x) <- "tvworkup"; x})
    class(strata.list) <- "tvworkup"
    
    #########################################################
    ##### CVs - Standard relative error calculation #####
    #########################################################
    n.station<- as.numeric(results.by.stratum.bias.corrected[,"Number of stations"])
    var.station<- as.numeric(results.by.stratum.bias.corrected[,"Observed variance"])
    mean.station<- as.numeric(results.by.stratum.bias.corrected[,"Mean burrow density (no./m2)"])
    
    n.station<- n.station[1:(length(n.station)-1)]
    var.station<- var.station[1:(length(var.station)-1)]
    mean.station<- mean.station[1:(length(mean.station)-1)]
    
    var.total<- sum(var.station*n.station)/sum(n.station)
    mean.total<- sum(mean.station*n.station)/sum(n.station)
    
    sd.total<- sqrt(var.total)/sqrt(sum(n.station))
    cv.total<- round(sd.total/mean.total,3)
    names(cv.total)<- f.u
    #Write CV
    write.table(as.matrix(cv.total), paste(getwd(), "/", "fishstats/", f.u, "_CV.csv", sep = ""), sep = ",", row.names=TRUE, col.names=FALSE)
    
    invisible( strata.list )
  }