Skip to contents

Introduction

The secpRod package provides a number of methods to calculate secondary production of populations, both those with clear cohort structure and those where cohort structure is not possible to discern. First, we walk through the available methods using simulated populations that mimic the data structure of real sampling regimes. The parameters to create these simulated populations is outlined in companion article, vignette(“Simulating the sampling of populations”). These simulated populations skip some of the steps of moving from initial sampling processes (e.g., length to mass conversions, sample subsetting, etc.) that are often part of the process. More complete examples that showcase helper functions within secpRod that deal with these processes can be found below in More complete examples. Here, we just showcase the secondary production methods available.

The methods within the package use two main objects, sampleInfo the sample-level information (e.g., site, date , replicate) and density and length or mass. An additional dataframe, taxaInfo that houses species-specific information such as length to mass conversions, and method specific information such as cohort production interval estimates (CPIs), production:biomass ratios (PBs), and growth equations. The columns contain information for the calculation of production. These include, but are not limited to:

  • taxonID: a character string that matches the name of taxonID from sampleInfo

  • massForm: a character string that is coercible to a formula for the conversion from length to mass (e.g., afdm_mg~(a*lengthClass^b))

  • a: numeric variable for the coefficients used in massForm

  • b: numeric variable for the coefficients used in massForm

  • percAsh: numeric integer 0–100

  • method: character string for the method to use. Must be one of the following:

    • ‘is’: increment-summation method. This cohort-based method calculates interval and annual production, PP, from field data as the sum of all interval production, i=1tN¯ΔW\sum_{i = 1}^{t} \overline{N} \Delta W, plus initial biomass, BiniitalB_{iniital}:

    P=Binitial+i=1tN¯ΔWP = B_{initial} + \sum_{i = 1}^{t} \overline{N} \Delta W

    • ‘rs’: removal-summation method is similar to the increment-summation method but instead calculates the prodcution lost during a sampling interval as the product of the decrease in density, ΔN\Delta N and the mean individual mass, W¯\overline{W} plus the increase in biomass, ΔB\Delta B.

    • ‘sf’: size frequency method

  • growthForm: a character string that is coercible to a formula. This formula is parsed and will reference existing information (e.g., size information) if variables are included in the formula. Further, an additional data

  • min.cpi: integer of the minimum estimated cohort production interval for adjusting annual production estimates using the size-frequency method

  • max.cpi: integer of the maximum estimated cohort production interval for adjusting annual production estimates using the size-frequency method

  • pb: numeric of the production to biomass (PB) ratio for the specific taxa. This can take three forms 1) a single value, 2) a vector of values the same length as bootNum, and 3) a string of a distribution to randomly sample, not including the n = (e.g., ‘rnorm(mean = 5, sd = 0.5)’, ‘runif(min = 3, max = 8)’). The function will automatically sample the bootNum values. Fair warning: We suggest you explicitly name the parameters in the function call string. Possible unwanted and unknown things may happen otherwise. The tests for this feature are still under development, it will not reject any nonsensical values (e.g, negative, NA, or Inf). If your distribution samples negative values you will get negative productions estimates–at least until the min.growth feature is implemented.

  • min.growth: a minimal value of production to assign when density is >0 and production estimates for a taxon are <0.

  • notes: notes for researcher use. This column will be maintained in output summaries.

  • …: additional information or data to pass to the function. This is mainly used to pass environmental data to inform growth rate information by included information such as temperature, food availability, etc.

taxaInfo <- data.frame(
  taxonID = c("sppX"),
  massForm = c("afdm_mg~(a*lengthClass^b)*percAsh"),
  # a is NA for this example because we simulated mass growth and don't need to convert from length to mass
  a = c(NA),
  # b is NA for this example because we simulated mass growth and don't need to convert from length to mass
  b = c(NA),
  # percAsh is NA for this example because we simulated mass growth and don't need to convert from length to mass
  percAsh = c(NA),
  # method can accept one or more values. This allow comparisons among different methods.
  method = c("is"),
  g.a = c(NA),
  growthForm = c("log(g_d) = 1 - 0.25*log(afdm_mg)"),
  min.cpi = c(335),
  max.cpi = c(365),
  pb = c("runif(min = 3, max = 8)"),
  min.growth = c(0.001),
  notes = c("This is here for your use. No information will be used, but this column will be maintained in some summaries. See *** for more information.")
)

A multi-method comparison with simulated data

First, we load and visualize the simulated cohort data set, which contains a single species and follows the progression of density and distribution of individual masses through it life cycle.

The data set is available within the secpRod package with:

data('singleCohortSim', package = 'secpRod')

The sampleInfo for this species includes:

#> # A tibble: 10 × 5
#>    taxonID dateID repID density afdm_mg    
#>    <chr>    <dbl> <int>   <dbl> <list>     
#>  1 sppX         1     1     463 <dbl [463]>
#>  2 sppX         1     2     458 <dbl [458]>
#>  3 sppX         1     3     479 <dbl [479]>
#>  4 sppX         1     4     422 <dbl [422]>
#>  5 sppX         1     5     639 <dbl [639]>
#>  6 sppX         1     6     575 <dbl [575]>
#>  7 sppX         1     7     518 <dbl [518]>
#>  8 sppX         1     8     677 <dbl [677]>
#>  9 sppX         1     9     566 <dbl [566]>
#> 10 sppX         1    10     468 <dbl [468]>

summary_stats <- singleCohortSim %>%
  unnest(afdm_mg, keep_empty = TRUE) %>%
  group_by(dateID) %>%
  dplyr::summarise(
    massMean = mean(afdm_mg, na.rm = TRUE),
    massSD= sd(afdm_mg, na.rm = TRUE),
    larvalDensityMean = mean(density, na.rm = TRUE)
  )

sim_plot = 
ggplot(summary_stats, aes(x = dateID)) +
  stat_halfeye(data = splitCohortSim, aes(x = dateID, y = density),
               color = 'green')+
  stat_halfeye(data = splitCohortSim %>% unnest(afdm_mg), aes(x = dateID, y = afdm_mg*100),
               color = 'red')+
  geom_path(aes(y = larvalDensityMean), color = 'green') +
  geom_path(aes(y = massMean * 100), color = 'red') +
  scale_y_continuous(
    name = "Larval Density",
    sec.axis = sec_axis(~./100, name = "Mean Mass (mg)"),
    
  ) +
  theme_minimal() +
  labs(title = "Larval Density and Mean Mass over Time", x = "Day")+
  theme(axis.title.y.right = element_text(color = 'red'))

More complete examples (#more-complete-examples)

Currently, the analysis is fully implemented for the size-frequency and production:biomass ratio methods. This quick tutorial will outline steps to calculate secondary production for whole communities. First, load the package with library(secpRod).

The package comes with a number of data sets of macroinvertebrate species and community data:

  1. simulated data set of a single univoltine species accessed with data("univoltine")

  2. a full community data set from Junker and Cross (2014) that can be accessed with data("wbtData")

The simulated data set is a single data frame of artificially sampled size-frequency data through time.

data("univoltine", package = 'secpRod')

#clean this data set to align 
head(univoltine)
#> # A tibble: 6 × 5
#>   taxonID dateID     repID lengthClass  n_m2
#>   <chr>   <date>     <chr>       <dbl> <dbl>
#> 1 sppA    2023-01-31 1             0.5     0
#> 2 sppA    2023-01-31 1             1       0
#> 3 sppA    2023-01-31 1             1.5     0
#> 4 sppA    2023-01-31 1             2       0
#> 5 sppA    2023-01-31 1             2.5     0
#> 6 sppA    2023-01-31 1             3       0

For this simple example, we can create the required ‘taxaInfo’ object:

# clean simulated data to workable form

taxaInfo <- data.frame(
  taxonID = c("sppA"),
  massForm = c("afdm_mg~(a*lengthClass^b)*percAsh"),
  a = c(0.0025),
  b = c(2.692),
  percAsh = c(0.958),
  method = c("sf"),
  g.a = c(NA),
  growthForm = c("log(g_d) = 1 - 0.25*log(afdm_mg) - "),
  min.cpi = c(335),
  max.cpi = c(365),
  pb = c("runif(min = 3, max = 8)"),
  min.growth = c(0.001),
  notes = c("This is here for your use. No information will be used, but this column will be maintained in some summaries. See *** for more information.")
)

This object represents the taxonomic information for calculating production. Each row contains information for a single species/taxon (in this case the only row because species == 1 ). The columns contain information for calculation of production. These include, but are not limited to:

  • taxonID: a character string that matches the name of taxonID from sampleInfo

  • massForm: a character string that is coercible to a formula for the conversion from length to mass (e.g., afdm_mg~(a*lengthClass^b))

  • a: numeric variable for the coefficients used in massForm

  • b: numeric variable for the coefficients used in massForm

  • percAsh: numeric integer 0–100

  • method: character string for the method to use. Must be one of the following: ‘sf’ (this will be updated as methods are finalized.)

  • growthForm: description forthcoming…

  • min.cpi: integer of the minimum estimated cohort production interval for adjusting annual production estimates using the size-frequency method

  • max.cpi: integer of the maximum estimated cohort production interval for adjusting annual production estimates using the size-frequency method

  • pb: numeric of the production to biomass (PB) ratio for the specific taxa. This can take three forms 1) a single value, 2) a vector of values the same length as bootNum, and 3) a string of a distribution to randomly sample, not including the n = (e.g., ‘rnorm(mean = 5, sd = 0.5)’, ‘runif(min = 3, max = 8)’). The function will automatically sample the bootNum values. Fair warning: We suggest you explicitly name the parameters in the function call string. Possible unwanted and unknown things may happen otherwise. The tests for this feature are still under development, it will not reject any nonsensical values (e.g, negative, NA, or Inf). If your distribution samples negative values you will get negative productions estimates–at least until the min.growth feature is implemented.

  • min.growth: description forthcoming…

  • notes: notes for researcher use. This column will be maintained in output summaries.

This ‘wbtData’ object is a list with two (2) elements:

  1. sampleInfo: a list with 32 data frames. One (1) for each taxonomic entity

  2. taxaInfo: a data frame with 32 rows and 13 columns. Each row represents the taxonomic information for a single taxonomic entity and the columns contain information for calculation of production as described above.

Examples

Single-species walkthrough

A quick walkthrough of the calculation of secondary production for a single species.

## isolate a single species data frame from sampleInfo object

acentrella <- wbtData[['sampleInfo']][['Acentrella sp.']]

## let's take a look at the data set
head(acentrella, 10)
#> # A tibble: 10 × 5
#>    taxonID        repID dateID              lengthClass  n_m2
#>    <chr>          <dbl> <dttm>                    <dbl> <dbl>
#>  1 Acentrella sp.     1 2008-07-17 00:00:00         0.5   0  
#>  2 Acentrella sp.     1 2008-07-17 00:00:00         1     0  
#>  3 Acentrella sp.     1 2008-07-17 00:00:00         2    62.5
#>  4 Acentrella sp.     1 2008-07-17 00:00:00         3    72.9
#>  5 Acentrella sp.     1 2008-07-17 00:00:00         4    72.9
#>  6 Acentrella sp.     1 2008-07-17 00:00:00         5    52.1
#>  7 Acentrella sp.     1 2008-07-17 00:00:00         6     0  
#>  8 Acentrella sp.     1 2008-07-17 00:00:00         7     0  
#>  9 Acentrella sp.     1 2008-07-17 00:00:00         8     0  
#> 10 Acentrella sp.     1 2008-07-17 00:00:00         9     0

The data contain all the replicates of density and body length distributions in long format.

The first step is to convert length to mass for estimating biomass patterns. Here we use the convert_length_to_mass() function, which adds a column of the individual masses based on the length-to-mass formula and coefficients in taxaInfo.

We can take a look at what this looks like:


acentrellaInfo = subset(taxaInfo, taxonID == "Acentrella sp.")
acentrellaInfo
#> # A tibble: 1 × 13
#> # Groups:   taxonID [1]
#>   taxonID  massForm      a     b percAsh method g.a   growthForm min.cpi max.cpi
#>   <chr>    <chr>     <dbl> <dbl>   <dbl> <chr>  <lgl> <chr>        <dbl>   <dbl>
#> 1 Acentre… afdm_mg… 0.0076  2.69     4.9 sf     NA    "log(g_d)…     335     365
#> # ℹ 3 more variables: pb <dbl>, min.growth <dbl>, notes <chr>

Information on its use can be viewed with ?convert_length_to_mass


acentrellaMass = convert_length_to_mass(taxaSampleList = acentrella,
                                         taxaInfo = wbtData[['taxaInfo']])

head(acentrellaMass)
#> # A tibble: 6 × 6
#>   taxonID        repID dateID              lengthClass  n_m2 afdm_mg
#>   <chr>          <dbl> <dttm>                    <dbl> <dbl>   <dbl>
#> 1 Acentrella sp.     1 2008-07-17 00:00:00         0.5   0   0.00577
#> 2 Acentrella sp.     1 2008-07-17 00:00:00         1     0   0.0372 
#> 3 Acentrella sp.     1 2008-07-17 00:00:00         2    62.5 0.240  
#> 4 Acentrella sp.     1 2008-07-17 00:00:00         3    72.9 0.716  
#> 5 Acentrella sp.     1 2008-07-17 00:00:00         4    72.9 1.55   
#> 6 Acentrella sp.     1 2008-07-17 00:00:00         5    52.1 2.83

From here you can view the size frequency histograms using the plot_cohorts() function. Check out the function options with ?plot_cohorts.


plot_cohorts(taxaSampleListMass = acentrellaMass,
             param = 'length',
             massClass = 'afdm_mg')

These figures can be helpful for identifying cohort structures and getting ballpark cohort production intervals (CPI) for species when estimating production using the size frequency method. In the future, I hope to implement an approach to delineate cohorts and estimate growth parameters from size-frequency data, but this is a feature still under development.

The next step is to estimate production. Again, only size-frequency is fully operational but I will work to update the other methods soon!

The function calc_production() is the workhorse function that will estimate community production.

To apply it to our single species example, we input the sample information and taxa information along with how many bootstraps we would like, as:


calc_production(
  taxaSampleListMass = acentrellaMass,
  infoCols = c(1:3),
  taxaInfo = acentrellaInfo,
  bootNum = 10,
  wrap = TRUE,
  taxaSummary = 'full'
)
#> $P.boots
#>            [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]    
#> P.ann.samp 7894.01  7079.587 6613.092 8324.897 7190.444 6737.9   8233.515
#> B.ann.samp 260.319  249.3877 264.5394 281.8144 246.7723 220.5484 290.1389
#> N.ann.samp 665.1926 663.8497 581.4014 701.705  596.9849 596.5369 792.4515
#>            [,8]     [,9]     [,10]   
#> P.ann.samp 8605.992 7213.12  6350.352
#> B.ann.samp 262.7356 249.5579 229.219 
#> N.ann.samp 754.534  649.3973 508.9685
#> 
#> $taxaSummary
#> $taxaSummary$summaryType
#> [1] "full"
#> 
#> $taxaSummary$taxonID
#> [1] "Acentrella sp."
#> 
#> $taxaSummary$method
#> [1] "sf"
#> 
#> $taxaSummary$P.ann.samp
#> [1] 7370.212
#> 
#> $taxaSummary$P.uncorr.samp
#> [1] 7067.326
#> 
#> $taxaSummary$cpi
#> [1] 350
#> 
#> $taxaSummary$pb
#> [1] 29.03778
#> 
#> $taxaSummary$meanN
#> [1] 657.723
#> 
#> $taxaSummary$meanB
#> [1] 253.8146
#> 
#> $taxaSummary$meanIndMass
#> [1] 0.3858989
#> 
#> $taxaSummary$Nmean
#>                   0.5         1         2         3         4         5
#> 2008-07-17   0.000000  11.45833  48.95833  33.33333 39.583333 30.208333
#> 2008-08-13   0.000000  12.35703  62.91970 105.36194 56.619625  7.291667
#> 2008-09-23  53.028394 399.27515 545.08029 132.82990  7.291667 33.333333
#> 2008-10-26  78.012266 582.35127 666.79595 173.86885  3.882576  0.000000
#> 2008-11-23   1.041667 290.29353 555.27152 123.57480 15.136670  2.083333
#> 2008-12-17 208.701307 593.97944 287.82091 103.54200 14.204545  2.083333
#> 2009-02-03  20.014881 250.22631 142.14720  27.80258  4.166667  2.083333
#> 2009-04-07  90.277778 239.72446 332.67915 102.71991 27.272727  4.924242
#> 2009-05-16   0.000000  33.44907 125.79807  32.69121 18.750000  5.208333
#> 2009-06-23   0.000000  19.50921  73.99330  80.57125 23.639759  1.041667
#>                    6         7        8     9 10 11 12 13 14 15 16 17 18 19 20
#> 2008-07-17 12.500000  2.083333 0.000000 0.000  0  0  0  0  0  0  0  0  0  0  0
#> 2008-08-13  2.083333  0.000000 0.000000 0.000  0  0  0  0  0  0  0  0  0  0  0
#> 2008-09-23 30.208333 20.833333 3.125000 3.125  0  0  0  0  0  0  0  0  0  0  0
#> 2008-10-26  2.840909  0.000000 0.000000 0.000  0  0  0  0  0  0  0  0  0  0  0
#> 2008-11-23  0.000000  0.000000 0.000000 0.000  0  0  0  0  0  0  0  0  0  0  0
#> 2008-12-17  2.083333  0.000000 0.000000 0.000  0  0  0  0  0  0  0  0  0  0  0
#> 2009-02-03  1.041667  0.000000 0.000000 0.000  0  0  0  0  0  0  0  0  0  0  0
#> 2009-04-07  6.875000  0.000000 0.000000 0.000  0  0  0  0  0  0  0  0  0  0  0
#> 2009-05-16  5.208333  0.000000 1.041667 0.000  0  0  0  0  0  0  0  0  0  0  0
#> 2009-06-23  8.798759  0.000000 0.000000 0.000  0  0  0  0  0  0  0  0  0  0  0
#>            21 22 23 24 25 26 27 28 29 30
#> 2008-07-17  0  0  0  0  0  0  0  0  0  0
#> 2008-08-13  0  0  0  0  0  0  0  0  0  0
#> 2008-09-23  0  0  0  0  0  0  0  0  0  0
#> 2008-10-26  0  0  0  0  0  0  0  0  0  0
#> 2008-11-23  0  0  0  0  0  0  0  0  0  0
#> 2008-12-17  0  0  0  0  0  0  0  0  0  0
#> 2009-02-03  0  0  0  0  0  0  0  0  0  0
#> 2009-04-07  0  0  0  0  0  0  0  0  0  0
#> 2009-05-16  0  0  0  0  0  0  0  0  0  0
#> 2009-06-23  0  0  0  0  0  0  0  0  0  0
#> 
#> $taxaSummary$Nsd
#>                   0.5         1         2         3         4         5
#> 2008-07-17   0.000000  29.64635  45.01854  27.25197 46.791335 31.614388
#> 2008-08-13   0.000000  24.55841 130.78970 151.69187 58.651399  9.882118
#> 2008-09-23  99.904573 640.40507 820.70373 166.25512  7.030714 34.652708
#> 2008-10-26 210.757419 434.06890 508.27227 144.97113  9.218577  0.000000
#> 2008-11-23   3.294039 195.29622 348.76995 107.84317 43.107151  4.392052
#> 2008-12-17 388.902788 748.41641 349.28272 114.55605 16.517014  6.588078
#> 2009-02-03  44.522062 535.29716 316.10626  53.51328  7.283395  6.588078
#> 2009-04-07 222.149872 399.20601 492.05437 170.97389 58.965796  9.319089
#> 2009-05-16   0.000000 102.16768 255.46931  57.94451 32.125093 11.251286
#> 2009-06-23   0.000000  31.63802  93.27213 107.03519 32.020550  3.294039
#>                    6         7        8        9 10 11 12 13 14 15 16 17 18 19
#> 2008-07-17 15.372183  4.392052 0.000000 0.000000  0  0  0  0  0  0  0  0  0  0
#> 2008-08-13  4.392052  0.000000 0.000000 0.000000  0  0  0  0  0  0  0  0  0  0
#> 2008-09-23 56.299275 32.200064 5.031728 7.030714  0  0  0  0  0  0  0  0  0  0
#> 2008-10-26  8.983743  0.000000 0.000000 0.000000  0  0  0  0  0  0  0  0  0  0
#> 2008-11-23  0.000000  0.000000 0.000000 0.000000  0  0  0  0  0  0  0  0  0  0
#> 2008-12-17  6.588078  0.000000 0.000000 0.000000  0  0  0  0  0  0  0  0  0  0
#> 2009-02-03  3.294039  0.000000 0.000000 0.000000  0  0  0  0  0  0  0  0  0  0
#> 2009-04-07 21.740659  0.000000 0.000000 0.000000  0  0  0  0  0  0  0  0  0  0
#> 2009-05-16 11.251286  0.000000 3.294039 0.000000  0  0  0  0  0  0  0  0  0  0
#> 2009-06-23 24.384816  0.000000 0.000000 0.000000  0  0  0  0  0  0  0  0  0  0
#>            20 21 22 23 24 25 26 27 28 29 30
#> 2008-07-17  0  0  0  0  0  0  0  0  0  0  0
#> 2008-08-13  0  0  0  0  0  0  0  0  0  0  0
#> 2008-09-23  0  0  0  0  0  0  0  0  0  0  0
#> 2008-10-26  0  0  0  0  0  0  0  0  0  0  0
#> 2008-11-23  0  0  0  0  0  0  0  0  0  0  0
#> 2008-12-17  0  0  0  0  0  0  0  0  0  0  0
#> 2009-02-03  0  0  0  0  0  0  0  0  0  0  0
#> 2009-04-07  0  0  0  0  0  0  0  0  0  0  0
#> 2009-05-16  0  0  0  0  0  0  0  0  0  0  0
#> 2009-06-23  0  0  0  0  0  0  0  0  0  0  0
#> 
#> $taxaSummary$Bmean
#>                    0.5          1         2         3         4         5
#> 2008-07-17 0.000000000  0.4267083  11.77357  23.86831 61.470320 85.519414
#> 2008-08-13 0.000000000  0.4601757  15.13101  75.44435 87.926563 20.642617
#> 2008-09-23 0.305806299 14.8690067 131.08164  95.11277 11.323480 94.366250
#> 2008-10-26 0.449884307 21.6867611 160.35199 124.49868  6.029385  0.000000
#> 2008-11-23 0.006007126 10.8105310 133.53244  88.48566 23.506256  5.897891
#> 2008-12-17 1.203547188 22.1197943  69.21556  74.14118 22.058727  5.897891
#> 2009-02-03 0.115422630  9.3184280  34.18375  19.90802  6.470560  5.897891
#> 2009-04-07 0.520617562  8.9273390  80.00313  73.55253 42.352756 13.940469
#> 2009-05-16 0.000000000  1.2456435  30.25209  23.40852 29.117520 14.744727
#> 2009-06-23 0.000000000  0.7265230  17.79401  57.69299 36.710995  2.948945
#>                     6        7        8        9 10 11 12 13 14 15 16 17 18 19
#> 2008-07-17  57.799582  14.5857  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2008-08-13   9.633264   0.0000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2008-09-23 139.682324 145.8570 31.33826 43.02554  0  0  0  0  0  0  0  0  0  0
#> 2008-10-26  13.136269   0.0000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2008-11-23   0.000000   0.0000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2008-12-17   9.633264   0.0000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2009-02-03   4.816632   0.0000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2009-04-07  31.789770   0.0000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2009-05-16  24.083159   0.0000 10.44609  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2009-06-23  40.685167   0.0000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#>            20 21 22 23 24 25 26 27 28 29 30
#> 2008-07-17  0  0  0  0  0  0  0  0  0  0  0
#> 2008-08-13  0  0  0  0  0  0  0  0  0  0  0
#> 2008-09-23  0  0  0  0  0  0  0  0  0  0  0
#> 2008-10-26  0  0  0  0  0  0  0  0  0  0  0
#> 2008-11-23  0  0  0  0  0  0  0  0  0  0  0
#> 2008-12-17  0  0  0  0  0  0  0  0  0  0  0
#> 2009-02-03  0  0  0  0  0  0  0  0  0  0  0
#> 2009-04-07  0  0  0  0  0  0  0  0  0  0  0
#> 2009-05-16  0  0  0  0  0  0  0  0  0  0  0
#> 2009-06-23  0  0  0  0  0  0  0  0  0  0  0
#> 
#> $taxaSummary$Bsd
#>                  0.5          1         2         3        4         5
#> 2008-07-17 0.0000000  1.1040302  10.82612  19.51375 72.66387 89.499937
#> 2008-08-13 0.0000000  0.9145551  31.45248 108.61887 91.08177 27.976152
#> 2008-09-23 0.5761338 23.8486848 197.36394 119.04687 10.91824 98.101384
#> 2008-10-26 1.2154045 16.1647259 122.23000 103.80648 14.31585  0.000000
#> 2008-11-23 0.0189962  7.2728313  83.87267  77.22104 66.94258 12.433845
#> 2008-12-17 2.2427404 27.8710270  83.99598  82.02779 25.64984 18.650768
#> 2009-02-03 0.2567516 19.9344664  76.01766  38.31815 11.31063 18.650768
#> 2009-04-07 1.2811029 14.8664317 118.32990 122.42575 91.57001 26.382223
#> 2009-05-16 0.0000000  3.8047243  61.43560  41.49113 49.88816 31.852250
#> 2009-06-23 0.0000000  1.1782000  22.43021  76.64248 49.72581  9.325384
#>                    6         7        8        9 10 11 12 13 14 15 16 17 18 19
#> 2008-07-17  71.08046  30.74935  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2008-08-13  20.30870   0.00000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2008-09-23 260.32597 225.43702 50.45939 96.80008  0  0  0  0  0  0  0  0  0  0
#> 2008-10-26  41.54053   0.00000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2008-11-23   0.00000   0.00000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2008-12-17  30.46305   0.00000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2009-02-03  15.23153   0.00000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2009-04-07 100.52808   0.00000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2009-05-16  52.02557   0.00000 33.03342  0.00000  0  0  0  0  0  0  0  0  0  0
#> 2009-06-23 112.75457   0.00000  0.00000  0.00000  0  0  0  0  0  0  0  0  0  0
#>            20 21 22 23 24 25 26 27 28 29 30
#> 2008-07-17  0  0  0  0  0  0  0  0  0  0  0
#> 2008-08-13  0  0  0  0  0  0  0  0  0  0  0
#> 2008-09-23  0  0  0  0  0  0  0  0  0  0  0
#> 2008-10-26  0  0  0  0  0  0  0  0  0  0  0
#> 2008-11-23  0  0  0  0  0  0  0  0  0  0  0
#> 2008-12-17  0  0  0  0  0  0  0  0  0  0  0
#> 2009-02-03  0  0  0  0  0  0  0  0  0  0  0
#> 2009-04-07  0  0  0  0  0  0  0  0  0  0  0
#> 2009-05-16  0  0  0  0  0  0  0  0  0  0  0
#> 2009-06-23  0  0  0  0  0  0  0  0  0  0  0
#> 
#> $taxaSummary$datesInfo
#>        dateID  N n_m2_mean afdm_mg_m2_mean
#> 1  2008-07-17 10  178.1250        255.4436
#> 2  2008-08-13 10  246.6333        209.2380
#> 3  2008-09-23 10 1228.1304        706.9621
#> 4  2008-10-26 10 1507.7518        326.1530
#> 5  2008-11-23 10  987.4015        262.2388
#> 6  2008-12-17 10 1212.4149        204.2700
#> 7  2009-02-03 10  447.4826         80.7107
#> 8  2009-04-07 10  804.4733        251.0866
#> 9  2009-05-16 10  222.1467        133.2977
#> 10 2009-06-23 10  207.5540        156.5586
#> 11 2009-07-16 NA  192.8395        206.0011

This will output a lot of information though it has a simple structure:

  1. P.boots–vectors of bootstrapped estimates of annual production, annual biomass, and annual abundance.

  2. taxaSummary–this comes in ‘full’ and ‘short’ versions. As you see above the full can be a lot. The ‘short’ is a paired down version of this


calc_production(
  taxaSampleListMass = acentrellaMass,
  infoCols = c(1:3),
  taxaInfo = acentrellaInfo,
  bootNum = 10,
  wrap = TRUE,
  taxaSummary = 'short'
)
#> $P.boots
#>            [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]    
#> P.ann.samp 6577.237 7689.972 8032.055 7666.049 7863.849 7524.68  7949.663
#> B.ann.samp 240.9445 270.1543 284.5143 232.0916 291.3699 249.4617 274.3334
#> N.ann.samp 588.0506 694.8745 691.2424 639.103  809.3192 647.6252 729.6995
#>            [,8]     [,9]     [,10]   
#> P.ann.samp 6265.627 7721.312 7122.196
#> B.ann.samp 214.3935 269.127  252.255 
#> N.ann.samp 451.6489 691.9201 685.241 
#> 
#> $taxaSummary
#> $taxaSummary$summaryType
#> [1] "short"
#> 
#> $taxaSummary$taxonID
#> [1] "Acentrella sp."
#> 
#> $taxaSummary$method
#> [1] "sf"
#> 
#> $taxaSummary$P.ann.samp
#> [1] 7370.212
#> 
#> $taxaSummary$cpi
#> [1] 350
#> 
#> $taxaSummary$pb
#> [1] 29.03778
#> 
#> $taxaSummary$meanN
#> [1] 657.723
#> 
#> $taxaSummary$meanB
#> [1] 253.8146
#> 
#> $taxaSummary$meanIndMass
#> [1] 0.3858989
#> 
#> $taxaSummary$datesInfo
#>        dateID  N n_m2_mean afdm_mg_m2_mean
#> 1  2008-07-17 10  178.1250        255.4436
#> 2  2008-08-13 10  246.6333        209.2380
#> 3  2008-09-23 10 1228.1304        706.9621
#> 4  2008-10-26 10 1507.7518        326.1530
#> 5  2008-11-23 10  987.4015        262.2388
#> 6  2008-12-17 10 1212.4149        204.2700
#> 7  2009-02-03 10  447.4826         80.7107
#> 8  2009-04-07 10  804.4733        251.0866
#> 9  2009-05-16 10  222.1467        133.2977
#> 10 2009-06-23 10  207.5540        156.5586
#> 11 2009-07-16 NA  192.8395        206.0011

You can also turn this off with ‘none’.

An example using the PB method

The pb method allows for multiple options to set the production to biomass ratio:

  1. A single value. This will use the same value for each bootstrap sample. The variability is derived entirely from resampling the sample units (e.g. surbers).

  2. a character string for a distribution from which to sample excluding the n=, such as rnorm(mean = 5, sd = 0.5). The function will automatically sample bootNum values.

  3. a numeric vector the same length as bootNum. This is a more explicit option to set all the values of pb that will be used. In the future, if the length \neqbootNum an error will be returned for the taxa.


acentrellaInfo$method = "pb"
acentrellaInfo$pb = "runif(min = 3, max = 8)"

calc_production(
  taxaSampleListMass = acentrellaMass,
  infoCols = c(1:3),
  taxaInfo = acentrellaInfo,
  bootNum = 10,
  wrap = TRUE,
  taxaSummary = 'short'
)
#> $P.boots
#>            [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]    
#> P.ann.samp 1088.466 1712.322 1742.765 1987.615 1809.562 2060.79  1218.325
#> B.ann.samp 270.4973 226.5819 251.4287 259.6744 259.8879 290.8709 283.018 
#> N.ann.samp 686.2019 607.9946 633.1745 629.4606 639.5113 736.3917 806.5242
#>            [,8]     [,9]     [,10]   
#> P.ann.samp 1831.69  1489.918 1000.54 
#> B.ann.samp 269.8947 256.5597 261.3129
#> N.ann.samp 720.5662 697.1245 774.2984
#> 
#> $taxaSummary
#> $taxaSummary$summaryType
#> [1] "short"
#> 
#> $taxaSummary$taxonID
#> [1] "Acentrella sp."
#> 
#> $taxaSummary$method
#> [1] "pb"
#> 
#> $taxaSummary$P.ann.samp
#> [1] 1968.429
#> 
#> $taxaSummary$pb
#> [1] 7.755384
#> 
#> $taxaSummary$meanN
#> [1] 657.723
#> 
#> $taxaSummary$meanB
#> [1] 253.8146
#> 
#> $taxaSummary$meanIndMass
#> [1] 0.3858989
#> 
#> $taxaSummary$datesInfo
#>        dateID  N n_m2_mean afdm_mg_m2_mean
#> 1  2008-07-17 10  178.1250        255.4436
#> 2  2008-08-13 10  246.6333        209.2380
#> 3  2008-09-23 10 1228.1304        706.9621
#> 4  2008-10-26 10 1507.7518        326.1530
#> 5  2008-11-23 10  987.4015        262.2388
#> 6  2008-12-17 10 1212.4149        204.2700
#> 7  2009-02-03 10  447.4826         80.7107
#> 8  2009-04-07 10  804.4733        251.0866
#> 9  2009-05-16 10  222.1467        133.2977
#> 10 2009-06-23 10  207.5540        156.5586
#> 11 2009-07-16 NA  192.8395        206.0011

Applying this to multiple taxa

The current internal code to do this with just calc_production is still buggy. But you can run production by splitting the full sampleInfo into a list by taxa and using apply or functions from the (purrr package)[https://purrr.tidyverse.org/]

The code below shows how to accomplish this. Again, only size frequency and pb are accepted currently.

# load needed packages
library(magrittr)
library(purrr)

# The sampleInfo is already a list by taxonID
str(wbtData[['sampleInfo']], list.len = 5)
#> list<tibble[,5]> [1:32] 
#> $ Acentrella sp.    : tibble [3,100 × 5] (S3: tbl_df/tbl/data.frame)
#>  ..$ taxonID    : chr [1:3100] "Acentrella sp." "Acentrella sp." "Acentrella sp." "Acentrella sp." ...
#>  ..$ repID      : num [1:3100] 1 1 1 1 1 1 1 1 1 1 ...
#>  ..$ dateID     : POSIXct[1:3100], format: "2008-07-17" "2008-07-17" ...
#>  ..$ lengthClass: num [1:3100] 0.5 1 2 3 4 5 6 7 8 9 ...
#>  ..$ n_m2       : num [1:3100] 0 0 62.5 72.9 72.9 ...
#> $ Ameletus spp.     : tibble [3,038 × 5] (S3: tbl_df/tbl/data.frame)
#>  ..$ taxonID    : chr [1:3038] "Ameletus spp." "Ameletus spp." "Ameletus spp." "Ameletus spp." ...
#>  ..$ repID      : num [1:3038] 1 1 1 1 1 1 1 1 1 1 ...
#>  ..$ dateID     : POSIXct[1:3038], format: "2008-07-17" "2008-07-17" ...
#>  ..$ lengthClass: num [1:3038] 0.5 1 2 3 4 5 6 7 8 9 ...
#>  ..$ n_m2       : num [1:3038] 0 0 0 0 0 0 0 0 0 0 ...
#> $ Attenella spp.    : tibble [3,038 × 5] (S3: tbl_df/tbl/data.frame)
#>  ..$ taxonID    : chr [1:3038] "Attenella spp." "Attenella spp." "Attenella spp." "Attenella spp." ...
#>  ..$ repID      : num [1:3038] 1 1 1 1 1 1 1 1 1 1 ...
#>  ..$ dateID     : POSIXct[1:3038], format: "2008-07-17" "2008-07-17" ...
#>  ..$ lengthClass: num [1:3038] 0.5 1 2 3 4 5 6 7 8 9 ...
#>  ..$ n_m2       : num [1:3038] 0 0 41.7 20.8 0 ...
#> $ Baetis spp.       : tibble [3,100 × 5] (S3: tbl_df/tbl/data.frame)
#>  ..$ taxonID    : chr [1:3100] "Baetis spp." "Baetis spp." "Baetis spp." "Baetis spp." ...
#>  ..$ repID      : num [1:3100] 1 1 1 1 1 1 1 1 1 1 ...
#>  ..$ dateID     : POSIXct[1:3100], format: "2008-07-17" "2008-07-17" ...
#>  ..$ lengthClass: num [1:3100] 0.5 1 2 3 4 5 6 7 8 9 ...
#>  ..$ n_m2       : num [1:3100] 0 0 20.8 31.2 93.8 ...
#> $ Barbaetis spp.    : tibble [3,100 × 5] (S3: tbl_df/tbl/data.frame)
#>  ..$ taxonID    : chr [1:3100] "Barbaetis spp." "Barbaetis spp." "Barbaetis spp." "Barbaetis spp." ...
#>  ..$ repID      : num [1:3100] 1 1 1 1 1 1 1 1 1 1 ...
#>  ..$ dateID     : POSIXct[1:3100], format: "2008-07-17" "2008-07-17" ...
#>  ..$ lengthClass: num [1:3100] 0.5 1 2 3 4 5 6 7 8 9 ...
#>  ..$ n_m2       : num [1:3100] 0 0 20.8 0 0 ...
#>  [list output truncated]
#> @ ptype: tibble [0 × 5] (S3: tbl_df/tbl/data.frame)
#>  ..$ taxonID    : chr(0) 
#>  ..$ repID      : num(0) 
#>  ..$ dateID     : 'POSIXct' num(0) 
#>  - attr(*, "tzone")= chr "UTC"
#>  ..$ lengthClass: num(0) 
#>  ..$ n_m2       : num(0)

# if you have a large data frame with all taxa the code below shows how to do this:
## install a package to split into a named_list
# devtools::install_github('jimjunker1/junkR')
# library(junkR)
## Convert to data frame
# sampleInfoDf = wbtData[['sampleInfo']] %>% bind_rows
## Convert back to named list split by taxonID
# sampleInfoList = sampleInfoDf %>% junkR::named_group_split(taxonID)

sampleInfoList = wbtData[['sampleInfo']]

sampleInfoListMass = sampleInfoList %>% 
  purrr::map(
    ~convert_length_to_mass(.x,
                            taxaInfo = taxaInfo)
  )

fullProduction = sampleInfoListMass %>% 
  purrr::map2(., list(taxaInfo),
                     ~calc_production(.x,
                                      infoCols = c(1:3),
                                      taxaInfo = .y, 
                                      bootNum = 10, 
                                      wrap = TRUE,
                                      taxaSummary = 'short'
                                      )
                     )

str(fullProduction, list.len = 5)
#> List of 32
#>  $ Acentrella sp.    :List of 2
#>   ..$ P.boots    :List of 30
#>   .. ..$ : num 6968
#>   .. ..$ : num 244
#>   .. ..$ : num 612
#>   .. ..$ : num 7217
#>   .. ..$ : num 257
#>   .. .. [list output truncated]
#>   .. ..- attr(*, "dim")= int [1:2] 3 10
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : chr [1:3] "P.ann.samp" "B.ann.samp" "N.ann.samp"
#>   .. .. ..$ : NULL
#>   ..$ taxaSummary:List of 10
#>   .. ..$ summaryType: chr "short"
#>   .. ..$ taxonID    : chr "Acentrella sp."
#>   .. ..$ method     : chr "sf"
#>   .. ..$ P.ann.samp : num 7370
#>   .. ..$ cpi        : num 350
#>   .. .. [list output truncated]
#>  $ Ameletus spp.     :List of 2
#>   ..$ P.boots    :List of 30
#>   .. ..$ : num 796
#>   .. ..$ : num 41.6
#>   .. ..$ : num 31.7
#>   .. ..$ : num 736
#>   .. ..$ : num 41.9
#>   .. .. [list output truncated]
#>   .. ..- attr(*, "dim")= int [1:2] 3 10
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : chr [1:3] "P.ann.samp" "B.ann.samp" "N.ann.samp"
#>   .. .. ..$ : NULL
#>   ..$ taxaSummary:List of 10
#>   .. ..$ summaryType: chr "short"
#>   .. ..$ taxonID    : chr "Ameletus spp."
#>   .. ..$ method     : chr "sf"
#>   .. ..$ P.ann.samp : num 660
#>   .. ..$ cpi        : num 350
#>   .. .. [list output truncated]
#>  $ Attenella spp.    :List of 2
#>   ..$ P.boots    :List of 30
#>   .. ..$ : num 2617
#>   .. ..$ : num 112
#>   .. ..$ : num 71.7
#>   .. ..$ : num 2814
#>   .. ..$ : num 125
#>   .. .. [list output truncated]
#>   .. ..- attr(*, "dim")= int [1:2] 3 10
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : chr [1:3] "P.ann.samp" "B.ann.samp" "N.ann.samp"
#>   .. .. ..$ : NULL
#>   ..$ taxaSummary:List of 10
#>   .. ..$ summaryType: chr "short"
#>   .. ..$ taxonID    : chr "Attenella spp."
#>   .. ..$ method     : chr "sf"
#>   .. ..$ P.ann.samp : num 2534
#>   .. ..$ cpi        : num 350
#>   .. .. [list output truncated]
#>  $ Baetis spp.       :List of 2
#>   ..$ P.boots    :List of 30
#>   .. ..$ : num 3568
#>   .. ..$ : num 139
#>   .. ..$ : num 221
#>   .. ..$ : num 3821
#>   .. ..$ : num 133
#>   .. .. [list output truncated]
#>   .. ..- attr(*, "dim")= int [1:2] 3 10
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : chr [1:3] "P.ann.samp" "B.ann.samp" "N.ann.samp"
#>   .. .. ..$ : NULL
#>   ..$ taxaSummary:List of 10
#>   .. ..$ summaryType: chr "short"
#>   .. ..$ taxonID    : chr "Baetis spp."
#>   .. ..$ method     : chr "sf"
#>   .. ..$ P.ann.samp : num 3833
#>   .. ..$ cpi        : num 350
#>   .. .. [list output truncated]
#>  $ Barbaetis spp.    :List of 2
#>   ..$ P.boots    :List of 30
#>   .. ..$ : num 1787
#>   .. ..$ : num 66.4
#>   .. ..$ : num 92.4
#>   .. ..$ : num 1495
#>   .. ..$ : num 52.7
#>   .. .. [list output truncated]
#>   .. ..- attr(*, "dim")= int [1:2] 3 10
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : chr [1:3] "P.ann.samp" "B.ann.samp" "N.ann.samp"
#>   .. .. ..$ : NULL
#>   ..$ taxaSummary:List of 10
#>   .. ..$ summaryType: chr "short"
#>   .. ..$ taxonID    : chr "Barbaetis spp."
#>   .. ..$ method     : chr "sf"
#>   .. ..$ P.ann.samp : num 1512
#>   .. ..$ cpi        : num 350
#>   .. .. [list output truncated]
#>   [list output truncated]
Junker, J. R., and W. F. Cross. 2014. Seasonality in the trophic basis of a temperate stream invertebrate assemblage: Importance of temperature and food quality. Limnology and Oceanography 59:507–518.