Incorrect estimates of integrated density from `integrate_over_size_distribution`

[Becca-90]

library(plant)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(ggplot2)
library(purrr)

Let’s assume a known density wrt height, to test against. Let

$$n(H) = a \, H ^{−b}$$

Then we can directly integrate to compare results against ‘plant’ integration

$$N= \int_{1}^{20} n(H) \, dH = a \bigg[\frac{H^{1-b}}{1-b}\bigg]^{20}_{1}$$ $$H_{T}= \int_{1}^{20}N(H) \, H \, dH = a \bigg[\frac{H^{2-b}}{2-b}\bigg]^{20}_{1}$$

Now generate data to tests

## set density = a*H^-b, height sequence 1:20, b = 3/2, a = starting density per unit area
## create stand of individuals with heights obtained and resulting density over time
## so we have constant height growth increase, therefore assume time steps are non-linear

a = 10
b = 3/2
plant_stand <- 
  tibble(
    time = 1,
    step=2,
    patch_density = 1, 
    species = 1, 
    height = seq(20,1, by = -0.05)
    ) %>% 
  mutate(
    cohort = seq_len(n()) %>% rev(),
    density = a*height^(-b)
  )

Now compare solutions

# integrate with plant inbuilt solver
stand_integrate <- integrate_over_size_distribution(plant_stand)

# Analyitcal soluntions
# N: a/(1-b) H ^ (1-b)
# H: a/(2-b) H ^ (2-b)

N_analytical <- a/(1-b) *20 ^ (1-b) - a/(1-b) *1 ^ (1-b)
N_trapezium <- -plant:::trapezium(plant_stand$height, plant_stand$density)
N_plant<- stand_integrate$individuals

c(N_analytical, N_trapezium, N_plant)

## [1]  15.52786  15.53099 241.21153

H_analytical <- a/(2-b) *20 ^ (2-b) - a/(2-b) *1 ^ (2-b)
H_trapezium <- -plant:::trapezium(plant_stand$height, plant_stand$height*plant_stand$density)
H_plant <- stand_integrate$height

c(H_analytical, H_trapezium, H_plant)

## [1] 69.44272 69.44375 69.44375

Hav_analytical <- H_analytical/N_analytical
Hav_trapezium <- H_trapezium/N_trapezium
Hav_plant <- H_plant/N_plant

c(Hav_analytical, Hav_trapezium, Hav_plant)

## [1] 4.4721360 4.4713033 0.2878956

So it turns out integrate_over_size_distreibution is getting the number of individuals wrong which is flowing onto average height

Here’s a fix:

integrate_over_size_distribution2 <- function (tidy_species_data) 
{
    tidy_species_data %>% dplyr::select(-.data$cohort) %>% stats::na.omit() %>% 
        dplyr::filter(.data$step > 1) %>% dplyr::group_by(.data$step, 
        .data$time, .data$patch_density, .data$species) %>% 
    dplyr::summarise(
      density_integrated = -plant:::trapezium(.data$height, .data$density), 
      min_height = min(.data$height),
      dplyr::across(where(is.double) & !c(density, density_integrated, min_height) , ~-plant:::trapezium(height, density * .x)), 
        .groups = "drop") %>% 
    rename(density = density_integrated)
}

integrate_over_size_distribution2(plant_stand) %>% 
  mutate(Hav = height/density)

## # A tibble: 1 × 8
##    step  time patch_density species density min_height height   Hav
##   <dbl> <dbl>         <dbl>   <dbl>   <dbl>      <dbl>  <dbl> <dbl>
## 1     2     1             1       1    15.5          1   69.4  4.47

[dfalster]

@pzylstra @fjrrobinson @aornugent @itowers1

This issue may have effected all of you. Have you used the function integrate_over_size_distribution to calculate either density of individuals per unit area, or average size (by any metric). If so, numbers seem to be incorrect, due to incorrect calculation of total number of individuals. But it seems total amounts of area or biomass (e.g. leaf area, biomass) are OK.

We have a proposed fix. Just waiting to

1. Add test for the function
2. See if this improves @Becca-90 's self-thinning plots

Thanks @Becca-90 for investigating

[Becca-90]

So integrate_over_size_density will be changed in the develop branch, and then I'll redo my self-thinning stands? Or should I just test directly first with the function on my end? I've got one example stand which seems to look as you'd expect (yay!) but something is going on when I do the same with my trait-varied stands... but looks like the integrate function is fine at least!!

image

[Becca-90]

Actually, that's on a log scale so the diameter is still getting to a huge size I think...will check it over

[dfalster]

addressed via #345 Test added via 974b17be395f662a0f78f66f0e86802463438a69