Bugs in makeLifeTable, dEntropy, and kEntropy

[patrickbarks]

I've noticed a few potential bugs in makeLifeTable, some of which also apply to dEntropy, and kEntropy.

Applies to makeLifeTable, dEntropy, and kEntropy

1. I think the loop to calculate age-specific survivorship

  matUtemp = matU
  # ...
  for (o in 1:nSteps) {
    survivorship[o,] = colSums(matUtemp %*% matU)
    matUtemp = matUtemp %*% matU
  }

should instead be

  matUtemp = matU
  # ...
  for (o in 1:nSteps) {
    survivorship[o,] = colSums(matUtemp)
    matUtemp = matUtemp %*% matU
  }

As written, the code effectively skips over the first age class, for which survivorship should be colSums(matU).

Caswell's Matlab code for lx is (Caswell 2001, pg 120):

  for x = 1:150
    surv(x,:) = sum(T^x);

I think the current R code is effectively doing this:

  for x = 1:150
    surv(x,:) = sum(T^(x+1));

Applies to makeLifeTable and dEntropy

2. In the loops to calculate age-specific fertility (fertMatrix) and clonality (clonMatrix), there are a few potential issues (these only seem to cause problems when startLife > 1):

• diag(t(e) %*% matUtemp2) is supposed to create a square matrix with vector t(e) %*% matUtemp2 along the diagonal, but because t(e) %*% matUtemp2 returns a 1-row matrix, diag(...) instead pulls out the single 'diagonal' element of that matrix, [1,1]. This could be fixed using as.numeric(), or by using colSums(matUtemp2) instead of the matrix multiplication.

• The single instance of scalar multiplication should be changed to matrix multiplication (Caswell 2001 pg. 120).

• The matrix inverse fails anytime a column within matU contains all zeros (even if it's only the last column). That whole operation is really just dividing columns of matUtemp2 by their sum, which can be coded without the matrix inverse (Caswell 2001 pg. 123).

Applies only to makeLifeTable

3. In the calculation of Lx, I think the substract operation should instead be addition, to get the average survivorship between age x and x+1. L[x] = (l[x] + l[x+1]) / 2

4. I think the calculation for Tx should be sum(L[x]:L[maxAge]), instead of sum(L[1]:L[x]) (cumsum is doing the latter).

[patrickbarks]

Few more issues:

Applies to dEntropy

1. The final calculation in dEntropy is:

abs(sum(lxmx * log(lxmx)) / sum(lxmx))

When lxmx == 0, the lxmx * log(lxmx) part yields NaN (i.e. 0 * -Inf). To get around this, instances of lxmx == 0 are first converted to 1. This trick works fine in the numerator, because the limit of x * log(x) as x approaches 0 is 0, and correspondingly 1 * log(1) == 0.

But the extra 1s added to lxmx are also getting carried into the denominator (R0 = sum(lxmx)), which is not correct. Suggest editing so that denominator is calculated before 0s are converted to 1s in lxmx.

2. There's a lot of code duplication in the final three blocks that calculate entropy with respect to the different types of reproduction: Fec, Clo, and FecClo. Suggest pulling the duplicated sections into a separate utility function that calculates entropy given vectors lx and fx (where fx is one of mx, cx, or mxcx).

Applies to kEntropy

3. When argument trapeze == TRUE, I think the kEntropy function currently implements the 'midpoint approximation' to a definite integral rather than the 'trapezoidal approximation', as advertised.

Some details on the two approximations below:
http://tutorial.math.lamar.edu/Classes/CalcII/ApproximatingDefIntegrals.aspx
https://www.r-bloggers.com/the-trapezoidal-rule-of-numerical-integration-in-r/

4. The kEntropy function returns NaN if any values of lx are 0 (due to the same log(0) issue described for dEntropy). Suggest adding code to omit late age-classes at/after which survivorship has fallen to 0.