In MortenVinther/FishStomachs: Compilation of stomach contents data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  options(rmarkdown.html_vignette.check_title = FALSE),
   fig.width=7, fig.height=8
)

FishStomachs provides functions for estimating the mean and variance of the stomach contents from bootstrapping. First, data from Example 03 are loaded.

library(FishStomachs)
suppressMessages(library(dplyr))
suppressMessages(library(readr))

load(file=file.path(system.file( package = "FishStomachs"),"extdata","ex03half.Rdata"),verbose=TRUE)

Options for bootstrapping are provided by the bootstrapping control (see below). It is a list, where boots_strata defines the strata for bootstrapping, in this case, data are made for each combination of year, quarter and predator. Within such stratum, bootstrap samples are drawn from the pool of boots_id, in this case defined as individual trawl hauls. The variable predator_size is not used for boots_strata to reflect that several predator lengths are often sampled at the same haul. If you want to draw replicates from the pool of samples (haul, predator and predator size) for each predator size, the boots_strata should include the predator size as well (boots_strata = expression(paste(year, quarter, pred_name, pred_size, sep='-')).

############### add variables for bootstrapping
s<-edit_control(s, bootstrapping = list(
        boots_id = expression(paste(ship,rectangle,year,quarter,station,haul,sep='_')),
        boots_strata = expression(paste(year,quarter,pred_name,sep='-'))
        ))

print(get_control(s),'bootstrapping')
s<-bootstrap_addVar(s)
names(s[['PRED']])

The call to bootstrap_addVar() adds the variables for bootstrapping to data.

The resulting stratification with number of strata and number of samples is shown below.

bootstrap_show(s,show=c("strata",'sample')[1],vari=c("stomach","sample")[1])
bootstrap_show(s,show=c("strata",'sample')[1],vari=c("stomach","sample")[2])

The tables above show that e.g. cod in quarter 1 has 18 "hauls" with in total 42 stomachs for predator size 0100-0149. The table below shows that each "haul" includes a number of stomachs from various predator sizes.

head(bootstrap_show(s,show=c("strata",'sample')[2],vari=c("stomach","sample")[1])[[1]])

We are now ready to do the bootstrapping, however, to speed up, tasks common for all replicates are done first. That is the definition of the included prey species and the specification of stratifications for the calculation of population diet.

# a very complex way of getting a list of prey names!
keep_prey <-(read_csv(file.path(system.file( package = "FishStomachs"),"extdata","species_info.csv"),
                      col_types = cols())  %>% 
               dplyr::filter(prey_sp)  %>% select(species) %>% unique())$species
keep_prey #input is just a vector of prey names

s<-edit_control(s,sel_preys=keep_prey)

from_to <-make_from_to_species(inp_file=
        file.path(system.file( package = "FishStomachs"),"extdata","from_to_species.csv"))

s<-edit_control(s,
                calc_sub_strata=list(
                  # transform into relative weight before data are compiled
                  relative_weight=FALSE,      
                  # use number of stomach in the sample as weighting factor
                  weighting_factor=expression(n_tot),
                  # set file to NA when weighting_factor is given
                  weighting_factor_file=NA
                ),
                calc_strata=list(
                  relative_weight=FALSE,
                  weighting_factor=expression(sqrt(mean_cpue)),
                  weighting_factor_file=NA
                ),
                calc_total=list(
                  relative_weight=FALSE,
                  # set to NA when input is provided in a file
                  weighting_factor=NA,
                  # use data from a specified file (in those case, a file provided by the FishStomachs)
                  weighting_factor_file=file.path(system.file( package = "FishStomachs"),"extdata","weighting_total.csv")
                )
)

The processing of observations into population diet is put into a function, do_boots, which is user defined. In this case, the code included in do_boots is copied from Example 03.

do_boots<-function(s) {

# cat('*') # to show that something is happening  

# to keep information on bootstrap replicate  
rep_id<-unlist(s[['PRED']][1,'rep_id'])

# by sample
s<-redist_unidentified_prey_sp(s,dist_time=stratum_time,dist_area=sample_id,
        dist_pred_size=pred_size, do_only=c(1,2),from_to_species=from_to,
        by_prey_size=FALSE,remains_to_other = FALSE)

# by sub_area (=rectangel)
s<-redist_unidentified_prey_sp(s,dist_time=stratum_time,dist_area=stratum_sub_area,
        dist_pred_size=pred_size, do_only=c(1,2),from_to_species=from_to,
        by_prey_size=FALSE,remains_to_other = FALSE)

# by area (=roundfish area)
s<-redist_unidentified_prey_sp(s,dist_time=stratum_time,dist_area=stratum_area,
        dist_pred_size=pred_size, do_only=c(1,2,3),from_to_species=from_to,
        by_prey_size=FALSE,remains_to_other = FALSE)

# all areas
s<-redist_unidentified_prey_sp(s,dist_time=stratum_time,dist_area='All',
        dist_pred_size=pred_size, do_only=c(1,2,3),from_to_species=from_to,
        by_prey_size=FALSE,remains_to_other = FALSE)

# all areas and half-year (dist_time is changed)
s<-redist_unidentified_prey_sp(s,dist_time=paste(year,ifelse(quarter %in% c(1,2),'S1','S2')),
        dist_area='All',dist_pred_size=pred_size, do_only=c(1,2,3),from_to_species=from_to,
        by_prey_size=FALSE,remains_to_other = FALSE)

# the same as above, but unallocated remains to "other food".
s<-redist_unidentified_prey_sp(s,dist_time=paste(year,ifelse(quarter %in% c(1,2),'S1','S2')),
        dist_area='All',dist_pred_size=pred_size, do_only=c(1,2,3),from_to_species=from_to,
        by_prey_size=FALSE,remains_to_other = TRUE)

# by statum_area (=roundfish area)
s<-redist_unidentified_prey_lengths(s,dist_time=stratum_time,dist_area=stratum_area,
      dist_pred_size=pred_size,
      remains_to_other = FALSE, # do not add records with missing length allocation to "other food"
      others_to_other=TRUE      # do add other species that selected prey species to "other food"
      )

# within all areas
s<-redist_unidentified_prey_lengths(s,dist_time=stratum_time,dist_area='All',
      dist_pred_size=pred_size,remains_to_other = FALSE)

#  within all areas, and year (ignoring quarter)
s<-redist_unidentified_prey_lengths(s=s,dist_time=substr(stratum_time,1,4),dist_area='All',
      dist_pred_size=pred_size,remains_to_other = FALSE)

# as above, but with the conversion of unallocated into "other"
s<-redist_unidentified_prey_lengths(s=s,dist_time=substr(stratum_time,1,4),dist_area='All',
      dist_pred_size=pred_size,remains_to_other = TRUE)

d<-calc_population_stom(s)
d[['PRED']]$rep_id<-rep_id
return(d)

}

Data for bootstrapping is made by the function bootstrap_data(), which creates a new STOMobs object including the bootstrap replicate data. The result of the call, bd, is a list of STOMobs objects.

make_dataSet<-FALSE  # just a varibale to create data set for the package 

if (make_dataSet) n_boot<-1:500 else n_boot<-1:3  # make 3 bootstrap replicates

bd<-lapply(n_boot,function(x) bootstrap_data(s,seed=x,rep_id=x))

print(bd[[1]],show_attributes=FALSE)

print(bd[[2]],show_attributes=FALSE)

The calculation of the population diet is done for each STOMobs object as shown below. The result is a list of STOMdiet objects.

#convert observations into population diet
bt<-lapply(bd,do_boots)

bt[[1]] # STOMdiet object 

bt[[2]] # STOMdiet object

The STOMdiet can be changed. In this case, predator and prey names are changed into English names and prey weights into relative weights

bt2<-lapply(bt,function(x){
 x<-from_to_species_diet(x,pred_from_to=c("species","short"),prey_from_to=c("species","short"),
      sp_info_file=file.path(system.file( package = "FishStomachs"),"extdata",'species_info.csv'))
 x<-diet_relative(x)
 return(x)

})

# change to English names in the source data for bootstrap, for a later merge
 source<-from_to_species_diet(s,pred_from_to=c("species","short"),prey_from_to=c("species","short"),
    sp_info_file=file.path(system.file( package = "FishStomachs"),"extdata",'species_info.csv'))

The estimated diets for each bootstrap replicates are different as shown below, where the two first replicates are compared. The size of the bubbles indicates the difference, d1 minus d2, between the two diets. Red is negative, green is positive and blue shows that data are missing in one of the two sets.

# comparison of two replicates
plotdif(d1=bt2[[1]],d2=bt2[[2]],show_plot=TRUE,cut_pred_size=c(1,4),cut_prey_size=c(2,4),addTitle=TRUE,
        tAngle=90, relative=FALSE,maxDif=3,
        byVar=c('year-quarter','year','quarter','none')[4])

Missing prey and prey size combination in the bootstrap replicates can be filled in with a selected value, in this case with zero.

bt3<-add_missing_boots(bt2, mis_value=0)

plotdif(d1=bt3[[1]],d2=bt3[[2]],show_plot=TRUE,cut_pred_size=c(1,4),cut_prey_size=c(2,4),addTitle=TRUE,
        tAngle=90, relative=FALSE, maxDif=3,
        byVar=c('year-quarter','year','quarter','none')[4])

Plots of the bootstrapped diet are available. In the case below, a subset of data the replicates are first made, and the distribution of the individual replicates is shown, by prey species and their size.

# sub set for plot
bt4<-lapply(bt3,function(x){
  subset(x,pred_name=='Cod' &  pred_size=='0400-0499' & stratum_time=='1991-Q1')
})

plotboots.size(b=bt4,show_plot=TRUE,cut_pred_size=c(1,10),cut_prey_size=c(2,4),addTitle=TRUE,tAngle=90)

In another plot, prey sizes are combined, and the plot shows the distribution of prey weight by prey species and predator size.

# sub set for plot
bt4<-lapply(bt3,function(x){
  subset(x,pred_name=='Cod' & stratum_time=='1991-Q1' & !(pred_size %in% c('0100-0149','0150-0199')))
})


plotboots(b=bt4,show_plot=TRUE,cut_pred_size=c(1,4),addTitle=TRUE,tAngle=90)

The bootstrap mean and variance from replicates are calculated by a call to bootsMean(). First, the raw observations are loaded.

# load diet data (point estimate) from example 03
load(file=file.path(system.file( package = "FishStomachs"),"extdata","ex03d.Rdata"),verbose=TRUE)

# change from latin to English names (used later on) 
 pointEst<-from_to_species_diet(d,pred_from_to=c("species","short"),prey_from_to=c("species","short"),
    sp_info_file=file.path(system.file( package = "FishStomachs"),"extdata",'species_info.csv'))

With more replicates, loaded from the package data set, simple mean and variance, and Dirichlet parameters are calculated. See help(bootsMean) for more information.

if (make_dataSet) {
  # sub set for plot
  bt5<-lapply(bt4,function(x){
    subset(x,pred_name=='Cod' &  pred_size=='0400-0499' & stratum_time=='1991-Q1')
  })
  save(bt5, file=file.path(system.file( package = "FishStomachs"),"extdata","bt5.Rdata")) 
}

load(file=file.path(system.file( package = "FishStomachs"),"extdata","bt5.Rdata"),verbose=TRUE)

bb<-bootsMean(b=bt5,pointEst,by_prey_size=FALSE,do_Diri=TRUE)

# summary statistics for predator cod, size 0400-0499, 1991, Q1,  
select(bb,prey_name, phi,param,mu,mean_w,sd_w,prey_w)

plotboots(b=bt5,cut_pred_size=c(1,4),addTitle=TRUE,tAngle=90,statistics=bb)