R Notebook for Lotka-Volterra competition model

Lotka-Volterra_competition.Rmd

---
title: "R Notebook"
output: html_notebook
---

```{r setup}
library(ggplot2)
library(cmdstanr)

## Lotka-Volterra competitive equation
lv_fun <- function(x, r, K, alpha) {
  dx1_dt <- r[1] * x[1] * (1 - (x[1] + alpha[1] * x[2]) / K[1])
  dx2_dt <- r[2] * x[2] * (1 - (x[2] + alpha[2] * x[1]) / K[2])
  new_x1 <- max(0, x[1] + dx1_dt)
  new_x2 <- max(0, x[2] + dx2_dt)
  c(new_x1, new_x2)
}

## function to plot two curves
plot_lv <- function(r, K, alpha, init1 = c(1, 5), init2 = c(5, 1),
                    Nt = 50, name = rep(c("1", "2"), each = Nt),
                    max_axis = 30) {
  x <- matrix(0, ncol = 2, nrow = Nt * 2)
  x[1, ] <- init1
  for (t in 2:Nt) {
    x[t, ] <- lv_fun(x = x[t - 1, ], r = c(r1, r2),
                     K = c(K1, K2), alpha = c(alpha12, alpha21))
  }
  x[Nt + 1, ] <- init2
  for (t in (Nt + 2):(Nt * 2)) {
    x[t, ] <- lv_fun(x = x[t - 1, ], r = c(r1, r2),
                     K = c(K1, K2), alpha = c(alpha12, alpha21))
  }
  
  ggplot(data.frame(name, x)) +
    geom_path(aes(x = X1, y = X2, colour = name), size = 1) +
    geom_point(aes(x = X1[1], y = X2[1]), colour = 2) +
    geom_point(aes(x = X1[Nt + 1], y = X2[Nt + 1]), colour = 3) +
    geom_segment(aes(x = K1, y = 0, xend = 0, yend = K1 / alpha12),
                 size = 0.2, colour = "blue") +
    geom_segment(aes(x = 0, y = K2, xend = K2 / alpha21, yend = 0),
                 size = 0.2, colour = "purple") +
    xlim(0, max_axis) + ylim(0, max_axis) +
    labs(x = "x1", y = "x2") +
    coord_fixed() +
    theme_bw() +
    theme(legend.position = "none")
}

## common settings
r1 <- 0.3
r2 <- 0.4
Nt <- 100
```
## Lotka-Volterra competiton model

$$
\frac{dx_1}{dt} = r_1 x_1 \left[1 - \left(\frac{x_1 + \alpha_{12} x_2}{K_1}\right)\right] \\
\frac{dx_2}{dt} = r_2 x_2 \left[1 - \left(\frac{x_2 + \alpha_{21} x_1}{K_2}\right)\right]
$$

## Example 1

$K_1 = 25, K_2 = 15, \alpha_{12} = 0.95, \alpha_{21} = 0.95$

このとき

$\frac{K_1}{\alpha_{12}} > K_2, \frac{K_2}{\alpha_{21}} < K_1$

```{r example1, echo=FALSE}
K1 <- 25
K2 <- 15
alpha12 <- 0.95
alpha21 <- 0.95

fig1 <- plot_lv(r = c(r1, r2), K = c(K1, K2), alpha = c(alpha12, alpha21),
                Nt = Nt)
print(fig1)
ggsave("fig1.png", fig1,
       width = 512, height = 512, units = "px", dpi = 96)
```

## Example 2

$K_1 = 15, K_2 = 25, \alpha_{12} = 0.95, \alpha_{21} = 0.95$

```{r example2, echo=FALSE}
K1 <- 15
K2 <- 25
alpha12 <- 0.95
alpha21 <- 0.95

fig2 <- plot_lv(r = c(r1, r2), K = c(K1, K2), alpha = c(alpha12, alpha21),
                Nt = Nt)
print(fig2)
ggsave("fig2.png", fig2,
       width = 512, height = 512, units = "px", dpi = 96)
```

## Example 3

$K_1 = 25, K_2 = 25, \alpha_{12} = 1.2, \alpha_{21} = 1.5$

```{r example3, echo=FALSE}
K1 <- 25
K2 <- 25
alpha12 <- 1.2
alpha21 <- 1.5

fig3 <- plot_lv(r = c(r1, r2), K = c(K1, K2), alpha = c(alpha12, alpha21),
                Nt = Nt)
print(fig3)
ggsave("fig3.png", fig3,
       width = 512, height = 512, units = "px", dpi = 96)
```

## Example 4

$K_1 = 15, K_2 = 20, \alpha_{12} = 0.6, \alpha_{21} = 0.9$

```{r example4, echo=FALSE}
K1 <- 15
K2 <- 20
alpha12 <- 0.6
alpha21 <- 0.9

fig4 <- plot_lv(r = c(r1, r2), K = c(K1, K2), alpha = c(alpha12, alpha21),
                Nt = Nt)
print(fig4)
ggsave("fig4.png", fig4,
       width = 512, height = 512, units = "px", dpi = 96)
```

## Fit

### Data gen.

```{r}
set.seed(1234)

r1 <- 0.3
r2 <- 0.4
K1 <- 15
K2 <- 20
alpha12 <- 0.4
alpha21 <- 0.9
Nt <- 25

init1 <- c(3, 1)
x <- matrix(0, ncol = 2, nrow = Nt)
x[1, ] <- init1
for (t in 2:Nt) {
  x[t, ] <- lv_fun(x = x[t - 1, ], r = c(r1, r2),
                   K = c(K1, K2), alpha = c(alpha12, alpha21))
}
for (k in 1:2)
  x[, k] <- rlnorm(Nt, log(x[, k]), 0.01)

p <- ggplot(data.frame(x1 = x[, 1], x2 = x[, 2])) +
  geom_point(aes(x = x1, y = x2)) +
  xlim(0, 15) + ylim(0, 15) +
  coord_fixed() +
  theme_bw()
print(p)
ggsave("fig5.png", p,
       width = 512, height = 512, units = "px", dpi = 96)
```

Stan

```{r cache=TRUE, warning=FALSE, message=FALSE}
stan_data <- list(N = Nt - 1, ts = 1:(Nt - 1), y_init = x[1, ], y = x[2:Nt, ])

model <- cmdstan_model("Lotka-Volterra_competition.stan")
fit <- model$sample(data = stan_data,
                    chains = 4, parallel_chains = 4,
                    iter_warmup = 2000, iter_sampling = 2000, refresh = 1000)
```

Summary

```{r}
fit$summary(variables = c("r", "K", "alpha", "sigma"))
```

```{r}
z_mean <- rbind(fit$summary(variables = c("z_init"))$mean,
                matrix(fit$summary(variables = c("z"))$mean, ncol = 2))
p +
  geom_path(data = data.frame(z_mean),
            mapping = aes(x = X1, y = X2),
            colour = "red")
ggsave("fig6.png",
       width = 512, height = 512, units = "px", dpi = 96)
```


## References

- Bob Carpenter (2018) Predator-Prey Population Dynamics: the Lotka-Volterra model in Stan. https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html#mechanistic-model-the-lotka-volterra-equations