The Mathematics

As we defined, population-level recycling \(R_t\) can be defined by the following equation:

\[ R_t =\int E_{z,t} \cdot P_{z,t} \cdot N_t \cdot dz \] Since the population-level excretion \(R_t\) at a given time \(t\) is the integral of individual excretion \(E_{z,t}\), size-structure \(P_{z,t}\), and density \(N_t\) at that time, the proportional change in excretion across sampling step \(t\) to \(t+1\) is

\[ \frac{R_{t+1}}{R_t} = \frac{\int E_{z,t+1} \cdot P_{z,t+1} \cdot N_{t+1} \cdot dz}{\int E_{z,t} \cdot P_{z,t} \cdot N_t \cdot dz} \]

This can be decomposed into the product of four contributions

\[ \frac{R_{t+1}}{R_t} = A \cdot S \cdot D \cdot I_{A,S} \] where \(A\) is proportional contribution of a change in individual excretion rate, \(S\) is the contribution of population density, \(D\) is the contribution of density, and \(I_{A,S}\) is the contribution of the interaction between density and size structure.

\[ A = \frac{\int E_{z,t+1} \cdot P_{z,t} \cdot N_{t} \cdot dz}{\int E_{z,t} \cdot P_{z,t} \cdot N_{t} \cdot dz} = \frac{\int E_{z,t+1} \cdot P_{z,t} \cdot dz}{\int E_{z,t} \cdot P_{z,t+1} \cdot dz} \] \[ S = \frac{\int E_{z,t} \cdot P_{z,t+1} \cdot N_{t} \cdot dz}{\int E_{z,t} \cdot P_{z,t} \cdot N_{t} \cdot dz} = \frac{\int E_{z,t} \cdot P_{z,t+1} \cdot dz}{\int E_{z,t} \cdot P_{z,t} \cdot dz} \] \[ D = \frac{\int E_{z,t} \cdot P_{z,t} \cdot N_{t+1} \cdot dz}{\int E_{z,t} \cdot P_{z,t} \cdot N_{t} \cdot dz} = \frac{N_{t+1}}{N_t} \] Finally, because excretion depends on size, there is an interaction between size structure and size-specific excretion, which can be calculated by multiplying the total change caused by A and P togeter, by the inverse of the respective contributions of A and S.

\[ I_{A,S} = \frac{\int E_{z,t+1} \cdot P_{z,t+1} \cdot dz}{\int E_{z,t} \cdot P_{z,t} \cdot dz} \cdot A^{-1} \cdot S^{-1} = \frac{\int E_{z,t+1} \cdot P_{z,t+1} \cdot dz}{\int E_{z,t} \cdot P_{z,t} \cdot dz} \cdot \frac{\int E_{z,t} \cdot P_{z,t} \cdot dz}{\int E_{z,t+1} \cdot P_{z,t} \cdot dz} \cdot \frac{\int E_{z,t} \cdot P_{z,t} \cdot dz}{\int E_{z,t} \cdot P_{z,t+1} \cdot dz} \]

Function Implemenation

In order to make things easier, I have created a function that calculates and ‘budles up’ the three main components needed for the MC simulation of population excretion, namely:

This is what function pop_bundle conveniently extracts, given a stream and monthsince introduction:

A <- pop_bundle(loc = 'CA', month_int = 14, nsims = 1000)
str(A)
## List of 3
##  $ coefs   : num [1:1000, 1:2] 0.2559 -0.0267 0.8103 2.9604 0.5185 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "(Intercept)" "log(size)"
##  $ size_str: num [1:1999, 1:2] 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "size" "prop"
##  $ pop_dens: Named num 1.93
##   ..- attr(*, "names")= chr "n"
B <- pop_bundle(loc = 'CA', month_int = 25, nsims = 1000)
str(B)
## List of 3
##  $ coefs   : num [1:1000, 1:2] 2.76 2.74 3.1 3.12 3.33 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "(Intercept)" "log(size)"
##  $ size_str: num [1:1999, 1:2] 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "size" "prop"
##  $ pop_dens: Named num 2.03
##   ..- attr(*, "names")= chr "n"

With these two bundles of parameters, we can now use the function change_decomp to do the actual decomposition of the difference in population recycling into its contributirs:

decomp <- change_decomp(A, B)
head(decomp)
##     change         A        S        D         I
## 1 16.32289 12.524650 1.232076 1.048295 1.0090439
## 2 25.39134 18.636890 1.212414 1.048295 1.0719574
## 3 22.91801 14.432357 1.276766 1.048295 1.1864368
## 4  8.94860  3.744827 1.311532 1.048295 1.7380446
## 5 19.54881 14.420625 1.377349 1.048295 0.9388766
## 6 18.95967 11.674981 1.312907 1.048295 1.1799320

The results shows a matrix with columns for the actual change (in ratio), as well as the 4 contributions \(A\), \(S\), \(D\), and \(I_{A,S}\), and as many columns as MC simulations. This allows for the calculation of the error or 95% CI:

apply(decomp, 2, quantile, c(0.025, 0.5, 0.975))
##          change        A        S        D         I
## 2.5%   8.178792  4.54620 1.205764 1.048295 0.7323777
## 50%   20.691820 14.02271 1.311689 1.048295 1.0716280
## 97.5% 53.309012 42.21942 1.430875 1.048295 1.5204008

And as a quick quality check, we know that the change should equal the product if the four components.

product <- apply(decomp[,c('A','S','D','I')], 1, prod)

plot(decomp[,'change'] ~ product,
     ylab = 'Proportional Change',
     xlab = 'Product of Components',
     pch = 16, cex = 0.8)
abline(0, 1, lty = 2)