
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/)

All Posts

Collection of all posts written to date

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