
In recent years, the concept of "win-win" solutions for environmental conservation and profitable farming has gained traction as a seemingly ideal way to balance competing interests. The idea is that by finding mutually beneficial solutions, conservation goals can be achieved while also promoting farm profitability. But is it really possible to balance conservation and food production in a way that truly benefits both? [A recent meta-analysis](https://www.sciencedirect.com/science/article/pii/01013) spanning 43 studies across 18 countries found that the diversification of food production systems led to win-wins in only 23% of cases. In this blog post, we will examine the merits and limitations of win-win solutions for conservation and farming, going beyond the rhetoric to explore the realities of this approach.

## What are win-wins for farming and conservation?

Win-wins are situations where environmental conservation, as well as farming businesses can benefit. They are often presented as mutually beneficial solutions where one goal is not achieved at the expense of the other. These solutions aim to balance competing interests by identifying opportunities for both conservation and farming to coexist and thrive. Some examples methods that aim to improve conservation and farming outcomes are reviewed in the following.

### Agroforestry

This practice involves integrating trees into farming systems, which can improve biodiversity, soil health, water management, and carbon sequestration. Additionally, agroforestry can provide additional income for farmers through products such as fruits, nuts, and non-timber forest products. [A meta analysis](https://link.springer.com/article/10.1007/s13593-019-0589-8) has shown that agroforestry can improve ecosystem services without sacrificing crop productivity. In contrast, many farming systems are deliberately simplified, such as broad acre cropping, in order to make use of highly specialised farming equipment.

Barriers to implementing agroforestry are the significant investment of time, money, and resources to establish and maintain a more complex and diverse farming system. Additionally, agroforestry may require changes to traditional farming practices and may not be well-suited to certain types of land or crops.

### Conservation tillage

This practice involves reducing or eliminating the disturbance of soil through tillage, which can improve soil health, water retention, and biodiversity. Additionally, conservation tillage can reduce costs for farmers by reducing fuel and labour inputs. In Australia, no-till cropping began to gain popularity in the 1970s and 1980s as a way to address soil erosion and improve crop yields. Today, no-till cropping is widely adopted in Australia, and [has been widely acknowledged](https://www.tandfonline.com/doi/abs/10.1300/J064v23n03_06) to improve soil health, increase crop yields, and reduce the need for synthetic fertilisers. In contrast, conventional farming systems that rely on heavy tillage can lead to soil erosion and compaction which can negatively impact biodiversity and soil health.

Barriers to adopting conservation tillage include novel pests, such as earwigs, millipedes, slugs, and slaters, favoured by stubble retention. Thick stubble loads can also clog machinery used during the crop sowing.

### Cover crops

Planting cover crops between cash crops can improve soil health, reduce erosion, and increase biodiversity. Conventional farming systems, in contrast, often leave fields barren between cash crops which can lead to erosion, soil compaction and loss of biodiversity. Cover crops can provide additional income for farmers through products such as forage, biofuel, and compost. Cover crops have been claimed to improve soil health, increase crop yields, and reduce the need for synthetic fertilizers. However, [new research](https://onlinelibrary.wiley.com/doi/10.1111/gcb.16489) may suggest cover cropping can reduce corn and soybean yields in the US, and could lead to expanded cultivation and habitat loss to make up for the yield reductions.

Barriers to adopting cover cropping include uncertainty around improved farming outcomes, potential water loss through increased evapotranspiration, and additional management complexity surrounding multiple crops.

### Integrated pest management (IPM)

This approach to pest control involves using a combination of methods to reduce pest populations, including biological control, cultural control, and minimal use of chemical pesticides. Conventional farming systems, in contrast, are heavily reliant on chemical pesticides that harm beneficial insects, and negatively impact biodiversity. Integrated pest management can reduce costs for farmers, while improving environmental outcomes through a reduced reliance on pesticide use. [A meta analysis](https://www.mdpi.com/2075-4450/6/1/152) of 85 IPM projects spanning 24 countries found that mean yield increased across projects and crops by 40.9% (SD 72.3), while pesticide use declined to 30.7% (SD 34.9) compared with baseline. A total of 35 of 115 (30%) crop combinations resulted in a transition to zero pesticide use.

Barriers to the adoption of IPM include the increased complexity and knowledge requirements in managing diverse beneficial and pest species, as well as the uncertainty surrounding required beneficial abundances that will mitigate economic losses in production.

### Buffer strips

This practice involves planting strips of vegetation, such as grasses or wildflowers, along the edges of fields and streams. In many conventional farming situations, cropping verges are left bare or support weeds. These buffer strips can improve water quality by reducing run off and erosion, and can also provide habitats for a variety of wildlife species.

Barriers to adoption include the added complexity of managing additional crops, and the [inconsistent effects of non-crop habitats on pest and beneficial communities](https://www.pnas.org/doi/abs/10.1073/pnas.1800042115).

### Adaptive multi-paddock grazing

This practice involves moving livestock to different pastures at different times to allow for rest and recovery of the land. Conventional practices, typically involve large paddock sizes, lower herd densities, and less frequent rotation. In adaptive multi-paddock grazing, multiple paddocks are used per herd with the goal of managing grazing resources to improve ecological function. This can result in improved soil health, biodiversity, and water management. Additionally, this method of grazing can improve animal health and increase milk and meat production. A overview of the approach and its benefits are [available here](https://www.frontiersin.org/articles/10.3389/fsufs.2020.534187/full).

Barriers to adoption include the increased complexity surrounding management of small paddocks and more frequent movement of stock. There is also [significant controversy surrounding the scientific merits](https://www.hindawi.com/journals/ijbd/2014/163431/) of similar grazing programs.

### Sustainable intensification

This approach involves increasing crop yields and productivity through technologies such as precision agriculture, robotics, or genetic modification, while reducing environmental impacts such as pollution, deforestation, and biodiversity loss. As noted by [Cassman and Grassini](https://www.nature.com/articles/s41893-020-0507-8), while sustainable intensification is necessary to address environmental conservation goals, it is not sufficient because success in conserving natural habitat also requires good governance, appropriate legal frameworks for land tenure, and international agreements to ensure that progress towards sustainable intensification on existing farmland achieves desired environmental outcomes.

Barriers to wider adoption include the economic costs and knowledge required to implement technological solutions in farming.

### Organic farming

Organic farming is an approach to agriculture that avoids the use of synthetic fertilizers and pesticides widely used in conventional farming systems. Organic farming can improve environmental outcomes by reducing synthetic inputs and pollution, and through a greater reliance on biological solutions. Additionally, organic farming can improve farm profitability by reducing costs and increasing market demand for organic products. While [studies have found organic farming is frequently associated with lower yields, they tend to be more profitable due to the higher prices organic produce attracts](https://www.pnas.org/doi/10.1073/pnas.1423674112). [A recent meta study](https://www.frontiersin.org/articles/10.3389/fsufs.2019.00082/full) concluded that organic farms reliably promote biotic abundance, biotic richness, soil carbon, and profitability compared with conventional farms, but conventional farms generally produced higher yields.

Barriers to wider adoption include the lower typical yields achieved and difficulties in sourcing organic solutions to problems like pests, weeds, and diseases.

### Social licence

Social license can be a potential win-win for farming and conservation if it is effectively implemented and balanced. Social license refers to the acceptance and support of a specific activity or project by the local community and stakeholders. It is also the public acknowledgement that farming land comprises a large proportion of our environment. In the context of farming and conservation, it means that farmers and conservationists can work together to find mutually beneficial solutions that balance the needs of food production with the protection of natural resources. For example, farmers may adopt conservation practices that improve the health of natural resources and biodiversity, while conservationists may support the development of sustainable food production systems.

However, achieving a true win-win scenario through social license is not always easy, it requires a lot of work and good faith from all parties involved. It requires a deep understanding of the needs and priorities of the different stakeholders, and a willingness to negotiate and find common ground. It also requires adequate resources, such as funding and technical support, to implement the necessary changes. In addition, social license may not always be a win-win, as there may be inherent trade-offs and compromises required to find a balance between conservation and food production.

### Dam biodiversity

When dams are used as water storage for livestock, there can be risks associated with the water quality that may decrease livestock health. Encouraging dam biodiversity can help to improve water quality by supporting the growth of aquatic plants, which can help to filter pollutants and improve the overall health of the water. Vegetation can also improve the flow of water and reduce sedimentation. Together, these changes can result in increased livestock health, improved weight gain, and higher profitability [according to recent research](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0256089). Besides livestock health, dam biodiversity can help to support a wide range of aquatic and terrestrial species, which can increase the overall environmental values of the area.

Despite these benefits, encouraging dam biodiversity can be a costly process, requiring significant investments in infrastructure, equipment, and personnel. This may be a barrier for some organisations and communities that lack the resources to invest in these efforts. However, the return on investment may be positive [in many Australian situations](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0256089).

## Can farming save the environment?

Compared to farms, natural ecosystems are far more biodiverse, with richer ecological interactions and processes, and a greater resilience to environmental changes. In contrast, farms are generally designed to support the production of specific crops and domesticated animals, which reduces the overall diversity and complexity from the natural ecosystems they have replaced. While sustainable food production is worthy goal, I don't think it can (or should) replace environmental conservation efforts. Taking a monocultural and chemically sterilised crop from one species to multi-species crop with improved soil biodiversity and beneficial insects is a big improvement. However, even the most sustainable farm is highly simplified in terms of biodiversity and function when compared to a natural ecosystem.

I think the concept of win-wins may be more accurately framed as trade-offs. recently read a [manuscript by McShane et al. (2011)](https://www.sciencedirect.com/science/article/abs/pii/S0006320710001849) who wrote:

> Win–win scenarios, where both natural resources are conserved and human well-being is improved in specific places over time, have been difficult, if not almost impossible, to realize. Compromise, contestation and conflict are more often the norm. Often hard choices need to be made between different kinds of conservation, and between conservation and human well-being, and these should be explicitly acknowledged. Not to do so leads to unrealized expectations and ultimately unresolved conflict. Conservationists need to make greater efforts to more effectively consider trade-off choices between different points-of-view, determine at what levels biodiversity loss is acceptable, mitigate human costs, and broaden the decision-making process. At the same time, development that ignores the benefits that humans derive from ecosystems and natural resources will ultimately prove unsustainable. The challenge for the conservation and development community is to engage in a social process that allows for compromise and the explicit acknowledgement of risks and costs, while at the same time gaining ever more clarity and purpose regarding those things that should not be traded off.”

So can farming save the environment? I guess it depends on the kind of environment you want to save. In [another blog post](/posts/20230107-environmental-ethic), I wrote about the importance of acknowledging that different people want to conserve different things according to what they value.

Nevertheless, I think the extent to which sustainable farming can contribute to biodiversity conservation is an important issue that I will try to better address in another blog post.


Win-win solutions for conservation and farming

What are environmental conservationists really conserving?


We can visually show the difference in a Bayesian vs a frequentist approach to statistical inference using a simple coin toss example. 

Here, we flip a coin 200 times with $p=0.5$. 

```r
n <- 200
set.seed(1234)
outcomes <- cumsum(rbinom( n, 1, 0.5 ))
samples <- 1:n
p_hat = outcomes/samples
sd_p =  sqrt( 0.5 * (1 - 0.5) / samples )
```

Under a frequentist perspective, we might see that the more coin flips we make, the closer our sample estimate converges to the true proportion of 0.5 and the smaller the sd of our sample estimate $\hat{p}$.

```r
ggplot() + 
    geom_line(aes(samples, p_hat)) + 
    geom_ribbon(aes(samples, ymin=p_hat - sd_p, 
                             ymax=p_hat + sd_p), alpha=0.3) + 
    geom_hline(yintercept = 0.5, linetype=2) +
    dark_theme_bw()
```

![The sample estimate converges to the true parameter with increasing n](https://github.com/jamesmaino/20220810-bayesian-lm-and-glm-with-stan/blob/master/plots/plot1.png?raw=true)

But another way of seeing the coin is as probability density function where $p$ changes as we observe more data. The probability density starts wide, and gradually concentrates around the true value of $p$ as more coin tosses are observed. 

```r 
p <- seq(0, 1, length = n)
tibble(samples, outcomes) %>% 
    group_by(samples, outcomes) %>%
    expand_grid(p) %>%
    rowwise %>%
    mutate(log_lik = sum(log(dbinom(outcomes, samples, p)))) %>% 
    mutate(lik = exp(log_lik)) %>%
    group_by(samples) %>%
    mutate(lik = lik/sum(lik)) %>%
    ggplot() + 
        geom_line(aes(p, lik, group=samples, color='white'), alpha=0.2) + 
        guides(color = "none") + 
        geom_vline(xintercept = 0.5, linetype=2) +
        theme_bw() + 
        dark_theme_gray()
```

![Bayesian perspective on increasing observations of a coin toss](https://github.com/jamesmaino/20220810-bayesian-lm-and-glm-with-stan/blob/master/plots/plot2.png?raw=true)

The Bayesian approach was a bit more computationally intensive, as you would notice if you tried to make these plots. This is due to estimating the probability density of each model parameter. 

In practice, we do not need to estimate the probability density of all conceivable parameter values, rather there are more clever ways to prioritise searching parameter space in areas with higher relative probability. This is where approaches like Gibbs sampling and Hamiltonian Monte Carlo come into play. Have a read of [Statistical Rethinking by Richard McElreath](https://xcelab.net/rm/statistical-rethinking/) for a great intro. 

## Fitting a glm with glm

Here we use the `iris` data set bundled with R, which can be visualised as follows.

```r
data(iris)
d <- as.data.frame(iris)

p <- ggplot(d, aes(Sepal.Length, Petal.Length, color = Species)) +
    geom_point() +
    dark_theme_grey()
p
```
![Plotting the iris data](https://github.com/jamesmaino/20220810-bayesian-lm-and-glm-with-stan/blob/master/plots/plot3.png?raw=true)

First, we start with a simple linear regression but use the `glm` function with the argument `family=gaussian` to make to approach more general other model families (binomial, poisson, etc.). I find it is always handy to know multiple ways to fit a model.

```r
m1 <- glm(Petal.Length ~ Sepal.Length + Species,
    family = gaussian,
    data = d
)

m1_pred <- predict(m1, se.fit = TRUE) %>%
    as_tibble() %>%
    mutate(
        lower = fit - 1.96 * se.fit,
        upper = fit + 1.96 * se.fit
    )
p +
    geom_line(aes(y = m1_pred$fit)) +
    geom_ribbon(aes(ymin = m1_pred$lower, ymax = m1_pred$upper, fill = Species),
        alpha = 0.3,
        color = NA
    )
```

![linear model prediction](https://github.com/jamesmaino/20220810-bayesian-lm-and-glm-with-stan/blob/master/plots/plot4.png?raw=true)

```r
> summary(m1)

Call:
glm(formula = Petal.Length ~ Sepal.Length + Species, family = gaussian, 
    data = d)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.76390  -0.17875   0.00716   0.17461   0.79954  

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -1.70234    0.23013  -7.397 1.01e-11 ***
Sepal.Length       0.63211    0.04527  13.962  < 2e-16 ***
Speciesversicolor  2.21014    0.07047  31.362  < 2e-16 ***
Speciesvirginica   3.09000    0.09123  33.870  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.07984347)

    Null deviance: 464.325  on 149  degrees of freedom
Residual deviance:  11.657  on 146  degrees of freedom
AIC: 52.474

Number of Fisher Scoring iterations: 2
```

## Fitting a glm with rstan

One of the downsides of using `rstan` is the manual data preparation required for model fitting. While there are handy packages like `rstanarm` which help to reduce this overhead, it is not too difficult to start from scratch if we use a couple of tricks.

Using the `model.matrix` function is a handy way to build a covariate matrix which will save us having to the update the model defined in `.stan` and data passed to the `stan()` each time we tweak the model. 

```r
> model.matrix(m1)
    (Intercept) Sepal.Length Speciesversicolor Speciesvirginica
1             1          5.1                 0                0
2             1          4.9                 0                0
3             1          4.7                 0                0
4             1          4.6                 0                0
5             1          5.0                 0                0
//...
146           1          6.7                 0                1
147           1          6.3                 0                1
148           1          6.5                 0                1
149           1          6.2                 0                1
150           1          5.9                 0                1
```

Another handy tool is the `lookup()` function from `stan`, which can help find a suitable probability density function in `stan`. 

```r
> lookup(dnorm)
          StanFunction
374 normal_id_glm_lpmf
375 normal_id_glm_lpmf
376      normal_id_glm
379        normal_lpdf
380             normal
                                                       Arguments ReturnType
374   (vector y , matrix x, real alpha, vector beta, real sigma)       real
375 (vector y , matrix x, vector alpha, vector beta, real sigma)       real
376                                                            ~       real
379                            (reals y , reals mu, reals sigma)       real
380                                                            ~       real
```

Here we select the `normal` function and construct our `.stan` file as follows:

```
//lm_model.stan

data {
  int <lower = 0> N; 
  vector [N] y ; 
  int<lower=0> K;   
  matrix[N, K] X;  
}

parameters {
  vector [K] beta;
  real <lower=0> sigma;
}

model {
    y ~ normal_lpdf(beta*X, sigma);
}

```

We then construct a list of the data that matches the structure defined in the stan file and fit the stan model.

```r
stan_data <- list(
    N = nrow(mat),
    y = d$Petal.Length,
    K = ncol(mat),
    X = mat
)

m2 <- stan(
    file = "lm_model.stan",
    data = stan_data,
    warmup = 500,
    iter = 3000,
    chains = 4,
    cores = 4,
    thin = 1,
    seed = 123
)
```
If we check the model parameter estimates, they are very close to those estimated above by `lm`.

```r
> print(m2, probs = c(0.025, 0.975))
Inference for Stan model: lm_model.
4 chains, each with iter=3000; warmup=500; thin=1; 
post-warmup draws per chain=2500, total post-warmup draws=10000.

          mean se_mean   sd   2.5%  97.5% n_eff Rhat
beta[1]  -1.71    0.00 0.23  -2.16  -1.26  2765    1
beta[2]   0.63    0.00 0.05   0.54   0.72  2564    1
beta[3]   2.21    0.00 0.07   2.07   2.35  3282    1
beta[4]   3.09    0.00 0.09   2.90   3.27  2956    1
sigma     0.28    0.00 0.02   0.25   0.32  5616    1
lp__    -25.02    0.03 1.54 -28.81 -22.93  3724    1

Samples were drawn using NUTS(diag_e) at Fri Aug 12 10:02:49 2022.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
```

Constructing the prediction interval is a little different to above. Because we don't have direct estimates to the uncertainty in parameter estimates, it must be accessed via sampling the posterior distribution. So we take a bunch of samples and calculate the appropriate quantiles of interest (e.g. 95% credible intervals). 

```r
e <- rstan::extract(m2)
str(e)

dim(e$beta)
dim(mat)
pred_post <- e$beta %*% t(mat)
pred <- apply(pred_post, 2, function(x) quantile(x, c(0.025, 0.5, 0.975)))
m2_pred <- d %>%
    transmute(
        lower = pred[1, ],
        mu    = pred[2, ],
        upper = pred[3, ]
    )
p +
    geom_line(aes(y = m2_pred$mu)) +
    geom_ribbon(aes(ymin = m2_pred$lower, ymax = m2_pred$upper, fill = Species),
        alpha = 0.3,
        color = NA
    )
```

![prediction interval](https://github.com/jamesmaino/20220810-bayesian-lm-and-glm-with-stan/blob/master/plots/plot5.png?raw=true)

While this sounds like more work, in practice it is a flexible approach that can allow us to very quickly make inference on any derived or compound parameters. 


## Extensions to other model families (binomial model)

Thanks to a little bit of linear algebra used with `model.matrix()` and the $X\textbf{b}$ operation in stan, the code is relatively easy to extend to other model families, such as binomial models. This is because if covariates change, we really only need to update the model matrix, rather than individually specifying covariates. 

For the binomial data we will predict the plant species *setosa* from the `Sepal.Length` variable.

```r
d2 <- d %>%
    mutate(setosa = as.integer(Species == "setosa"))
```

There are just a few changes to make. 

Fit a binomial `glm`.

```r
glm2 <- glm(setosa ~ Sepal.Length, data = d2, family = binomial(link = "logit"))
```

The fitted values can be transformed back to the real scale using the inverse logit function 

```r
ilogit = function(x) exp(x)/(1+exp(x))
```

Update the outcome variable (1 = setosa; 0 otherwise).

```r
stan_data2 <- list(
    # ...
    y = d2$setosa,
    # ...
)
```

Update the stan file. 

```
// logit_model.stan 
data {
  int <lower = 0> N; 
  int<lower=0,upper=1> y[N];
  int<lower=0> K;   
  matrix[N, K] X;  
}

parameters {
  vector [K] beta;
}

model {
    y ~ bernoulli_logit(X*beta);
}
```

And that's it. 

![binomial data](https://github.com/jamesmaino/20220810-bayesian-lm-and-glm-with-stan/blob/master/plots/plot8.png?raw=true)

[Source code can be found here.](https://github.com/jamesmaino/20220810-bayesian-lm-and-glm-with-stan/)

Fitting simple Bayesian glm models with stan in R


Since around the start of 2022, I have been taking some time to learn some basic web development fundamentals to improve my ability to better communicate information and roll out small projects and web apps. I have a background in statistics and quantitative ecology and, like many ecologists, have spent most of my coding hours in the R language for scientific computation. 

After several years of tinkering around with web apps through frameworks, like R shiny, I became a bit frustrated by my lack of understanding of the underlying "magic" that made these things work. Playing around with some CSS, to change the default <text style={{color:"tomato", fontWeight:600}}>text color</text> of a widget was about as deep into the weeds I went. 

```html
<div style="color:red;"></div>
```
## Lessons from my first blog

**I stopped using it**. [In my old site](https://github.com/jamesmaino/jmaino_website), I wrote about some of the projects I had been working on and things I had learnt along the way. It was a fun experience, but I eventually stopped using it for a number of reasons, which are worth reflecting on here.

**It stopped being fun**. Blog posts became more of an after thought following something new I had learned. I am hoping in this new site, I will have a post in development while I am learning, so that the post can ride the same wave of excitement.  

**The compiler broke**. R does not have a great out-of-the-box versioning system for installed dependencies (but see [Packrat](https://rstudio.github.io/packrat/)), so after a few migrations to new computers the burden of getting the whole beast to compile finally beat me. I was using [blogdown](https://bookdown.org/yihui/blogdown/), which is essentially a wrapper for the static site generator [Hugo](https://gohugo.io/getting-started/quick-start/). 
 
**Interactive web apps** 
Perhaps most importantly, lately I have been finding that developing interactive, user-facing interactive applications has become the most satisfying component of my work. If you know the R programming language, R shiny is a great framework to get you started, but as you start remedying the limitations with your own custom css, and javascript, it does not take long until your codebase is a dog's breakfast. Also, javascript is *waaaaay* more ubiquitous than R. I found the css and js hacks were not a great way to learn solid fundamentals for web development. 

## The purpose of this website

1. A place I can rapidly publish content related to topics I am learning (i.e. [Posts](/posts))
2. A playground for exploring ideas, data sets, or new web components through short blog posts.
2. My portfolio of projects (i.e. [Projects](/projects))
3. Somewhere I can store and organise things I want to remember or access easily (i.e. [Resources](/resources)). 

## Ecological computation and science

I am planning on partitioning my content between two sites: one personal and one professional site (still under construction). The professional site will include topics surrounding quantitative and computational ecology, which I see as an area I will be continuing to work for some time. Any content that doesn't coherently fit this theme will go into the personal site, along with my general contact information and cv. So this web site will likely develop into a hodge-podge of thoughts, information, and coding experiments. 


Adventure into web development


[Part 1 of this post is available here.](./20181128-bayesian-dose-response-curves)

The dose at which some proportion of mortality occurs (e.g. 50%) has classically been a common summary statistic as it can be convenient to talk about mortality in terms of doses rather than model coefficients. This can be computed algebraically from the regression equation for each population as $logit(y) = a + b x_{dose}$, which can be solved as $x_{dose} = (log(\frac{y}{1-y}) - a)/b$. The lethal dose required for 50% mortality occurs when $y = 0.5$ or $x_{LD50} = -a/b$. 


```r
-coef(m1)[1]/coef(m1)[2]
```

This is a compound parameter of two model parameters that each have their own variance and co-variance and it possible to [estimate the distribution of this quantity analytically](https://cran.r-project.org/web/packages/modmarg/vignettes/delta-method.html), but not straightforward (at least not for me). 

## The delta method

The delta method can be used to estimate variance of $x_{LD50}$ by assuming the quantity is multivariate normally distributed. The `dose.p` function in the `MASS` package can do this for us using the following code. 


```r
library(MASS)
# we can simply run dose.p(m1, cf = c(1,3) or more verbosely
cf = c(1, 3) # The terms in the coefficient vector giving the intercept and coefficient of (log-)dose
p = 0.5 # Probabilities at which to predict the dose needed.
eta <- family(m1)$linkfun(p)
b <- coef(m1)[cf]
x.p <- (eta - b[1L])/b[2L]
names(x.p) <- paste("p = ", format(p), ":", sep = "")
pd <-  -cbind(1, x.p)/b[2L]
SE <- sqrt(((pd %*% vcov(m1)[cf, cf]) * pd) %*% c(1, 1))
```
```r
> # LD50
> exp(c(x.p - 1.96 * SE, x.p + 1.96 * SE))
[1]  7.294287 13.992988
```

## Using the posterior distribution

I find the solution to this problem using the Bayesian framework more intuitive. The parameter covariation structure is captured in the posterior distribution, so we can sample parameter values from the posterior distribution and then do some algebra (solving for $x_{LD50}$) as before and estimate the LC50 for the control population

```r
exp(c(x.p - 1.96 * SE, x.p + 1.96 * SE))
post <- extract(fit_mat)
logx50 <- -post$beta[, 1] / post$beta[, 3]
exp(mean(logx50))
```

```r
> exp(quantile(logx50, c(0.025, 0.975)))
     2.5%     97.5% 
 7.330893 13.955303 
```

In general, solving for the distribution of arbitrary compound parameters is much more straightforward in a Bayesian context.

## The ratio of two LC50 estimates

But perhaps the most compelling case for a Bayesian approach with dose responses is when it comes to estimating whether or not the LC50s of both populations are different (and then how different the responses are expressed as a resistance factor). You see many papers getting around this issue by somewhat informally inspecting the overlap of confidence intervals, which is likely to be an overly conservative test. The resistance factor is usually provided as the ratio of the LC50 values, which usually lacks a confidence interval.  

```r
# LC50 using maximum liklihood
# control
con = dose.p(m1, c(1, 3), p = 0.5)
exp(con)
t95 <- qnorm(0.025, lower.tail = FALSE)
exp(c(con - t95 * attr(con, "SE"), con + t95 * attr(con, "SE")))

# resistant
res = dose.p(m1, c(2, 4), p = 0.5)
res
exp(c(res - t95 * attr(con, 'SE'), res + t95 * attr(con, 'SE')))
```

```r
# resistance factor
> exp(res[1])/exp(con[1])
4.369501 
```

When using this test statistic, I circumvent the issue of estimating a resistance ratio by performing a log-likelihood ratio test on the full model with a simplified models that assume universal slope, or intercept parameters, but this approach is a little less direct than ideal. 


```r
m1_full = glm(cbind(dead, alive) ~ pop*log(dose), family = binomial('logit'), 
         data = d )
m1_no_pop = glm(cbind(dead, alive) ~ log(dose), family = binomial('logit'), 
         data = d )
anova(m1_full, m1_no_pop, test = 'LRT')
```

```r
> anova(m1_full, m1_no_pop, test = 'LRT')
Analysis of Deviance Table

Model 1: cbind(dead, alive) ~ pop * log(dose)
Model 2: cbind(dead, alive) ~ log(dose)
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        10     28.145                          
2        12     59.612 -2  -31.468 1.469e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```

But the issue of confidence bounds around the resistance factor remains.

Sampling from our posterior distribution makes this all look like a walk in the park.   

```r
# control LC50
exp(mean(-post$beta[, 1] / post$beta[, 3]))

# resistant LC50
exp(mean(-post$beta[, 2] / post$beta[, 4]))

# resistance factor with 95% credible intervals
x50_diff <- exp(-post$beta[, 2] / post$beta[, 4] - 
    -post$beta[, 1] / post$beta[, 3])
```

```r
> mean(x50_diff)
[1] 4.551865
```

```r
> quantile(x50_diff, c(0.025, 0.975))
    2.5%    97.5% 
2.737778 7.041024 
```








Estimating the ratio of two parameters and confidence bounds


Next.js is a web development framework that supports React components, folder based routing,   static page  rendering, and server-side requests, to name a few. I recently made my [personal website using Next.js for the first time](/posts/20220512-website-goals). 

MDX is an extension of markdown, a terse markup language for html, but extends support for file imports, exports, and jsx (e.g. for incorporating React components).

After reading [this great blog post](https://www.joshwcomeau.com/blog/how-i-built-my-blog/) by Josh Comeau, I decided I would try using mdx with with Next.js to build my personal blog. I then went of a bit of an adventure trying different packages and strategies.

## mdx rendering options in Next.js

There are several packages for incorporating mdx files into Next.js including:
- The official `@next/mdx`
- Hashicorp's `next-mdx-remote`
- Hashicorp's `next-mdx-enhanced`
- Kent C Dodd's `mdx-bundler`

See [this post by Tyler Smith](https://dev.to/tylerlwsmith/quick-comparison-of-mdx-integration-strategies-with-next-js-1kcm) for an overview of the strengths and weaknesses of each approach. 

I ended up flip-flopping between `@next/mdx` and `next-mdx-remote` but ultimately chose `@next/mdx` after a long journey that centered around building a Table of Contents like you see in the full-screen version of this page. At one point I even had my code base supporting `@next/mdx` and `next-mdx-remote` for different parts of my website, which shows they do not conflict. 

**My summary of the two**

|                     | `@next/mdx` | `next-mdx-remote`  |
| :-----------------  | :---------: | :----------------: |
|  Export content as React components           | ✓ |    |
|  Offical support from Next.js (Vercel)        | ✓ |    |
|  `.mdx` can be rendered as pages              | ✓ |    |
|  Supports `import` and `export`               | ✓ |    |
|  Automatic refresh on save                    | ✓ |    |    
|  Remote (serialised) markdown, e.g on servers |   |  ✓ |
|  Built-in frontmatter support                 |   |  ✓ |


## Some useful resources

If you want to make a blog, there are some good existing resources including [this comprehensive article](https://blog.logrocket.com/create-next-js-mdx-blog/) by Elijah Agbonze.

Another honorable mention is an mdx extension of Vercel's [official blog starter example on github](https://github.com/vercel/next.js/tree/canary/examples/blog-starter) by [John Polacek](https://github.com/johnpolacek/nextjs-mdx-blog-starter), but it is a bit out of date. 

## Making the Table of Contents

In summary, I use API requests to the server which parses the serialised front matter, extracting the titles similar to [Josh Comeau's solution](https://twitter.com/joshwcomeau/status/1334649414812659715?lang=en), but using `@next/mdx` to render mdx, which allows us to keep some useful features of .mdx files. 

The `MDXProvider` wraps the page in `_app.js` to create a context that allows the mdx files access to the supplied components. The `wrapper` component is a special component that will wrap the mdx content.

```javascript
// _app.js

import './markdown-wrapper.js'

const MDXComponents = {
  // ...
  wrapper: ({components, ...props}) => <MarkdownWrapper components={components}{...props} />,
  // ...
};

function MyApp({ Component, pageProps }) {
  return (
    <>
      <MDXProvider components={MDXComponents}>
          <Component {...pageProps} />
      </MDXProvider>
    </>
  )
}

export default MyApp
```

In `markdown-wrapper.js` I create a state for headings, make an API request to retrieve headings, and then add it to the page.

```javascript
// markdown-wrapper.js

import classes from './markdown-wrapper.module.css'
import React, { useEffect, useState } from 'react';
import {useRouter} from 'next/router';
import TableOfContents from '../../components/posts/table-of-contents';

const MarkdownWrapper = (props) => {
  const [headings, setHeadings] = useState('');

  const router = useRouter()
  
  const loadContents = async () => {
    const slug = router.asPath.replace(/\/(\w+)\//, '')
    const postsFolder = router.xxasPath.match(/\w+/)[0]
    const res = await fetch(`/api/post-toc?slug=${slug}&postsFolder=${postsFolder}`);
    const data = await res.json();
    setHeadings(data)
  };
  
  useEffect(() => {
    loadContents()
  }, [])

  return (
    <>
      <div className={classes.container}>
        {headings && <TableOfContents headings={headings}/> }
        <div id="introduction"></div>
        <div className={classes.postContainer} {...props}></div>
      </div>
    </>
  );
}
```
To make the url link to ech page sections work, I assign a lower case id to each heading based on the title name and then update the component styles

```javascript 
// _app.js

import './markdown-wrapper.js'
  const Heading2 = ({ children }) => {
    const idText = children.replace(/ /g, "_").toLowerCase();
    return <h2 id={idText} className={classes.heading2}>{children}</h2>;
};

const MDXComponents = {
  // ...
  h2: Heading2,
  wrapper: ({components, ...props}) => <MarkdownWrapper components={components}{...props} />,
  // ...
};

//...
```

Note that I manually insert the introduction in `markdown-wrapper.js` which I then hid in the styles module.

## Conclusion

To be honest, this whole activity, was a bit of a learning excercise, which was not efficient from the perspective of making a blog site. But I am looking forward to tinkering with other aspects of the page as I realise just how flexible and powerful a Next.js + mdx approach to blogging is. 

Making a Table of Contents for mdx posts with Next.js


Dose response are used in a variety of situations where the impact of exposure, dose, or concentration on some measurement endpoint is of interest. 

## Dose responses in ecology
For a given species, a simple mortality response to environmental conditions can be represented with the probabilistic outcome (death), which occurs with probability $p$. This simple process is know as a Bernoulli random variable.  

A motivating example is how a pest responds to increasing doses of a pesticide. Invertebrate [pests cause 10-20% of yield losses](https://www.cambridge.org/core/journals/journal-of-agricultural-science/article/crop-losses-to-pests/AD61661AD6D503577B3E73F2787FE7B2) in modern food systems. While cultural practices such as crop rotation and biological control through beneficial insects increasingly form a core component of effective and sustainable management, pesticides remain a widely used tool. The following data set from Maino et al. ([2018](http://www.bioone.org/doi/abs/10.1071/CP17327)) shows the mortality response of two populations to the pesticide omethoate.  

```r
library(tidyverse)
library(rstan)
```

```r
d = read_csv('data/DoseResponseMaino2018.csv', col_types = 'ccniii') %>%
  mutate(dead = dead + incapacitated) %>%
  dplyr::select(-incapacitated) %>%
  group_by(treatment, pop, dose) %>%
  summarise_all(sum) %>%
  mutate(total = dead + alive) %>%
  filter(dose != 0)
```

```R
> d                                                                           
# A tibble: 14 × 6
# Groups:   treatment, pop [2]
   treatment pop           dose alive  dead total
   <chr>     <chr>        <dbl> <int> <int> <int>
 1 omethoate control      0.007    52     0    52
 2 omethoate control      0.07     49     0    49
 3 omethoate control      0.7      53     3    56
 4 omethoate control      7        36    13    49
 5 omethoate control     70         0    51    51
 6 omethoate control    700         0    48    48
 7 omethoate control   7000         0    51    51
 8 omethoate resistant    0.007    49     0    49
 9 omethoate resistant    0.07     47     0    47
10 omethoate resistant    0.7      50     0    50
11 omethoate resistant    7        51     0    51
12 omethoate resistant   70         9    37    46
13 omethoate resistant  700         4    45    49
14 omethoate resistant 7000         0    52    52
```

## Logistic regression and the Bernoulli probability distribution
The Bernoulli mortality probability $p$ changes with the set of conditions represented by the design matrix $X$, where rows are observations and columns are covariates. It is convenient to assume $p$ follows the relationship:

$$p = \frac{1}{1 + e^{-X\boldsymbol{b}}}$$

where $b$ is a vector of coefficients. In this equation, $p$ is constrained to values between 0 and 1. The $Xb$ component is just a linear equation of the form $b_1x_1 + b_2x_2 + ... +b_nx_n$ which can be solved through rearrangement:

$$\ln \frac{p}{1 - p} = logit(p) = X\boldsymbol{b}$$

With the logit link function and the linear equation specifying the relationship of binomial parameter $p$ with its covariates, we now have the ingredients for a general linear model. For $X$ we allow a unique intercept and dose slope for each population, which is specified using the following formula notation.

```r
m1 = glm(cbind(dead, alive) ~ pop/log(dose) - 1, family = binomial('logit'), 
         data = d )
```

We can print the model matrix can be observed using using the `model.matrix` function.

```r
> model.matrix(m1)
   popcontrol popresistant popcontrol:log(dose) popresistant:log(dose)
1           1            0           -4.9618451              0.0000000
2           1            0           -2.6592600              0.0000000
3           1            0           -0.3566749              0.0000000
4           1            0            1.9459101              0.0000000
5           1            0            4.2484952              0.0000000
6           1            0            6.5510803              0.0000000
7           1            0            8.8536654              0.0000000
8           0            1            0.0000000             -4.9618451
9           0            1            0.0000000             -2.6592600
10          0            1            0.0000000             -0.3566749
11          0            1            0.0000000              1.9459101
12          0            1            0.0000000              4.2484952
13          0            1            0.0000000              6.5510803
14          0            1            0.0000000              8.8536654
```

The summary output tells us that one population called "resistant" appears to have a higher intercept.
```r
> summary(m1)

Call:
glm(formula = cbind(dead, alive) ~ pop/log(dose) - 1, family = binomial("logit"), 
    data = d)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.66626  -0.37482  -0.02094   0.29396   2.08559  

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
popcontrol              -3.7440     0.6145  -6.093 1.11e-09 ***
popresistant            -5.7053     0.9141  -6.242 4.33e-10 ***
popcontrol:log(dose)     1.6188     0.2548   6.354 2.09e-10 ***
popresistant:log(dose)   1.5064     0.2253   6.686 2.29e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 817.127  on 14  degrees of freedom
Residual deviance:  28.145  on 10  degrees of freedom
AIC: 50.208

Number of Fisher Scoring iterations: 7
```

This can be confirmed visually by plotting the model with 95% confidence intervals.

```r
m1pred = expand.grid(
  pop = c('control','resistant'),
  dose = 10^seq(-2, 4, length = 100)
)
m1pred$logdose = log(m1pred$dose)
pred = predict(m1, newdata = m1pred, type = 'response' , se.fit = TRUE)
m1pred$p = pred$fit
m1pred$se = pred$se.fit

ggplot() + 
  geom_point(data = d, aes(dose, dead/total, colour = pop)) +
  geom_line(data = m1pred, aes(dose, p, colour = pop)) +
  geom_ribbon(data = m1pred, 
              aes(dose, ymin = p - 1.96*se, max = p + 1.96*se, fill = pop), 
              alpha = 0.5) +
  mydarktheme +
  xlab('omethoate (mg/L)') +
  ggtitle('Frequentist') +
  scale_x_log10()

```

![Dose response curve using lm and logit link in R](https://github.com/jamesmaino/20181128-bayesian-dose-response-curves/blob/main/plots/plot1.png?raw=true)

We can also rearrange the data into a Bernoulli format (each row is an individual trial) using the following$^1$.

```r
db = d %>%
  rowwise() %>%
  mutate(outcome = list(rep(c(0, 1), c(alive, dead)))) %>%
  dplyr::select(-alive,-dead, -total) %>%
  unnest()

```
```R
> head(db)
# A tibble: 6 × 4
# Groups:   treatment, pop [1]
  treatment pop      dose outcome
  <chr>     <chr>   <dbl>   <dbl>
1 omethoate control 0.007       0
2 omethoate control 0.007       0
3 omethoate control 0.007       0
4 omethoate control 0.007       0
5 omethoate control 0.007       0
6 omethoate control 0.007       0
```

The model coefficients and standard errors are the same despite our trial size changing to $n = 1$.

```r
m1 = glm(outcome ~ pop/log(dose) - 1, family = binomial('logit'), 
         data = db )
```

```R
> summary(m1)

Call:
glm(formula = outcome ~ pop/log(dose) - 1, family = binomial("logit"), 
    data = db)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.89083  -0.06234  -0.00392   0.04577   2.94434  

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
popcontrol              -3.7440     0.6145  -6.093 1.11e-09 ***
popresistant            -5.7053     0.9141  -6.242 4.33e-10 ***
popcontrol:log(dose)     1.6188     0.2547   6.355 2.09e-10 ***
popresistant:log(dose)   1.5064     0.2253   6.686 2.29e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 970.41  on 700  degrees of freedom
Residual deviance: 181.42  on 696  degrees of freedom
AIC: 189.42

Number of Fisher Scoring iterations: 8
```

## Bayesian dose response model 

For anyone looking to dive into the Bayesian world I cannot overstate how good the course, and similarly titled book, **Statistical Rethinking** is by Richard McElreath. I initially used his `rethinking` package to fit these models, but more recently have moved over to `rstan`, which has a steeper learning curve but I thought it was about time I learnt how to write models in `stan`.

So let's refit the same model using a Bayesian approach. 

The data format required by `stan` is a bit more fiddly than `lm`, as the models are very flexible but carefully typed, e.g. integers, arrays, matrices, etc. This is required so that `stan` know what kind of data it will expect so that the model can be precompiled for computational efficiency. Computational efficiency is important as the Markov chain Monte Carlo methods for sampling from the probability distribution is a more intensive numerical procedure compared with maximum likelihood optimisation methods in frequentist statistics. 

```r
stan_data <- list(
    N = nrow(db),
    y = db$outcome,
    pop = ifelse(db$pop == "control", 1, 2),
    dose = db$dose
)

```
Here is the stan code saved in `model.stan` specifying the model. 

```
data {

  int <lower = 0> N; 
  array [N] int <lower = 0, upper = 1> y ;
  array [N] int <lower = 1> pop;
  vector <lower = 0> [N] dose; 

}

parameters {

  vector [2] a;
  vector [2] b;
  
}

model {
     // Prior models for the unobserved parameters
    a ~ normal(0, 10);
    b ~ normal(1, 10);

    y ~ bernoulli_logit(a[pop] + (b[pop] .* log(dose)));

}

```

```R
fit <- stan(
    file = "model.stan",
    data = stan_data,
    warmup = 500,
    iter = 3000,
    chains = 4,
    cores = 4,
    thin = 1,
    seed=123
)
```
The parameter estimates are quite similar to those estimated above. 

```
> print(fit, probs=c(0.025, 0.975))
Inference for Stan model: model-202208041918-1-68bf99.
4 chains, each with iter=3000; warmup=500; thin=1; 
post-warmup draws per chain=2500, total post-warmup draws=10000.

       mean se_mean   sd   2.5%  97.5% n_eff Rhat
a[1]  -3.92    0.01 0.64  -5.32  -2.80  3377    1
a[2]  -5.93    0.02 0.94  -8.06  -4.28  3407    1
b[1]   1.70    0.00 0.27   1.24   2.32  3321    1
b[2]   1.57    0.00 0.23   1.15   2.08  3435    1
lp__ -93.02    0.03 1.44 -96.65 -91.22  3101    1

Samples were drawn using NUTS(diag_e) at Thu Aug 04 19:18:11 2022.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
```
## The posterior distribution
One of the cool things about Bayesian approaches is that each parameter is assumed to have a probability distribution (rather than an error). What's more, the posterior probability distributions can take a range of forms, i.e., they are not necessarily normally distributed. 

```
e <- rstan::extract(fit)
a <- tibble(control = e$a[, 1], resistant = e$a[, 2]) %>%
    mutate(par = "a", id = 1:n())
b <- tibble(control = e$b[, 1], resistant = e$b[, 2]) %>%
    mutate(par = "b", id = 1:n())

dl <- bind_rows(a, b) %>%
    pivot_longer(-c(par, id))

ggplot(dl, aes(x = value, color = name, fill = name)) +
    geom_density(alpha = 0.5) +
    facet_wrap(~par, scales = "free") +
    mydarktheme
```
![Posterior distribution of dose response parameters](https://github.com/jamesmaino/20181128-bayesian-dose-response-curves/blob/main/plots/plot2.png?raw=true)

Here we can see how the posterior distributions of the parameters covary with the other parameter estimates. The density of the points corresponds to the probability distribution of the parameters. 

```R

dl %>%
    pivot_wider(names_from = par, values_from = value) %>%
    ggplot(aes(a, b, color = name)) +
    geom_point(alpha = 0.1) +
    facet_wrap(~name, scales = "free") +
    mydarktheme +
    guides(color = "none")
```
![The posterior distribution of the dose response parameters and their correlation structure](https://github.com/jamesmaino/20181128-bayesian-dose-response-curves/blob/main/plots/plot3.png?raw=true)

Retrieving model predictions from the posterior distribution can be a bit fiddly as there is no analytical expression based on parameter estimates. Rather, parameter samples from the posterior distribution must be transformed into model predictions and the quantiles must be estimated numerically. There are likely packages that can take care of this for you, but I haven't looked into it just yet. 

```R

pred <-
    expand.grid(dose = 10^seq(-2, 4, by = 0.1), pop = c(1, 2)) %>%
    as_tibble()

pred$pred <- NA
pred$predu95 <- NA
pred$predl95 <- NA
for (i in 1:nrow(pred)) {
    pop <- pred$pop[i]
    dose <- pred$dose[i]
    logit <- e$a[, pop] + e$b[, pop] * log(dose)
    p <- 1 / (1 + exp(-logit))
    qs <- quantile(p, c(0.025, 0.5, 0.975))
    pred$predl95[i] <- qs[1]
    pred$pred[i] <- qs[2]
    pred$predu95[i] <- qs[3]
}

pred %>%
    mutate(pop = c("control", "resistant")[pop]) %>%
    ggplot(aes(x = dose, y = pred, color = pop, fill = pop)) +
    geom_point(data = d, aes(dose, dead / total, color = pop)) +
    geom_line() +
    geom_ribbon(aes(ymin = predl95, ymax = predu95), alpha = 0.5) +
    mydarktheme +
    xlab("omethoate (mg/L)") +
    ggtitle("Bayesian") +
    scale_x_log10()

```

## Comparing the two model predictions
Hover your mouse over the plot below to compare the two predictions. 

<html>
  <style>
    {`
    .imageContainer  {
      cursor:pointer;
    }
    .imageContainer img:nth-child(1) {
      display:none
    }
    .imageContainer:hover img:nth-child(2) {
      display:none
    }
    .imageContainer:hover img:nth-child(1) {
      display:inline;
    }
    `}
  </style>
  <body>
    <div className = 'imageContainer'>
      <img width='1000px' src='https://github.com/jamesmaino/20181128-bayesian-dose-response-curves/blob/main/plots/plot1.png?raw=true'/>
      <img width='1000px' src='https://github.com/jamesmaino/20181128-bayesian-dose-response-curves/blob/main/plots/plot4.png?raw=true'/>
    </div>
    
  </body>
</html>

$^1$ Upon reading a bit more I found out that stan can indeed handle aggregated binomial data.

[Source code available here](https://github.com/jamesmaino/20181128-bayesian-dose-response-curves/)

[Part 2 of this post is available here.](./20181129-bayesian-dose-response-curves-part2)


Testing differences in dose responses in two populations


What if a growing population gets eaten by another population?

In a previous post I showed why we might [expect a population to grow exponentially when not resource limited](/posts/20181106-population-growth-modelling). We then extended this to the case where a [population reaches some carrying capacity](/posts/20181107-logistic-growth-modelling) (using a simple and non-mechanistic logistic function). 

But population growth can also be curtailed through interactions with another species population, such as a predator. In my area of study, we deal with a lot of herbivorous pests of agricultural systems. Despite the predominance of pesticides as a pest management strategy, there is a growing appreciation that "beneficial" arthropods, such as predators, will suppress pests. 

## Lotka-Volterra
Most ecologists will have heard the conjugation of Lotka-Volterra in reference to predator-prey dynamics, and here we will have a look at the equations that have become synonymous with their names. 

To introduce another population, we require another variable to keep track of it's state (another state variable). Let's go with $P$ for predator numbers. If our unconstrained prey population growth $\frac{dN}{dt}$  rate can be described as  $rN$ where $r$ is the difference between per capita reproduction and mortality rates, it is reasonable to suppose that a predator population will have an negative effect on $\frac{dN}{dt}$ as a consequence of the amount of prey ingested per time. The ingestion rate per predator might depend linearly on the number of prey available or $aN$ where $a$ is a predator's consumption rate of prey per available prey. Multiplying this by the total predators $aNP$ is the total prey consumed per time and prey population growth can be described as:

$\frac{dN}{dt} = rN  - aNP$

As with the prey, the predator population growth $\frac{dP}{dt}$, will be proportional to the reproduction rate minus mortality. Per capita predator reproduction can be assumed to depend on the amount of prey ingested multiplied by some efficiency constant $e$ or $eaN$, with per capita mortality simply a constant $m$ which together makes:  

$\frac{dP}{dt} = (eaN - b)P$ 

## Solving the equations numerically

As these populations depend on each other, both equations are needed to describe how the system changes through time. This is called a system of differential equations. As before, we can use the `deSolve` R package to numerically solve the situation were $r = 0.3$ (prey pop growth rate), $a = 0.1$ (predator ingestion rate per prey), $b = 0.4$ (predator mortality rate), and $e = 0.5$ (conversion efficiency of prey to predators).  

```r
library(tidyverse)
library(deSolve)

parameters <- c(r = 0.3, a = 0.1, b = 0.4, e = 0.5)

state <- c(N = 2, P = 2)

lotka_volterra <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <-  r*N - a*N*P
    dP <-  e*a*N*P - b*P
    list(c(dN, dP))
  })
}

times <- seq(0, 100, by = 0.2)

sol = as.data.frame(
  ode(y = state, times = times, func = lotka_volterra, parms = parameters)
)

```

```r
sol2 = gather(sol, species, individuals, -time)

ggplot(sol2) + 
  geom_line(aes(time, individuals, colour = species))+
   theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(.2,.21,.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  )
```
![Cyclical Lotka-Volterra predator prey](https://raw.githubusercontent.com/jamesmaino/20181111-predator-prey-dynamics/main/plots/plot1.png)

Another way of visualising the solution to the ODE system is by plotting each population on each axis.

```r
ggplot(sol) + 
  geom_path(aes(N, P), colour = 'red', size = 1)+
   theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(.2,.21,.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  )
```
![Cyclical Lotka-Volterra predator prey trace](https://raw.githubusercontent.com/jamesmaino/20181111-predator-prey-dynamics/main/plots/plot2.png)

While you may have already guessed it, this plot shows that the population densities are cyclical.   


## Incorporating logistic growth

Let's now incorporate the logistic growth function in the prey population to see the effect (assuming a carrying capacity of $k = 15$).

```r
parameters <- c(r = 0.3, a = 0.1, b = 0.4, e = 0.5, K = 15)
state <- c(N = 2, P = 2)
lotka_volterra_K <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <-  (1 - N/K)*r*N - a*N*P
    dP <-  e*a*N*P - b*P
    list(c(dN, dP))
  })
}

times <- seq(0, 100, by = 0.2)
sol = as.data.frame(
  ode(y = state, times = times, func = lotka_volterra_K, parms = parameters)
)

```

```r
soll = gather(sol, species, individuals, -time)

ggplot(soll) + 
  geom_line(aes(time, individuals, colour = species))+
   theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(.2,.21,.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  )
```

![Cyclical Lotka-Volterra predator prey with logistic growth dampening](https://raw.githubusercontent.com/jamesmaino/20181111-predator-prey-dynamics/main/plots/plot3.png)

And now plotting $N$ vs. $P$ with rainbow time.

```r
ggplot(sol) + 
  geom_path(aes(N, P, colour = time), size = 1)+
   theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(.2,.21,.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  ) + 
  scale_color_gradientn(colours = rainbow(9))
```

![Cyclical Lotka-Volterra predator prey with logistic growth and dampening trace](https://raw.githubusercontent.com/jamesmaino/20181111-predator-prey-dynamics/main/plots/plot4.png)

Wow! Now we have dampening and a stable state!

```r
> sol[nrow(sol),]
    time        N        P
501  100 8.002624 1.398892
```


At first glance, these simulations look kind of reasonable. Our prey population grows exponentially as before, but now reaches a point of decline as the predator population builds up. Population numbers can be cyclical, like in our first example, or may approach a steady state, like in our second example. 

One of the side-effects of working with continuous variables like $N$ and $P$ is that we get "fractions" of individuals, which is not realistic in the strict sense. For example, what does it mean when the predator population dips below 1? Should the population go extinct? This highlights another other issue around having assumed a deterministic system, rather than a stochastic one. 

In this example, we used a concrete case where parameters, such as reproduction, or ingestion efficiency were fixed, but there is a large literature that explores the behaviour of such systems in abstract terms. What combination of parameters while lead to particular system behaviour? For example, we might want to solve the steady state of the system in general terms. For interested readers, I can highly recommend the book by the legendary John Maynard-Smith for a clear overview of some classic models, [Models in Ecology](https://books.google.com.au/books/about/Models_in_Ecology.html?id=Jgk4AAAAIAAJ).

[Source code is available here](https://github.com/jamesmaino/20181111-predator-prey-dynamics/)

Prey population growth constrained by predators (Lotka-Volterra)


Let's derive some more population growth functions!

In [a previous post](/posts/20181106-population-growth-modelling) we derived a function for population growth based on the vital rates of reproduction and mortality. We assumed that the growth rate was constant with respect to the number of individuals in the population or $\frac{{dt}}{{dN}} = rN$. This led to the unrealistic prediction that populations will grow indefinitely. Of course, populations will eventually run into resource problems (e.g. no food) or other density dependent issues (e.g. over-crowding and disease). 

One way to set un upper limit on population growth is to scale the growth rate $rN$ by a function that is 1 when $N$ is small and 0 when $N$ is big. Here is one such function: 

$$\frac{{dN}}{{dt}} = (1 - N/K)rN$$

When $N$ is small the growth rate is approximately $rN$ but as $N \to K$ the growth rate scales to zero. Thus $K$ can be thought of as the carrying capacity of the population. It is important to note that unlike the equation for unconstrained exponential growth, this is a function of convenience rather than having mechanistic interpretation based on how reproduction or mortality depend on population density - it is the simplest function that satisfies our requirements for density dependent growth. 

## Numerical solution
Before we solve this equation with maths (analytically), let's solve it using computers (numerically) for $r=0.3$, and $K=100$ assuming the population started with just 2 individuals. 


```r
library(tidyverse)
library(deSolve)
parameters <- c(r = 0.3, K = 100)
state <- c(N = 2)
logistic <- function(t, state, parameters) {{
  with(as.list(c(state, parameters)), {{
    dN <-  (1 - N / K) * (r * N )
    list(dN)
  }})
}}
times <- seq(0, 50, by = 0.2)
sol = as.data.frame(
  ode(y = state, times = times, func = logistic, parms = parameters)
)
 
```

```r
ggplot(sol) + 
  geom_line(aes(time, N), col = 'red')+
   theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(.2,.21,.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  )
```
![Logistic equation solved numerically](https://github.com/jamesmaino/20181106-logistic-population-growth/raw/main/plots/ode1.png)

Wow! That was easy. 

## Analytical solution
Now that we've seen the pretty curve, let's use our excitement to try and describe the curve with a single equation without derivatives (solve the differential equation for $N$ as a function of $t$ and get rid of the $\frac{{dN}}{{dT}}$). 

$$\frac{{dN}}{{dt}} = (1 - N/K)rN$$

First, separate the variables. 

$${{rdt}} = \frac{{1}}{{(1 - N/K)N}} {{dN}}$$

Now we want to split up the fraction into something more manageable so we assume there exists some $A$ and $B$ such that: 

$$\frac{{1}}{{(1 - N/K)N}} = \frac{{A}}{{(1 - N/K)}} + \frac{{B}}{{N}}$$

$$\frac{{1}}{{(1 - N/K)N}} = \frac{{AN+B(1 - N/K)}}{{(1 - N/K)N}}$$

$$1 = AN + B(1 - N/K)$$

When $N=K$:

$$1 = AK + B(1 - K/K)$$

$$A = 1/K$$

Now substituting $A$ and N = 0:

$$1 = (0)/K + B(1 - (0)/K)$$

$$B = 1, A = 1/K$$

Thus we have:

$${{rdt}} = \frac{{1/K}}{{(1 - N/K)}} + \frac{{1}}{{N}} {{dN}}$$

Now integrate both sides:

$$\int {{rdt}} = \int \frac{{1/K}}{{(1 - N/K)}} + \frac{{1}}{{N}} {{dN}}$$

$$rt = \int \frac{{1/K}}{{(1 - N/K)}}dN + \int \frac{{1}}{{N}} {{dN}}$$

$$rt = \int \frac{{1/K}}{{(1 - N/K)}}dN + \ln|N|$$

The part $\int \frac{{1}}{{(K - N)}}dN$ is a little tricky, but substituting $u = 1 - N/K$ and $du = -1/K dN$ will help. 


From here we can finish the original integration.  

$$rt = \int \frac{{1/K}}{{(1 - N/K)}}dN + \ln|N| + c$$

$$rt =  -\ln|1 - N/K| + \ln|N| + c$$

As N is strictly positive we can remove the absolute operators. Remembering your logarithm rules $\ln(ab) = \ln(a) + \ln(b)$ we get:

$$rt = \ln\frac{{cN}}{{1 - N/K}}$$

$$ce^{{rt}} = \frac{{N}}{{1 - N/K}}$$

We really just want one $N$ so that we can conveniently write the function but for now let's solve the constant $c$ at $t=0$ and $N=N_0$

$$ce^{{r(0)}} = c = \frac{{N_0}}{{1 - N_0/K}} = \frac{{KN_0}}{{K - N_0}}$$

Back to solving for $N$. 

$$ce^{{rt}} = \frac{{N}}{{1 - N/K}}$$

$$(1 - N/K)ce^{{rt}} ={{N}}$$

$$ce^{{rt}} - ce^{{rt}}N/K ={{N}}$$

$$ce^{{rt}} = {{N}} + ce^{{rt}}N/K$$

$$ce^{{rt}} = N({{1}} + ce^{{rt}}/K)$$

$$N = \frac{{ce^{{rt}}}} {{({{1}} + ce^{{rt}}/K)}}$$

Substitute in the constant we solved for earlier.

$$N = \frac{{\frac{{KN_0}}{{K - N_0}}e^{{rt}} }} {{ ( {{1}} + \frac{{KN_0}}{{K - N_0}}e^{{rt}}/K)}}$$

And simplify by multiplying through by $(K - N_0)e^{{rt}}$

$$N = \frac{{KN_0}} {{ {{N_0}} + (K - N_0)e^{{-rt}} }}$$

Let plot it!

```r
r = 0.3
K = 100
N0 = 2
time <- seq(0, 50, by = 0.2)
sol = data.frame(
  time,
  N = K*N0/(N0 + (K - N0)*exp(-r*time))
)
 
```

```r
ggplot(sol) + 
  geom_line(aes(time, N), col = 'red')+
   theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(.2,.21,.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  )
```

![Logistic equation solved analytically](https://github.com/jamesmaino/20181106-logistic-population-growth/raw/main/plots/ode2.png)

All that calculus for the same line `deSolve` gave us...

[Source code.](https://github.com/jamesmaino/20181106-logistic-population-growth/)

Modelling density dependent (logistic) population growth


Let's derive some population growth functions!

Populations grow. They grow positively, if rates of reproduction > mortality, or negatively, if reproduction < mortality. For an unconstrained population of size $N$ (see [here](/posts/20181107-logistic-growth-modelling) for the contained case), this difference in per capita reproduction and mortality is referred to as the *intrinsic growth rate*, $r$ and has the units individuals per individual per time $N.N^{-1}.t^{-1}$ . The value of $r$ is a constant if the age-distribution is constant (e.g. the proportion of reproductively active individuals in the population does not change). Taking $r$ as constant, the change in population size with time is simply: 

$$\frac{dN}{dt} = rN$$ 

This is a rudimentary first-order differential equation which can be solved in a few simple steps. 

$$\int 1/N dN = \int r dt$$ 
$$ln|N| = rt + C_0$$ 
$$N = C_1 e^{rt}$$

Solving $C_1$ for $N = N_0$ at $t = t_0$:

$$N = N_0 e^{rt}$$

While the parameter $r$ hopefully makes conceptual sense, it is not straightforward to measure for a given species, as it will depend on the age structure (all juveniles means no reproduction), resource availability (food for growth), and environmental temperature (slow cold ectotherms). The age structure is assumed to be in a steady state, i.e. the relative proportion of different age classes does not change. Likewise, resource availability and other environmental conditions such as temperature are also taken as fixed. Steady state can be illustrated with a simple transition matrix model. Although there are initially 1000 individuals in the first age class (e.g. a mass hatching event triggered by rain), the age classes eventually stabilise. 

```r
library(tidyverse)

trans = matrix(
  c(0.00, 0.50, 0.00,
    0.00, 0.00, 0.50,
    4.00, 0.00, 0.10),
  ncol = 3,byrow = T
  )

rownames(trans) = c('V1','V2','V3')

colnames(trans) = c('V1','V2','V3')

```

```R
> trans
   V1  V2  V3
V1  0 0.5 0.0
V2  0 0.0 0.5
V3  4 0.0 0.1
```

```r

Tmax = 200

N = matrix(rep(0, Tmax*3), nrow = 3)

N[1,1] = 1000

for(t in 2:Tmax) N[, t] = N[, t-1] %*% trans
```

```r
# plot

ages = as.data.frame(t(N))

ages$total = rowSums(ages)

ages$t = 1:Tmax

agesl = gather(ages, age.class, value, -t, - total)

ggplot(agesl) + 
  geom_area(aes(t, value/total, fill = age.class)) +
  ylab('proportion')+
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(0.2,0.21,0.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  )

```

![Stable age distribution](https://raw.githubusercontent.com/jamesmaino/20181106-population-growth-modelling/main/plots/plot1.png)

## Regressing population vs. time 
There are a number of ways that $r$ can be measured, but they all rely on the above assumptions.

The most direct way to access $r$ is to look at how an unconstrained growing population changes through time. Let's use the numbers from our previous matrix model simulation.  

```r
df = data.frame(
      t = 0:(Tmax-1),
      N = colSums(N)
     )
```

To estimate $r$ through regression, notice that a log tranformation of the growth function results in:

$$\ln N = \ln  (N_0 e^{rt})$$
$$\ln N = \ln N_0 + rt$$
This is a linear regression with an intercept $\ln N_0$ and slope $r$ which we can recover through fitting a simple linear model.


```r

lm1 = lm(log(N) ~ t , data = df)

coef(lm1)

```

```R
> coef(lm1)
(Intercept)           t 
 6.30299077  0.03407767 
```

```r
# plot 

df$pred = predict(lm1)

ggplot(df) +
  geom_line(aes(t,N), colour = 'red', alpha = 0.5) + 
  geom_line(aes(t,exp(pred)), colour = 'grey') + 
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(.2,.21,.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  ) +
  scale_y_log10()
  
```


![Regression for retrieving intrinsic growth rate](https://raw.githubusercontent.com/jamesmaino/20181106-population-growth-modelling/main/plots/plot2.png)


Notice that the recovered intercept `exp(coef(lm1)[1]) =` `r round(exp(coef(lm1)[1]))` was actually lower than the real values of 1000. This was due to the  age distribution not being in steady state, and the population not growing properly exponential. 

Now let's try and run the same matrix model simulation but with the 1000 individuals distributed across age classes at the steady state proportions.  

```r

N1_stable = 1000*N[,Tmax]/sum(N[,Tmax])

N = matrix(rep(0, Tmax*3), nrow = 3)

N[,1] = N1_stable 

for(t in 2:Tmax) N[, t] = N[, t-1] %*% trans

df = data.frame(
      t = 0:(Tmax-1),
      N = colSums(N)
     )

lm1 = lm(log(N) ~ t , data = df)

```

```r
> coef(lm2)
(Intercept)           t 
 6.90777973  0.03388836 
```
After exponentiating the intercept we are able to recover the initial population size. 

```r
> exp(coef(lm2)[1])
(Intercept) 
   1000.024 
```

## Measuring vital rates and building a life table
Old school population demographers also built life tables for a cohort of individuals followed through life, which include direct measurements of mortality, and reproduction to calculate $r$.

Let's use an example published by Carey and Bradley$^1$ of the Atlantic Spider mite. The below table shows for each day $x$ the proportion surviving $P(x)$ and mean daily reproduction rate per individual for a synchronised cohort through time starting at egg lay.   

```r
mite = read_csv('data/AtlanticSpiderMite_CareyBradley1982.csv')
```

```r
> mite
# A tibble: 26 × 3
       x `P(x)` `m(x)`
   <dbl>  <dbl>  <dbl>
 1     0      1   0   
 2     1      1   0   
 3     2      1   0   
 4     3      1   0   
 5     4      1   0   
 6     5      1   0   
 7     6      1   0   
 8     7      1   0.37
 9     8      1   3.53
10     9      1   7.43
# … with 16 more rows
# ℹ Use `print(n = ...)` to see more rows
```
Firstly, we estimate net generational reproduction $R_0 = \sum_{x} P(x)m(x) = 42.8$ . As not all individuals will reproduce at $P(x)$, i.e. some will die, this value is scaled by $m(x)$. Now let's plot mortality and reproduction data. But how fast does a population grow by this factor?

```r
# plot 

mite$`P(x).m(x)` = mite$`P(x)`*mite$`m(x)`

dw = mite %>% gather(variable, value, -x)

ggplot(dw) +
  geom_point(aes(x,value, colour = variable), size = 2) +  
  geom_line(aes(x,value, colour = variable)) +
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = rgb(.2,.21,.27)),
    text = element_text(colour = 'grey'), 
    axis.text = element_text(colour = 'grey'), 
    panel.grid = element_line(colour = 'grey')
  )
  
```

![Mite growth data from Carey and Bradley 1982](https://raw.githubusercontent.com/jamesmaino/20181106-population-growth-modelling/main/plots/plot3.png)

An analytical estimate of r from this data relies on the Lotka equation.

$$\sum_{x=\alpha}^\beta P(x)m(x)e^{-rx} = 1$$

Here we can set $x = T$ where 

$$T = \frac{\sum_{x}xP(x)m(x)}{\sum_{x}P(x)m(x)}$$

Substituting into the Lotka equation and recognising $\sum_{x}P(x)m(x) = R_0$ we get  :


$$e^{rT} \approx R_0$$

$$r \approx \ln R_0/ T$$

Alternatively, we can [solve this equation numerically.](https://scipython.com/book/chapter-8-scipy/examples/solving-the-euler-lotka-equation/)

[Source code available here](https://github.com/jamesmaino/20181106-population-growth-modelling)

$^1$ Carey, J.R. and Bradley, J.W., 1982. Developmental rates, vital schedules, sex ratios and life tables for Tetranychus urticae, T. turkestani and T. pacificus (Acarina: Tetranychidae) on cotton. Acarologia, 23(4), pp.333-345.


All Posts

Collection of all posts written to date

Win-win solutions for conservation and farming

Read more

What are environmental conservationists really conserving?

Read more

Fitting simple Bayesian glm models with stan in R

Read more

Adventure into web development

Read more

Estimating the ratio of two parameters and confidence bounds

Read more

Making a Table of Contents for mdx posts with Next.js

Read more

Testing differences in dose responses in two populations

Read more

Prey population growth constrained by predators (Lotka-Volterra)

Read more

Modelling density dependent (logistic) population growth

Read more

Modelling unconstrained population growth rates

Read more