Skip to contents

04 Create IC format from estimation results

Source: vignettes/v04-Create-IC-format-from-estimation-results.Rmd
v04-Create-IC-format-from-estimation-results.Rmd

Introduction

This document shows how to take the output from doEstimationForAllStrata() and format it into the InterCatch format. This is done by using the function exportEstimationResultsToInterCatchFormat() in the RDBEScore package.

Load the package

Step 1) load and prepare some test data 

# load an example dataset

# H1 directory
dirH1 <- "../tests/testthat/h1_v_20250211/"

fixProbs <- function(myH1RawObject){

  # Edit our data so that we have SRSWOR on each level and calculate the probs
  myH1RawObject[["VS"]]$VSselectMeth <- "SRSWOR"
  myH1RawObject[["VS"]]$VSincProb <- myH1RawObject[["VS"]]$VSnumSamp / myH1RawObject[["VS"]]$VSnumTotal
  myH1RawObject[["VS"]]$VSselProb <- 1/myH1RawObject[["VS"]]$VSnumTotal
  myH1RawObject[["FT"]]$FTselectMeth <- "SRSWOR"
  myH1RawObject[["FT"]]$FTincProb <- myH1RawObject[["FT"]]$FTnumSamp / myH1RawObject[["FT"]]$FTnumTotal
  myH1RawObject[["FT"]]$FTselProb <- 1/myH1RawObject[["FT"]]$FTnumTotal
  myH1RawObject[["FO"]]$FOselectMeth <- "SRSWOR"
  myH1RawObject[["FO"]]$FOincProb <- myH1RawObject[["FO"]]$FOnumSamp / myH1RawObject[["FO"]]$FOnumTotal
  myH1RawObject[["FO"]]$FOselProb <- 1/myH1RawObject[["FO"]]$FOnumTotal
  myH1RawObject[["SS"]]$SSselectMeth <- "SRSWOR"
  myH1RawObject[["SS"]]$SSincProb <- myH1RawObject[["SS"]]$SSnumSamp / myH1RawObject[["SS"]]$SSnumTotal
  myH1RawObject[["SS"]]$SSselProb <- 1/myH1RawObject[["SS"]]$SSnumTotal
  myH1RawObject[["SA"]]$SAselectMeth <- "SRSWOR"
  myH1RawObject[["SA"]]$SAincProb <- myH1RawObject[["SA"]]$SAnumSamp / myH1RawObject[["SA"]]$SAnumTotal
  myH1RawObject[["SA"]]$SAselProb <- 1/myH1RawObject[["SA"]]$SAnumTotal

  # Update our test data with some random sample measurements (it didn't include these)
  # set the random seed
  set.seed(1234)
  myH1RawObject[['SA']]$SAsampWtLive <- round(runif(n=nrow(myH1RawObject[['SA']]),min = 1, max = 100))
  myH1RawObject[['SA']]$SAsampWtMes <- round(runif(n=nrow(myH1RawObject[['SA']]),min = 1, max = 100))

  myH1RawObject

}


## Step 1) load and prepare some test data
myH1RawObject <- createRDBESDataObject(input  = dirH1)
#Filter our data for WGRDBES-EST TEST 1, 1965, H1
myCatchFraction <- "Lan"
mySpecies <- 1019159

myFields <- c("DEyear","DEhierarchy","DEsampScheme","DEstratumName","SAspeCode","SScatchFra")
myValues <- c(1965,1,"National Routine","DE_stratum1_H1",mySpecies,myCatchFraction)
myH1RawObject <- filterRDBESDataObject(myH1RawObject,
                                       fieldsToFilter = myFields,
                                       valuesToFilter = myValues )
myH1RawObject <- findAndKillOrphans(myH1RawObject)
myH1 <- fixProbs(myH1RawObject)

Step 2) Estimate for live weight for the data 

myTestData <- createRDBESEstObject(myH1, 1, stopTable = "SA")
# Get rid of rows that don't have an SA row
myTestData <- myTestData[!is.na(myTestData$SAid),]
# Estimate using the data
myStrataResults <- doEstimationForAllStrata(myTestData, "SAsampWtLive")

Step 3) Export the results to InterCatch format (just HI and SI for this data) 


# Get our estimated values for the PSU
psuEstimates <- myStrataResults[myStrataResults$recType == "VS",]

# This is the data we will export to IC format
dataToOutput <- data.frame(Country = myH1[["SD"]]$SDctry,
                           Year = myH1[["DE"]]$DEyear,
                           SeasonType = NA,
                           Season = NA,
                           Fleet= NA,
                           AreaType = "Stratum",
                           FishingArea = psuEstimates$stratumName,
                           DepthRange = "NA",
                           Species = mySpecies ,
                           Stock = mySpecies,
                           CatchCategory = substr(myCatchFraction, 1, 1),
                           ReportingCategory = "A",
                           Usage = "H",
                           SamplesOrigin = "O",
                           UnitCATON = "kg",
                           CATON = psuEstimates$est.total/1000.0,
                           varCATON = psuEstimates$var.total/1000000.0,
                           Sex = "NA",
                           PlusGroup = "NA",
                           MeanWeight = -9,
                           unitMeanWeight = "kg",
                           unitCANUM = "n",
                           UnitAgeOrLength = "year",
                           NumberCaught = -9,
                           MeanWeight = -9,
                           CANUMtype = "Age",
                           AgeLength = c(1,1,1,2,2,2,3,3,3),
                           NumSamplesAge = c(300,300,300,200,200,200,100,100,100))

tempIC <- exportEstimationResultsToInterCatchFormat(dataToOutput)
print(tempIC)
#>  [1] "HI,ZW,1965,NA,NA,NA,Stratum,VS_stratum1,NA,NA,-9,NA"                                                                         
#>  [2] "SI,ZW,1965,NA,NA,NA,Stratum,VS_stratum1,NA,1019159,1019159,L,A,NA,H,O,NA,kg, 4835.556,-9,3006072,NA,NA,NA"                   
#>  [3] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum1,NA,1019159,1019159,L,A,NA,Age,1,NA,-9,-9,-9,300,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"
#>  [4] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum1,NA,1019159,1019159,L,A,NA,Age,2,NA,-9,-9,-9,200,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"
#>  [5] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum1,NA,1019159,1019159,L,A,NA,Age,3,NA,-9,-9,-9,100,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"
#>  [6] "HI,ZW,1965,NA,NA,NA,Stratum,VS_stratum2,NA,NA,-9,NA"                                                                         
#>  [7] "SI,ZW,1965,NA,NA,NA,Stratum,VS_stratum2,NA,1019159,1019159,L,A,NA,H,O,NA,kg,15217.778,-9,8936549,NA,NA,NA"                   
#>  [8] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum2,NA,1019159,1019159,L,A,NA,Age,1,NA,-9,-9,-9,300,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"
#>  [9] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum2,NA,1019159,1019159,L,A,NA,Age,2,NA,-9,-9,-9,200,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"
#> [10] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum2,NA,1019159,1019159,L,A,NA,Age,3,NA,-9,-9,-9,100,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"
#> [11] "HI,ZW,1965,NA,NA,NA,Stratum,VS_stratum3,NA,NA,-9,NA"                                                                         
#> [12] "SI,ZW,1965,NA,NA,NA,Stratum,VS_stratum3,NA,1019159,1019159,L,A,NA,H,O,NA,kg,10311.111,-9,8937750,NA,NA,NA"                   
#> [13] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum3,NA,1019159,1019159,L,A,NA,Age,1,NA,-9,-9,-9,300,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"
#> [14] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum3,NA,1019159,1019159,L,A,NA,Age,2,NA,-9,-9,-9,200,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"
#> [15] "SD,ZW,1965,NA,NA,NA,Stratum,VS_stratum3,NA,1019159,1019159,L,A,NA,Age,3,NA,-9,-9,-9,100,-9,kg,n,year,NA,NA,-9,-9,-9,-9,-9,-9"

See also other package vignettes:

#END