Cumulative mortality does not increase monotonically with time in some cases.

[itowers1]

Mortality is non-monotonic in certain instances when looking at communities of two species. In the standard one species case, we find that all mortality~time curves are monotonic (i.e. monotonic object below is all TRUE). This was conducted in the tidy_patch branch.

devtools::load_all()

base_parameters <- function() {
  p0 <- scm_base_parameters("FF16", "FF16_Env")
}

p0 = base_parameters()
hyper_par_fn = make_FF16_hyperpar(latitude = 25)

traits =  trait_matrix(c(0.24), c("lma"))

p1 <- expand_parameters(traits, p0, hyper_par_fn, mutant = FALSE)

result <- build_schedule(p1)
result <- run_scm_collect(result)

monotonic <- rep(NA, ncol(result$species[[1]]["mortality",,]))

for(i in 1:ncol(result$species[[1]]["mortality",,])){
x<-result$species[[1]]["mortality",,i]
x[!is.na(x)] -> x
monotonic[i] <- all(x == cummax(x))
}

monotonic

However, if we inspect a two species case (and just species one in this case), we find at least ten cohorts, such as cohort 100, with non-monotonic curves.

base_parameters <- function() {
  p0 <- scm_base_parameters("FF16", "FF16_Env")
}

p0 = base_parameters()
hyper_par_fn = make_FF16_hyperpar(latitude = 25)

traits =  trait_matrix(c(0.07, 0.24), c("lma"))
p1 <- expand_parameters(traits, p0, hyper_par_fn, mutant = FALSE)

result <- build_schedule(p1)
result <- run_scm_collect(result)

result$species[[1]]["mortality",,100]

monotonic <- rep(NA, ncol(result$species[[1]]["mortality",,]))

for(i in 1:ncol(result$species[[1]]["mortality",,])){
x<-result$species[[1]]["mortality",,i]
x[!is.na(x)] -> x
monotonic[i] <- all(x == cummax(x))
}

Is this possibly related to the relative trait values of the species. If we hold species 1 constant at lma = 0.07, and compare to higher and lower trait values, it looks like the proporiton of cohorts with monotonic mortality~time curves declines as the difference in trait values increases

upper_trait_value <- seq(from=0.04, to=0.24,by= 0.01)

find_mortality_monotonic <- function(upper_trait_value) {
traits =  trait_matrix(c(0.07, upper_trait_value), c("lma"))
p1 <- expand_parameters(traits, p0, hyper_par_fn, mutant = FALSE)

result <- build_schedule(p1)
result <- run_scm_collect(result)

monotonic_sp1 <- rep(NA, ncol(result$species[[1]]["mortality",,]))

for(i in 1:ncol(result$species[[1]]["mortality",,])){
result$species[[1]]["mortality",,i] -> x
x[!is.na(x)] -> x
monotonic_sp1[i] <- all(x == cummax(x))
}

monotonic_sp1 <- sum(monotonic_sp1)/length(monotonic_sp1)


monotonic_sp2 <- rep(NA, ncol(result$species[[2]]["mortality",,]))

for(i in 1:ncol(result$species[[2]]["mortality",,])){
result$species[[2]]["mortality",,i] -> x
x[!is.na(x)] -> x
monotonic_sp2[i] <- all(x == cummax(x))
}
monotonic_sp2 <- sum(monotonic_sp2)/length(monotonic_sp2)

tibble(sp1=monotonic_sp1, sp2=monotonic_sp2)
}

In the figure below, the orange line represents when sp1 and sp2 have no difference in LMA (i.e. both lmas are 0.7)

out <- purrr::map(upper_trait_value, find_mortality_monotonic)
out %>% bind_rows() %>% mutate(species1_trait = 0.07, species2_trait = upper_trait_value) -> out

out %>%
  pivot_longer(cols=c(monotonic_sp1,monotonic_sp2), names_to = "species") %>%
  mutate(sp2_sp1_lma_diff = species2_trait-species1_trait) %>%
  ggplot() +
  geom_line(aes(x=sp2_sp1_lma_diff, y = value, group = species, col=species), size=1.5) +
  xlim(-0.05,0.15) +
  theme_classic() +
  geom_vline(xintercept=0, col="orange", size=1, linetype=2) +
  ylab("Proportion of cohorts with monotonic mortality~time curves")

[dfalster]

Very nicely documented @itowers1 ! The MWE is awesome and I can repeat behaviour here. In future plans also post a pic of output. Here is a figure generated from the code above:

This behaviour is curious. My initial sense is that cumulative mortality should never decrease. As per eq 24-25 in https://traitecoevo.github.io/plant/articles/demography.html, cumulative mortality is solved as an initial value problem, with initial value of $y(0) = - \ln\left(S_{\rm G} (x, H_0, E_{a0})\right)$ and $\frac{dy}{dt} = d(x, H_i(t) , E_t)$. So the only way mortality could decrease is

1. $d<0$ shouldn't happen but who knows
2. There's an error in exporting information
3. The solver makes an error

First let's check option 1 $d<0$

The equation for mortality suggests this isn't possible, unless a.

• very by recompiling with a pause / print if it occurs

(BTW - the code above needed a slight tweak to run .....

out <- purrr::map_df(upper_trait_value, find_mortality_monotonic) %>%
  mutate(species1_trait = 0.07, species2_trait = upper_trait_value)


out %>%
  pivot_longer(cols = c(sp1, sp2), names_to = "species") %>%
  mutate(sp2_sp1_lma_diff = species2_trait - species1_trait) %>%
  ggplot() +
  geom_line(aes(x = sp2_sp1_lma_diff, y = value, group = species, col = species), size = 1.5) +
  xlim(-0.05, 0.15) +
  theme_classic() +
  geom_vline(xintercept = 0, col = "orange", size = 1, linetype = 2) +
  ylab("Proportion of cohorts with monotonic mortality~time curves")

)

[dfalster]

I've confirmed that mortality rate never goes negative. I added the following code into compute_rates function in cohort.h here:

 if (plant.rate("mortality") < 0) {
    std::cout << plant.rate("mortality") <<" ";
  }

Nothing prints to screen!

(but it does if I ask it to print with mortality > 0, just to check print is working)