Aim

This vignette exposes the logic of the calculations to scale up from individual excretion to total population recycling, as well as the calculation of its uncertainty using Monte Carlo. simulation.

Individual Allometry

The first step is to estimate the allometric relationship between individual size \(x\) and nitrogen excretion rate \(E_N\), which can be simply done by fitting as a log-log linear model.

\[\log_{10}(E_N) = a + b \cdot \log_{10}(x) \Leftrightarrow E_N = a \cdot x^b \]

For example, for Upper La Laja’s second year:

load('../data/data_excretion.rda')
edata <- subset(data_excretion, stream == 'UL' & month_intro == 14)

op <- par(mfrow = c(1,2))
plot(N_rate ~ size, data = edata,
     pch = 16, xlab = 'Mass (g)', ylab ='Excretion Rate (gN/h)', main = 'natural scale')
plot(N_rate ~ size, data = edata, log = "xy",
     pch = 16, xlab = 'Mass (g)', ylab ='Excretion Rate (gN/h)', main = 'log-log scale')
#> Warning in xy.coords(x, y, xlabel, ylabel, log): 3 x values <= 0 omitted from
#> logarithmic plot

par(op)

This is implemented by function allo_fit of our package allopop. You can access the help typing ‘?allo_fit’ For example:

library(allopop)

fit <- allo_fit(response = N_rate, size = size, data = edata)

fit
#> $coef
#> (Intercept)   log(size) 
#>   3.2254379   0.6202035 
#> 
#> $vcov
#>             (Intercept)   log(size)
#> (Intercept) 0.015757805 0.006071351
#> log(size)   0.006071351 0.002952926
#> 
#> $fit
#> 
#> Call:
#> lm(formula = log(response) ~ log(size))
#> 
#> Coefficients:
#> (Intercept)    log(size)  
#>      3.2254       0.6202

Note that this function returns a list with two elements: a vector of coefficients (\(a\) and \(b\)), and a variance-covariance matrix of the estimates.


plot(N_rate ~ size, data = edata,
     pch = 16, xlab = 'Mass (g)', ylab ='Excretion Rate (gN/h)', log = "xy")
#> Warning in xy.coords(x, y, xlabel, ylabel, log): 3 x values <= 0 omitted from
#> logarithmic plot

Since we have the mean estimates of \(a\) and \(b\), as well as their variance-covariance matrix, we can simulate many pairs of coefficients using a multivariate (bivariate) normal distribution. To do so, we have created function allo_mc

mc <- allo_mc(fit, n = 1000)
head(mc$coefs)
#>      (Intercept) log(size)
#> [1,]    3.345309 0.6843555
#> [2,]    3.209217 0.6053073
#> [3,]    3.418213 0.6597732
#> [4,]    3.165267 0.6448216
#> [5,]    3.217008 0.5995975
#> [6,]    3.268599 0.6700945

, which can also create predictions with associated uncertainty if we set . Let’s show an example with only 10 simulations for the sake of illustration.

mc_ci <- allo_mc(fit, n = 10, predict = TRUE, xpred = seq(0, 1, 0.1))

mc_ci
#> $coefs
#>       (Intercept) log(size)
#>  [1,]    3.484668 0.7525020
#>  [2,]    3.389429 0.6858017
#>  [3,]    3.344982 0.6378859
#>  [4,]    3.084696 0.5220515
#>  [5,]    3.112757 0.6245347
#>  [6,]    3.047597 0.5615681
#>  [7,]    2.988688 0.5154992
#>  [8,]    3.349101 0.6821404
#>  [9,]    3.351012 0.6382618
#> [10,]    3.365344 0.6489440
#> 
#> $xpred
#>  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
#> 
#> $pred_ci
#>       [,1]     [,2]      [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
#> 2.5%     0 5.433712  8.296738 10.61703 12.43011 13.97828 15.38552 16.68547
#> 50%      0 6.086169  9.606722 12.75526 15.52487 17.98689 20.28603 22.45803
#> 97.5%    0 6.568876 10.207539 13.24609 16.27632 19.15515 21.90754 24.54806
#>           [,9]    [,10]    [,11]
#> 2.5%  17.90011 19.04480 20.13074
#> 50%   24.52684 26.50947 28.41859
#> 97.5% 27.09169 29.55356 31.94502

Note that you get a list of three objects:

This can be used, for example, to produce the fit plot.

plot(N_rate ~ size, data = edata,
     pch = 16, cex= 0.7, xlab = 'Mass (g)', ylab = 'Excretion Rate (gN/h)', log = 'xy'
     )
#> Warning in xy.coords(x, y, xlabel, ylabel, log): 3 x values <= 0 omitted from
#> logarithmic plot

polygon(x = 10^(c(mc$xpred, mc$xpred)), y = 10^(c(mc$pred_ci[1,], mc$pred_ci[3,])))

Demography

The first step is to obtain the population census data to extract the size distribution and the population density. This can be found in the package by loading the processed mark-recapture data data_population and applying function calc_demog.

load('../data/data_population.rda')
pdata <- subset(data_population, stream == 'UL' & month_intro == 14)

size_structure <- calc_sizestr(pdata)
head(size_structure)
#>       size prop
#> [1,] 0.001    0
#> [2,] 0.002    0
#> [3,] 0.003    0
#> [4,] 0.004    0
#> [5,] 0.005    0
#> [6,] 0.006    0

We can plot this to get a visual:

plot(size_structure[,1], size_structure[,2],
     type = 'h',
     xlab = 'Weight (g)',
     ylab = 'Proportion of the population',
     xlim = c(0, 1))

The population density (in individuals per benthic \(m^2\)) is already variable pop_dens in data_population

str(pdata)
#> tibble [748 × 7] (S3: tbl_df/tbl/data.frame)
#>  $ stream     : Factor w/ 4 levels "CA","LL","TY",..: 4 4 4 4 4 4 4 4 4 4 ...
#>  $ month_intro: num [1:748] 14 14 14 14 14 14 14 14 14 14 ...
#>  $ month_cal  : num [1:748] 4 4 4 4 4 4 4 4 4 4 ...
#>  $ area       : num [1:748] 482 482 482 482 482 ...
#>  $ n_ind      : Named int [1:748] 748 748 748 748 748 748 748 748 748 748 ...
#>   ..- attr(*, "names")= chr [1:748] "n" "n" "n" "n" ...
#>  $ size       : num [1:748] 0.633 0.851 0.579 0.633 0.635 0.492 0.654 0.829 0.747 0.141 ...
#>  $ pop_dens   : Named num [1:748] 1.55 1.55 1.55 1.55 1.55 ...
#>   ..- attr(*, "names")= chr [1:748] "n" "n" "n" "n" ...

Scaling up to the Population Effect

The scaling up of from individual-level excretion \(E_{z,t}\) at time \(t\) and size \(z\) to population-level recycling \(R_t\) can be defined by the following equation:

\[ R_t =\int_0^\infty E_{z,t} \cdot P_{z,t} \cdot N_t \cdot dz \]

, where \(P_{z,t}\) is the proportion of individuals of size \(z\) at time \(t\); and \(N_t\) is the population size (or areal density) at time \(t\).

This is implememted in function ind2pop, which requires the three above elemets nicely bundled into a list of three elements. To do so, we can use pop_bundle, which delivers precisely that, given a stream and a month since introduction.

mybundle <- pop_bundle(loc = 'UL', month_int = 14, nsims = 1000)
Rt <- ind2pop(params = mybundle)

mean(Rt)
#> [1] 13.55228

The projected average recycling is therefore \(R_t =\) 13.5522781 \(gN \cdot h^{-1} \cdot m^{-2}\). Because the result of ind2pop is a vector of estimates resulting from the 10000 Monte Carlo simulations, we can plot the distribution of expected values and calculate the 95% confidenc intervals, for example.

plot(density(Rt), main = '',
     xlab = expression(gN/h/m^2)
     )
abline(v = mean(Rt), lty = 2)


quantile(Rt, c(0.025, 0.975))
#>     2.5%    97.5% 
#> 11.97860 15.29825