# Experimental design for dose-response experiments - a simulation

## Introduction

Ecotoxicologist face often the problem how to allocate resource when conducting a dose-response experiment to determine the $EC_{50}$ of species towards a substance.

*How many concentrations should I test?* *How many replicates should I run per concentration?*
*How many animals should I use per replicate?*

Often the answers are constrained by laboratory resources (e.g. How many beakers are available?) or the amount of work per replicate (e.g. How much time do I need to prepare I replicate?).

One way to address such questions is to run simulations of different experimental setups and to compare their performance. In Szöcs & Schäfer (2015) I used simulations to study the statistical power of hypotheses testing in common ecotoxicological setups. In this post, I will demonstrate how to run such simulations in the context of dose-response modelling.

## Simulation setup

First of all, we need an algorithm that generates data where we know the properties of. I will simulate data from a logistic regression with a binomial response variable

$y \sim Bin(N, p)$

Here `N`

is the total number of animals and `p`

is the probability of an animal to survive.
I simulate the probability `p`

as function of the concentration `x`

:

$logit(p) = \alpha + \beta~x$

where $\alpha$ and $\beta$ are the parameters of the relationship and the $EC_{50}$ is given by $- \frac{\alpha}{\beta}$.

## Simulation scenarios

`x`

values were evenly spaced in $0 \le x \le 1$, with varying numbers of treatments:

$n_{x} = {4,6,8,10,12}$.

I will also modify the number of replicates per treatment:

$n_{rep} = {1, 3, 5, 7, 9}$

And the number of animals per replicate:

$N = {3, 6, 9, 12, 15, 18, 21, 24, 27, 30}$.

I chose $\alpha = 5$ and $\beta = -10$, giving a LC50 of 0.5 and a decent slope.

First we need a function that generates data, according to the above described setup:

The arguments of this function are all described above, let’s see how the generated data looks like:

Here I generate data with the default set-up ($\alpha = 5$, $\beta = -10$, $n_x = 6$, $n_{rep} = 3$, $N = 10$ and $n_{sim}=1$)

We can fit a log-logistic dose-response model (DRM) to this data using the `drc`

package to determine the LC50

We see that the simulated data corresponds to the $EC_{50}$ we aimed to simulate (=0.5).

As we will have a lot of simulation scenarios ($5 \cdot 5 \cdot 10~=~250$), we need a function to analyse each generated dataset. I am only interested in the deviation from the specified $EC_{50}$ value and therefore use the absolute error (AE) from 0.5 as a benchmark. Here is a function that fits a log-logistic drm to given simulated data and returns the AE from the specified value of 0.5:

This function fits a DRM to a given dataset and returns the absolute deviation from 0.5.
I wrapped some calls into `try()`

because for some data we might fail to fit a model (because of convergence problems, complete separation etc).

Let’s see the function in action with the `sim1`

object from above:

and a cross-check:

OK, now we have the parts we need - let’s get to work. First I define the scenarios I want to run through:

Each row corresponds to a specific scenario. I will now run through all scenarios, generate 100 datasets from each and calculated the AE of each:

[Important note: This code block takes some time to run! I actually run it on an Amazon EC2 instance…]

For convenience I provide the results as a file:

The bias of each of the 100 generated datasets for each scenario is stored in `ae`

.
Next I compute the mean of the AE for each scenario and add it to the `scenarios`

data.frame:

## Results

Let’s look at the results. Because we have 3 factors possibly affecting the error in $EC_{50}$ determination I will use a faceted plot:

We see that with increasing $N$, $n_{rep}$ and $n_x$ the absolute error decreases and levels off at something around 0.02.

Having only 4 tested concentrations (left-most column) is not feasible, as the AE never reaches this low level. Also, if you have only one replicate you need to have either lots of concentrations or animals (or both). However, having lots of replicates, animals and concentrations does not add much and is possibly a waste of resources.

The displayed plot matrix can help to decide how to allocate your resources. This simulation is of course only valid in the specified scenarios but could be easily modified to other scenarios.

## References

Szöcs, E., & Schäfer, R. B. (2015). Ecotoxicology is not normal: A comparison of statistical approaches for analysis of count and proportion data in ecotoxicology. Environmental Science and Pollution Research, 22(18), 13990–13999. http://doi.org/10.1007/s11356-015-4579-3