Sequential Experimental Design for Functional Response Experiments

07/04/2019 ∙ by Hayden Moffat, et al. ∙ qut 0

Understanding functional response within a predator-prey dynamic is essentially the cornerstone for most quantitative ecological studies. Over the past 60 years, the methodology for modelling functional response has gradually transitioned from the classic mechanistic models to more statistically oriented models. To obtain inferences on these statistical models, a substantial number of experiments need to be conducted. The obvious disadvantages of collecting this volume of data include cost, time and the sacrificing of animals. Therefore, optimally designed experiments are useful as they may reduce the total number of experimental runs required to attain the same statistical results. In this paper, we develop the first sequential experimental design method for functional response experiments. To make inferences on the parameters in each of the statistical models we consider, we use sequential Monte Carlo, which is computationally efficient and facilitates convenient estimation of important utility functions. It provides coverage of experimental goals including parameter estimation, model discrimination as well as a combination of these. The results of our simulation study illustrate that for functional response experiments the sequential experimental design is beneficial for achieving our experimental goals. Computer code for implementing the methodology is available via https://github.com/haydenmoffat/sequential_design.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The term functional response refers to the number of prey consumed per predator as a function of prey density (Solomon, 1949)

. Predators’ feeding behaviour can be classified according to the type of the functional response. In this task, when the consumption rate increases linearly with prey density up to a threshold level at which it remains constant, we speak of type I functional response, which is exclusive to filter feeders

(Jeschke et al., 2004). In type II functional response, the consumption rate continuously decreases, whereas in type III it follows a sigmoid curve (Holling, 1959). Both type II and III functional responses reach a plateau at high prey densities and seem to prevail in nature (Jeschke et al., 2004; Sarnelle and Wilson, 2008).

Mathematical modelling and statistical analysis of functional response plays a crucial role in ecology as it enables us to gain a better understanding of the predator-prey interactions. Biological invasions, species extinction, biological control practices, as well as the management of ecosystems are strongly related to predators’ functional response (Smith and Slatkin, 1973; Papanikolaou et al., 2011, 2014; Dick et al., 2014). It has been shown that type II functional response destabilise predator-prey dynamics, whereas at low prey densities type III functional response acts as a stabilising factor (Oaten and Murdoch, 1975). Thus, describing a predator-prey system in such a quantitative manner allows for more accurate prediction and simulations.

Functional response experiments are set up so that a single predator (or multiple predators) has access to fixed numbers of prey for a given period of time. The number of prey that are attacked out of the total that the predator has access to in that given time period are recorded. In order to gain inferential information from these experiments, several trials need to be conducted. The obvious disadvantages of collecting this inordinate volume of data and conducting these experiments include cost, time and the sacrificing of animals. Consequently, optimal experimental design has become beneficial to behavioural ecologists to reduce the number of experimental trials that need to be run. The experimental design involves optimising a particular measure for an experimental purpose or goal.

The current literature for optimal experimental design for functional response models is scarce with one paper by Zhang et al. (2018). The approach of Zhang et al. (2018) only considers static designs, which requires selecting the number of prey available to the predator(s) for each experiment in the study prior to any experimentation. If there is little prior information on functional response model parameters, then optimal designs may be inefficient. Furthermore, Zhang et al. (2018) consider optimal designs for the purpose of precise frequentist parameter estimation for a single, assumed true, model. However, Papanikolaou et al. (2016) demonstrate that there can be significant uncertainty in which functional response model might be responsible for data generation. Therefore, the ability to acknowledge model structure uncertainty and the ability to use optimal design to help discriminate between the models is highly desirable.

In this paper we develop the first sequential experimental design approach for functional response experiments. Unlike the static design framework used by Zhang et al. (2018), the sequential design set-up allows practitioners to update their information about model structure and parameter values as observations are collected sequentially. In the optimal design context, this is important as this additional information can lead to more efficient design choices for future observations. Moreover, in contrast to Zhang et al. (2018), our approach explicitly accommodates uncertainty in the model structure. We consider optimal experimental design utility functions for the purpose of parameter estimation and/or model discrimination. Our sequential design methodology uses the sequential Monte Carlo (SMC) approach of Drovandi et al. (2013) and Drovandi et al. (2014), which is computationally efficient and permits convenient estimation of utility functions for parameter estimation or model discrimination. We demonstrate how the total entropy criterion of Borth (1975), a dual-purpose utility function for the goals of parameter estimation and model discrimination, is easily computed with our approach.

The rest of the paper is outlined as follows. Sections 2 and 3 provide background information regarding the paramount concepts required to understand the paper. They describe the fundamentals of modelling functional response and sequential experimental design, respectively. Section 4 outlines a simulation study that was conducted to illustrate the methodology proposed in the paper. This simulation study enables us to gain insight into examples of a predator-prey interaction while also evaluating the performance of the algorithm and of resulting designs. The paper is then concluded in Section 5 with a discussion of the simulation study results, the limitations of our approach and possible future work.

2 Background on Functional Response Models and Experiments

Among current behavioural ecologists, the mechanistic equations developed by Holling (1959) are favoured when modelling functional response for predator-prey interactions (see, for example, Beddington, 1975; Okuyama, 2012

). The preference for these models stems from their simple structure, where parameters can be easily translated to physical phenomena such as consumption rate and handling time. The simplest of Holling’s equations is often referred to as the disc equation or Holling’s type II functional response model. Holling’s type II model is given by the ordinary differential equation,

(1)

The parameters and represent the attack rate, i.e. the per capita prey consumption at low prey densities, and the handling time, i.e. the time a predator spends subduing, pursuing and eating a prey item, respectively. denotes the prey density in a given area, and is a function of . Holling’s modelling approach for type II illustrates a functional response curve where the consumption increases with prey density at a decelerating rate, until it reaches a plateau/constant consumption rate. An extension of the disc equation is Holling’s type III model (Holling, 1959). The type III functional response model describes situations in which the functional response curve forms a characteristic “S” shape. That is, the consumption rate accelerates at low prey densities, decelerates at high prey densities and then reaches a plateau/constant consumption rate. Holling’s type III model is given by

(2)

Although there are more complex prey-dependent functional response models in the current literature (see, for example, Jeschke et al., 2002), this paper will solely focus on the implementation of Holling’s type II and type III models due to their popularity and simplicity. However, our method can easily accommodate any functional response model.

The primary interest in (1) and (2) is on the parameters and . To obtain more information on a predator-prey interaction, particularly inferences on the related parameters, experimental data with varying initial prey densities are collected. Figure 1 shows an example of functional response data collected from an experiment conducted by Papanikolaou et al. (2016).

Figure 1: An example functional response dataset from Papanikolaou et al. (2016).

For the study, independent runs of the predator-prey system are conducted. The initial prey density at run is denoted by for . This variable is controlled by the experimenter. At each run, the number of prey consumed in a fixed time period (in hours), , denoted by

, is observed and used as the response variable. In this paper, we link Holling’s type II and type III models to probabilistic models to help account for uncertainty in the observational data. We link them in such a way so that solutions to mechanistic models are used to determine the expected proportions of prey eaten in the probabilistic models we consider.

We define to be the number of prey consumed/eaten for a single experiment that has a fixed time period of . Given that is usually fixed across experiments, we write for notational simplicity. We consider two possible distributions for . Since for any fixed time period , each of the

prey is either dead or alive, a binomial distribution might be a reasonable assumption. Alternatively, in the case where the data seems to indicate overdispersion, which often arises in functional response data

(Trexler et al., 1988), the beta-binomial distribution may be more appropriate to describe the distribution of . Fenlon and Faddy (2006) and Zhang et al. (2018) use the beta-binomial distribution to capture the variability in the data in a similar context.

For the case where the number of prey consumed is modelled by the binomial distribution, we have a set-up that is similar to that of Papanikolaou et al. (2016):

where

is the probability that a single prey has been consumed by time

and is the solution of the differential equation of the type II or type III model. and both implicitly depend on the model parameters and , but we do not explicitly write them here for notational convenience.

For the beta-binomially distributed case, we have a set-up that is similar to Zhang et al. (2018). The probability mass function for a single observation and the expected value are given by

(3)
(4)

respectively. In (3) and (4), and represent the two parameters of the beta-binomial distribution and is the beta function.

To link the solutions of the mechanistic equations to the beta-binomial distribution, we re-parameterise the beta-binomial distribution in terms of the expected proportion, , and over-dispersion parameter, , such that

Therefore, we have that

3 Sequential Experimental Design

3.1 General Notation

The following section outlines the general notation that is used throughout the paper. Define to be the total number of candidate models. Define

to be a random variable that indicates which model is responsible for data generation.

can take the values . Let

denote a vector of all the observations up to experiment number

and represent a vector of all the selected design points up to experiment number . The likelihood of observing for model with a set of parameters is denoted by . Denote to be our prior distribution, that is, our knowledge of the parameter for model prior to the experiment. The posterior distribution of for model after experiments is given by

where is the evidence for model and is given by the prior predictive probability of the observed data:

For the remainder of this article, will be referred to as just for simplicity. In the context of functional response models, , , refers to a particular functional response model and its corresponding parameter; for example, for a beta-binomial type II or III functional response model.

We now define the notation relevant to the sequential experimental design aspect of the paper. We denote the proposed design point as and a possible value of the response after we have observed as . Define to be a set of all the possible design points for a single observation and be a set of all the possible responses. The utility for the design point at observation and for model based on the current data is denoted . The utility for the proposed design point, , can be obtained by taking the expectation over the model and observation space. Section 3.3 outlines the specific utility functions used in this paper.

When we collect a new observation, we can easily update the posterior, assuming independence among observations, by multiplying the current posterior by the likelihood of the next observation :

(5)

3.2 Sequential Optimal Design

In this section we discuss the relevant theory necessary to understand the proposed sequential optimal experimental design algorithm. The algorithm itself is presented in Section 3.5. Sequential experimental design involves the utilisation of previously collected data in conjunction with a utility function to improve future data collection. We collect data points one-at-a-time and make an informed decision on the next design point. This myopic approach to experimental design has many advantages over static designs. Sequential experimental designs are generally more efficient in the presence of parameter and model uncertainty (see, for example, Dror and Steinberg, 2008) and involve lower-dimensional design optimisation problems at each iteration. We conducted a formal comparison of the efficiency of the static and sequential design methods for functional response experiments and the results are displayed in the appendix. The results illustrate that sequential designs outperform static designs for the experimental goals of parameter estimation, model discrimination and a combination of the two. Another benefit is that the sequential nature of the experimental design is well suited to SMC. SMC has several benefits which will be discussed later in Section 3.4.

Optimal experimental design involves selecting design points such that the experimental goals are achieved in the minimum possible number of experimental runs. In this paper we consider the experimental goals of parameter estimation, model discrimination and a combination of these. The experimental goals can be captured by different utility functions which depend on the currently collected data. Define to be the optimal design point for the next observation. We obtain the optimal design point by maximising the utility over the design space :

The utility of design point , , is determined by taking the expectation of the user-specified utility function, , over the response and model space:

(6)

The quantity

is the posterior predictive probability of a future observation and is given by

3.3 Utility Functions

Selecting a utility function that adequately captures the goals of an experiment is an integral part of optimal experimental design. Our aim is to select design points in order to increase our certainty around the “true model” upon observation of the experimental outcomes. In this section, we outline the three utility functions used in our SMC algorithm, all of which correspond to a specific experimental goal.

An important component of all the utilities is the Kullback-Leibler divergence (KLD)

(Kullback and Leibler, 1951). The KLD is an information-based measure of disparity between two distributions. In our case, the KLD represents the information gain on the true data generating process.

3.3.1 Parameter Estimation Utility

If the objective of our sequential experimental design is to maximise the precision of model parameter posterior distributions, the KLD between the current and updated posterior distributions is a highly useful utility. For model , the KLD between the current posterior, , and the posterior based on the observation and the proposed design point , , is given by

(7)

where the dependency of the current and updated posterior on is omitted for brevity. Equation (7) is simplified to,

(8)

where again the dependency of the current and updated posterior as well as the likelihood of on is omitted. The value represents the evidence at experiment number for model if we observe the response at the next design point .

The utility given in (8) will allow us to optimally design an experiment for the goal of parameter estimation for all of the candidate models. The utility for design point , , is given by substituting (8) into (6).

3.3.2 Model Discrimination Utility

Alternatively, a model discrimination utility may be of interest. In this case, we use a utility which is based upon the mutual information between the model indicator, , and the predicted observation,

. The mutual information is mathematically equivalent to the KLD between the joint distribution of

and and the product of the marginal distributions of and . This utility was initially suggested by Box and Hill (1967) and has been recently implemented by Drovandi et al. (2014). The utility for model discrimination is given by

(9)

The utility for design point , , is given by substituting (9) into (6).

3.3.3 Dual-Purpose Utility

Similar to the design problems discussed by Dette et al. (2001) and Zen and Tsai (2004), we consider a dual-purpose experimental goal which combines parameter estimation and model discrimination using the total entropy criterion (Borth, 1975). Denote to be the parameter estimation utility from (8) and denote the model discrimination utility from (9). The dual-purpose utility for design point is given by

(10)

which is purely the sum of the parameter estimation and model discrimination utilities.

Through the process of simplifying and removing terms which do not depend on the design point , we arrive at a dual-purpose utility:

The posterior predictive distribution is determined by averaging over all the models.

Given the form of these utility functions, most of these utilities are analytically intractable and therefore must be estimated. Unfortunately, estimating these quantities is not a straightforward process. SMC enables us to form particle approximations to a number of intractable integrals contained within utility functions. In addition, the form of our utility functions is convenient for estimation through SMC. Sections 3.4 and 3.5 discuss these approximations in greater depth.

3.4 Sequential Monte Carlo

SMC samplers are a useful tool for assessing parameter and model uncertainty when conducting sequential experimental design. The main advantage of using an SMC framework for sequential design is that after collecting an observation, it enables us to obtain efficient approximations to the posterior and other quantities of interest. Consequently, this allows us to efficiently obtain approximations of utility functions (see Section 3.5) as data is collected and thus allows us to easily explore the design space, , for the next optimal design point.

SMC involves traversing a distinct set of weighted samples (particles) for each of our models through a sequence of slowly evolving target distributions by iteratively conducting re-weighting, resampling and move steps. We denote the set of particles representing the target for model at experiment number to be with the corresponding weights . We denote the unnormalised and normalised weights for the particle of model at experiment number as and , respectively. In this implementation of SMC, the sequence of distributions is formed through the process of data annealing. This process involves setting up a sequence of distributions by introducing data one-at-a-time to arrive at the posterior. Given a particular model, , the sequence of targets is given by

After an observation is collected, we initially re-weight the particles to reflect the new target distribution, which in our case is the updated posterior distribution. This is conducted through the process of importance sampling where the unnormalised weights for our new target are given by

After the re-weighting process is completed and the weights are normalised, the weights tend to become more skewed. This leads to the reduction in the effective sample size (ESS). The ESS is a measure of the efficiency of a particle set and refers to the number of independent samples (of equal weight) from the target distribution that the weighted sample is worth. The ESS for model

at experiment number can be estimated by

(11)

After an observation is introduced into each model, we check the condition that the ESS is greater than some threshold, , for example

. Once the ESS drops below this threshold, it indicates that the particles are less informative than a sample of E independent draws from the target distribution. Using such a sample can lead to estimates of integrals with very high or even infinite variance. Therefore, to tackle this problem and boost the ESS back up to

, we use a resampling algorithm. Although this improves the value of the ESS, the sample will contain many duplicates. Therefore, after conducting this resampling step, a move step is required.

The purpose of a move step is to diversify the set of particles whilst maintaining invariance for the current target distribution. A major benefit of SMC is that it conveniently enables us to form a Markov Chain Monte Carlo (MCMC) proposal distribution from our particle set for an MCMC move step. We are able to compute the sample covariance matrix and use it in a multivariate normal random walk proposal. A disadvantage of using an MCMC kernel is that movement of all the particles is not guaranteed. Therefore, one iteration may not be enough to diversify the particle set. An appropriate number of times to conduct the move step was proposed by

Drovandi and Pettitt (2011) and must satisfy

(12)

The value is our pre-specified probability that the particle will move and is the probability of acceptance at the MCMC move step. This acceptance probability, , is estimated by conducting one “probing” MCMC move step for each particle in the set and determining the overall proportion of particles which move.

A useful property of this algorithm is that for each model , we can approximate the log evidence, , using the particle weights. Del Moral et al. (2006) show that we can approximate the ratio of normalising constants, , and hence the posterior predictive distribution, , for each model at the current experimental number using

Since we know that and , we are able to easily approximate in our SMC algorithm. This can be achieved by adding the logarithm of the normalising constant ratio at each experimental number to the current log evidence: .

Section 3.5 discusses a major benefit of using SMC for sequential experimental design, which is that SMC produces convenient outputs that are used for estimation of utilities and other important quantities.

3.5 Estimation of Utility Functions

SMC provides an efficient way for estimating the utility functions for proposed design points. The quantities within the utility functions are estimated solely using the particle values and the corresponding weights and are computed in the same way regardless of the model and parameters chosen for the data. We now demonstrate how we can approximate and other relevant quantities such as posterior model probabilities appearing in Algorithm 1. Define and to be the updated unnormalised and normalised weights of the particle after observing response at design , respectively. We estimate the ratio of two normalising constants and hence the posterior predictive distribution by

(13)

Using the normalised weighted samples and Monte Carlo integration, we can approximate the integral within the parameter estimation utility in (8):

We estimate the posterior predictive distribution, , by

The posterior model probabilities at experiment number , , can be estimated by normalising the evidences (see (13)). Using these approximations together with (6), we can approximate the utilities.

The SMC algorithm for optimal sequential experimental design is presented in Algorithm 1.

INPUT: Number of samples for each model, , an appropriate ESS threshold, , the model prior distributions, , and the likelihood function, .

OUTPUT: The selected design points, and the responses observed at those design points, .

1:  Draw samples from model priors, , for and for .
2:  Initialise weights, , for and for .
3:  Initialise log evidences, , for .
4:  for  to  do
5:     Select design point to maximise some given utility .
6:     Generate/collect observation at the design point .
7:     for  to  do
8:        Compute the updated unnormalised weights, , for .
9:        Update the log evidence, .
10:        Normalise the weights, , for .
11:        Compute the effective sample size, .
12:        if  then
13:           Resample particle set to obtain .
14:           Set for .
15:           Set .
16:           Determine the parameters of the MCMC proposal using the current particles, .
17:           for  to  do
18:              Conduct a one iteration move step by moving the particle with an MCMC kernel of invariant distribution .
19:           end for
20:           Calculate acceptance probability, , and hence .
21:           for  to  do
22:              Move particle with an MCMC kernel of invariant distribution iterated times.
23:           end for
24:        else
25:           Set for .
26:        end if
27:     end for
28:  end for
Algorithm 1 SMC algorithm for sequential experimental design

4 Simulation Study

To demonstrate the methods highlighted in this paper, we now conduct a simulation study. The purpose of the simulation study is to determine whether our sequential experimental design methodology and the utility functions are beneficial for our experimental goals within functional response experiments. We do this by comparing our proposed methods of selecting design points to randomly generated design points. The advantages of conducting a simulation study rather than multiple physical experiments are low cost and time efficiency.

For our simulation study, we consider 40 different settings for the true data generating process, consisting of 4 different true models each with 10 different true parameter settings. The four different models are defined in Table 1.

Mechanistic Model
Distribution of Holling’s type II Holling’s type III
beta-binomial Model 1 Model 2
binomial Model 3 Model 4

Table 1: Models used for simulation study

We explore four different methods of selecting the next design point: generate the design points randomly using a uniform distribution over the discrete design space or find the optimal design with respect to one of the three optimal design utilities defined in Section

3.3. For each of the 40 different model/parameter configurations, and for each design selection method, we run Algorithm 1 30 times independently. The total number of observations collected in each run of the SMC algorithm is . The set of true parameter settings are drawn from the posterior distribution based on the dataset from Papanikolaou et al. (2016) (we use SMC to sample the posterior here). The dataset used is the same data from Figure 1.

For this simulation study, we consider the total exposure time of prey and predator to be hours. The design points are members of the set and the responses can assume the values for . The prior distribution for the parameters is given by , and . The parameter is only required in the beta-binomial models. In addition, the prior model probabilities are equal across the models.

To assess the effect of the different design point selection methods on the precision of the parameters, we compare the distribution of the Bayesian D-posterior precision (Drovandi et al., 2013) obtained by the simulation study across the different methods. The Bayesian D-posterior precision, which is estimated by taking the determinant of the weighted sample covariance matrix of the SMC particles, is a measure that allows us to quantify the precision of model parameters. The comparison of the Bayesian D-posterior precision is conducted for each true model. Similarly, to assess the model discrimination power of the algorithm, we can compare the distribution of posterior model probabilities for each model across the different design selection methods and different true models. Section 4.1 explores these comparisons and displays results from the simulation study.

4.1 Results

We begin our analysis of the simulation study results by examining the precision of our parameters in the true model after the SMC design process. Figure 2 displays the distribution of the log Bayesian D-posterior precision at experiment number across the different design selection methods for each true model. It is apparent from Figure 2 that there is a noticeable pattern in the results across the different true models. The parameter estimation design selection method outperforms the other methods, as expected. However, the total entropy (which incorporates model discrimination and parameter estimation) only performs marginally worse overall compared with the parameter estimation design. Furthermore, we find the model discrimination design does not perform well for parameter estimation.

(a) True Model 1
(b) True Model 2
(c) True Model 3
(d) True Model 4
Figure 2: Distribution of log Bayesian D-posterior precision of the true model across design selection methods and different true models. RG, PE, MD, and TE represent randomly generated, parameter estimation, model discrimination and dual-purpose (total entropy) design selection methods, respectively.

The distributions of posterior model probabilities at experiment number for the different design selection methods and each true model are shown in Figure 3. As we would expect, the model discrimination design selection method outperforms the other methods for all true models. The dual-purpose design is only marginally less effective. It is evident from Figures 3(a) and 3(b) that for the randomly generated and parameter estimation design selection methods, the algorithm struggles to discriminate between the beta-binomial models when the true model is beta-binomial. Similarly, Figures 3(c) and 3(d) show that for the same two design selection methods, the SMC design algorithm does not effectively discriminate between the binomial models when the true model is binomial. The parameter estimation design selection method does not perform well for the goal of model discrimination.

(a) True Model 1
(b) True Model 2
(c) True Model 3
(d) True Model 4
Figure 3: Distribution of posterior model probabilities across design selection methods for different true models. RG, PE, MD, and TE represent randomly generated, parameter estimation, model discrimination and dual-purpose (total entropy) design selection methods, respectively.

In both experimental goals the total entropy utility is only marginally less effective than the single-purpose utilities. This suggests that this utility is beneficial for efficient parameter estimation and model discrimination simultaneously.

Figure 4 displays the distribution of the log Bayesian D-posterior precision for each true model at several iterations in the SMC design process. This enables us to identify the experimental number at which the gain in parameter precision becomes negligible. The results in Figure 4 are generated from running SMC using the parameter estimation design selection method. Figure 4 suggests that the increase in precision for all true models seems to decay exponentially.

(a) True Model 1
(b) True Model 2
(c) True Model 3
(d) True Model 4
Figure 4: Distribution of log Bayesian D-posterior precision across multiple iterations of the SMC design process for different true models. Results are generated from running SMC using the parameter estimation design selection method.

Figure 5 shows the posterior model probabilities after 5, 10, 15, 20 and 25 experiments have been run. The results in Figure 5 are generated using the model discrimination design selection method. Similarly to Figure 4, this enables us to identify how much gain in model discriminative ability can be achieved with increasing the sample size for each true model. It is evident that after 25 experiments, Models 3 and 4 can already be identified with high probability.

(a) True Model 1
(b) True Model 2
(c) True Model 3
(d) True Model 4
Figure 5: Distribution of posterior model probabilities across multiple iterations of the SMC design process for different true models. Results are generated from running SMC using the model discrimination design selection method.
(a) Parameter Estimation Design Selection
(b) Model Discrimination Design Selection
(c) Dual-Purpose Design Selection
Figure 6: Distribution of design points after experiments for the different true models and the different design selection methods.

The distributions of design points for the different utility functions and different true models are displayed in Figure 6. Analysing these distributions separately provides insight into where we need to place design points to gain the most information about parameters and the true model. The designs for our three proposed design selection methods are predominantly bimodal with modes at boundaries of the design space. Therefore, it appears selecting design points at the boundaries is optimal for parameter estimation and model discrimination. The designs for parameter estimation are less concentrated in the lower third of the design space than the model discrimination designs. However, the parameter estimation designs are still heavily left-skewed in the lower third. The design point distribution for the dual-purpose design selection method seems to contain features of both the parameter estimation and model discrimination design distributions. This is expected as the dual-purpose utility is formed from the other two utilities. The model discrimination designs differ quite strongly between the binomial and beta-binomial models as shown in Figure 6. The beta-binomial designs have comparatively little mass at the upper boundary (almost all mass at the lower boundary), the binomial designs have much more mass at the upper boundary.

5 Discussion

Optimally designed functional response experiments are largely advantageous for quantitative ecological studies. Optimal designs in functional response can reduce cost, improve time efficiency and prevent the sacrificing of animals. The only approach for optimal experimental design within functional response experiments (Zhang et al., 2018) is a static design which is solely reliant on the information known prior to the experiment. Therefore, current methodology does not take into account the information collected during the experiment. This can lead to inefficient designs and outcomes which are not significantly informative. Using sequential designs over static designs enables practitioners to update their information on model and parameter uncertainty as observations are collected sequentially. Appendix A demonstrates that our sequential method outperforms the previous static design approach for our simulated functional response experiment.

This paper outlines the first sequential experimental design method for functional response experiments. It utilises SMC, which is computationally efficient and allows the convenient estimation of utility functions for parameter estimation or model discrimination. The results of our simulation study illustrate that our methodology is beneficial for parameter estimation and/or model discrimination within functional response experiments. We find that our dual-purpose design selection method enables us to effectively discriminate between functional response models and acquire precise parameter estimation simultaneously.

Computational issues with our algorithm occasionally arise when fitting the incorrect model to the data. After a new observation is collected, the particle weights for each model are updated in our SMC algorithm. If we are fitting the incorrect model to the data, there can be a large difference between the posterior at experiment number and the posterior at experiment number . This results in the ESS at becoming extremely low, even to the point where the ESS is approximately 1. If this occurs, only a few unique particles remain in the sample, even after the resampling and move steps are completed. Using this sample can lead to poor estimates of utility functions and other quantities. These computational issues could possibly be avoided by constructing a sequence of target distributions between the posterior at experiment number and at experiment number .

Another obvious limitation to this myopic approach is that looking only one observation ahead is not optimal (Borth, 1975). We could possibly investigate a two-step ahead dynamic programming algorithm and determine whether it has a significant effect on results.

The ability to use optimal design to help achieve the experimental goals of parameter estimation and model discrimination is highly desirable in any application. The conducted simulation study illustrates that our sequential experimental design algorithm is a useful tool for functional response experiments. Our algorithm enables a conclusion regarding model and parameter uncertainty to be reached at a much earlier stage than would be possible with a random or static design.

Acknowledgements

The authors would like to thank Dr. Theodore Kypraios of the University of Nottingham for providing feedback on an earlier draft. Computational resources and services used in this work were provided by the HPC and Research Support Group, Queensland University of Technology, Brisbane, Australia. HM was supported by a QUT Master of Philosophy Award. MH was funded by the Austrian Science Fund (FWF): J3959-N32. NEP gratefully acknowledges financial support from the Foundation for Education and European Culture under a postdoctoral research grant.

References

  • Beddington (1975) Beddington, J. R. (1975). Mutual interference between parasites or predators and its effect on searching efficiency. The Journal of Animal Ecology, 44(1):331–340.
  • Borth (1975) Borth, D. M. (1975). A total entropy criterion for the dual problem of model discrimination and parameter estimation. Journal of the Royal Statistical Society. Series B (Methodological), 37(1):77–87.
  • Box and Hill (1967) Box, G. E. and Hill, W. J. (1967). Discrimination among mechanistic models. Technometrics, 9(1):57–71.
  • Del Moral et al. (2006) Del Moral, P., Doucet, A., and Jasra, A. (2006). Sequential Monte Carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3):411–436.
  • Dette et al. (2001) Dette, H., Franke, T., et al. (2001). Robust designs for polynomial regression by maximizing a minimum of D-and D1-efficiencies. The Annals of Statistics, 29(4):1024–1049.
  • Dick et al. (2014) Dick, J. T., Alexander, M. E., Jeschke, J. M., Ricciardi, A., MacIsaac, H. J., Robinson, T. B., Kumschick, S., Weyl, O. L., Dunn, A. M., Hatcher, M. J., et al. (2014). Advancing impact prediction and hypothesis testing in invasion ecology using a comparative functional response approach. Biological Invasions, 16(4):735–753.
  • Dror and Steinberg (2008) Dror, H. A. and Steinberg, D. M. (2008). Sequential experimental designs for generalized linear models. Journal of the American Statistical Association, 103(481):288–298.
  • Drovandi et al. (2013) Drovandi, C. C., McGree, J. M., and Pettitt, A. N. (2013). Sequential Monte Carlo for Bayesian sequentially designed experiments for discrete data. Computational Statistics & Data Analysis, 57(1):320–335.
  • Drovandi et al. (2014) Drovandi, C. C., McGree, J. M., and Pettitt, A. N. (2014). A sequential Monte Carlo algorithm to incorporate model uncertainty in Bayesian sequential design. Journal of Computational and Graphical Statistics, 23(1):3–24.
  • Drovandi and Pettitt (2011) Drovandi, C. C. and Pettitt, A. N. (2011). Estimation of parameters for macroparasite population evolution using approximate Bayesian computation. Biometrics, 67(1):225–233.
  • Fenlon and Faddy (2006) Fenlon, J. S. and Faddy, M. J. (2006). Modelling predation in functional response. Ecological Modelling, 198(1):154–162.
  • Holling (1959) Holling, C. S. (1959). Some characteristics of simple types of predation and parasitism. The Canadian Entomologist, 91(07):385–398.
  • Jeschke et al. (2002) Jeschke, J. M., Kopp, M., and Tollrian, R. (2002). Predator functional responses: discriminating between handling and digesting prey. Ecological Monographs, 72(1):95–112.
  • Jeschke et al. (2004) Jeschke, J. M., Kopp, M., and Tollrian, R. (2004). Consumer-food systems: why type I functional responses are exclusive to filter feeders. Biological Reviews, 79(2):337–349.
  • Kullback and Leibler (1951) Kullback, S. and Leibler, R. A. (1951). On information and sufficiency. The Annals of Mathematical Statistics, 22(1):79–86.
  • Oaten and Murdoch (1975) Oaten, A. and Murdoch, W. W. (1975). Functional response and stability in predator-prey systems. The American Naturalist, 109(967):289–298.
  • O’Hagan and Forster (2004) O’Hagan, A. and Forster, J. J. (2004).

    Kendall’s advanced theory of statistics, volume 2B: Bayesian inference

    , volume 2.
    Arnold.
  • Okuyama (2012) Okuyama, T. (2012). Flexible components of functional responses. Journal of Animal Ecology, 81(1):185–189.
  • Overstall et al. (2018) Overstall, A. M., McGree, J. M., and Drovandi, C. C. (2018).

    An approach for finding fully Bayesian optimal designs using normal-based approximations to loss functions.

    Statistics and Computing, 28(2):343–358.
  • Overstall and Woods (2017) Overstall, A. M. and Woods, D. C. (2017). Bayesian design of experiments using approximate coordinate exchange. Technometrics, 59(4):458–470.
  • Papanikolaou et al. (2016) Papanikolaou, N., Williams, H., Demiris, N., Preston, S., Milonas, P., and Kypraios, T. (2016). Bayesian inference and model choice for Holling’s disc equation: a case study on an insect predator-prey system. Community Ecology, 17(1):71–78.
  • Papanikolaou et al. (2011) Papanikolaou, N. E., Martinou, A. F., Kontodimas, D. C., Matsinos, Y. G., and Milonas, P. G. (2011). Functional responses of immature stages of Propylea quatuordecimpunctata (Coleoptera: Coccinellidae) to Aphis fabae (Hemiptera: Aphididae). European Journal of Entomology, 108(3).
  • Papanikolaou et al. (2014) Papanikolaou, N. E., Milonas, P. G., Demiris, N., Papachristos, D. P., and Matsinos, Y. G. (2014). Digestion limits the functional response of an Aphidophagous Coccinellid (Coleoptera: Coccinellidae). Annals of the Entomological Society of America, 107(2):468–474.
  • Sarnelle and Wilson (2008) Sarnelle, O. and Wilson, A. E. (2008). Type III functional response in Daphnia. Ecology, 89(6):1723–1732.
  • Smith and Slatkin (1973) Smith, J. M. and Slatkin, M. (1973). The stability of predator-prey systems. Ecology, 54(2):384–391.
  • Solomon (1949) Solomon, M. (1949). The natural control of animal populations. The Journal of Animal Ecology, 18(1):1–35.
  • Trexler et al. (1988) Trexler, J. C., McCulloch, C. E., and Travis, J. (1988). How can the functional response best be determined? Oecologia, 76(2):206–214.
  • Zen and Tsai (2004) Zen, M.-M. and Tsai, M.-H. (2004). Criterion-robust optimal designs for model discrimination and parameter estimation in fourier regression models. Journal of Statistical Planning and Inference, 124(2):475–487.
  • Zhang et al. (2018) Zhang, J. F., Papanikolaou, N. E., Kypraios, T., and Drovandi, C. C. (2018). Optimal experimental design for predator–prey functional response experiments. Journal of The Royal Society Interface, 15(144).

Appendix A

For completeness, we have included a quantitative comparison of the efficiency of the optimal static and sequential design methodologies for functional response experiments. This comparison considers optimal designs for the purpose of precise parameter estimation and/or model discrimination. In this section, we will briefly discuss the methodology of optimal static design and display the results of a simulation study that was conducted to compare these two methodologies.

Static Optimal Design Methodology

Bayesian optimal static design involves using prior information to select multi-dimensional designs for a particular experimental goal. Obtaining an optimal design in this way is computationally challenging as it requires maximising an analytically intractable utility function over a high-dimensional design space. To determine these designs, we use the methodology of Overstall et al. (2018) to find a multi-dimensional near-optimal static design.

General Methodology

Recall the notation used within the optimal design section of this paper. For a given experiment, the total number of observations is . The user-specified utility function for the model , the design and the response is denoted . Here we do not explicitly include in the utility function because the utilities that we define in the next subsection and in the main paper have integrated out.

The optimal static design is found by maximising the expected utility over the design space. The expectation is taken with respect to the joint distribution of all unknown quantities, which include the models and experimental responses. The expectation of the user-specified utility, and thus the utility of the multi-dimensional design, , is given by

Since the expectation of the utility function typically does not have a closed-form solution, numerical methods such as Monte Carlo estimation are required to approximate it. The Monte Carlo approximation of the expected utility is given by

where represents the number of Monte Carlo draws and are samples from the joint distribution of . Therefore, to compute this approximation for the expected utility, we only require the ability to sample from the joint distribution of and to evaluate the user-specified utility function for each sample. The evaluation of the utility functions is the focus of the next subsection.

Utility Functions

The following utility functions can be used to encapsulate the goals of the experiment (e.g. parameter estimation, model discrimination and dual-purpose). The parameter estimation utility for model is given by the KLD between the prior distribution of , , and the posterior distribution of , . The parameter estimation utility is given by

(14)

The static design utility in (14) differs from the sequential parameter estimation utility function given in the main text in (8). In sequential design, data is collected one observation at a time. Therefore, utilities in these sequential designs are based on a singular design point rather than a mutli-dimensional design. Consequently, the sequential parameter estimation utility is given by the KLD between the current (after observations are collected) and updated (after observations are collected) posterior distributions rather than the prior (before any data is collected) and posterior (after all data is collected) distributions of the parameters.

The Box and Hill (1967) utility for model discrimination is given by

(15)

The dual-purpose utility suggested by Borth (1975) is given by the sum of the parameter estimation and model discrimination utilities. In all three utilities, posterior quantities are required to evaluate the utility function. Since the posterior distribution of and will not have a neat closed form, a further approximation is required.

Approximating Posterior Distributions and Other Quantities

To make inferences on models/parameters in the Bayesian framework, the joint posterior distribution of the parameters and model is of interest. This joint distribution can be represented in a convenient form: . For more information on Bayesian inference in the presence of parameter and model uncertainty, see O’Hagan and Forster (2004). The posterior distribution of parameters is given by

(16)

where the marginal likelihood is given by

The posterior model probabilities are given by

(17)

where the inside the summation represents . As displayed in (17), to compute the posterior model probability for , we require the marginal likelihood, , which can be analytically intractable. Overstall et al. (2018) proposed the use of normal-based Laplace approximations to approximate distributions of interest for static Bayesian designs. Using this approach, the marginal likelihood can be constructed by using only the posterior mode and the Fisher information. The posterior mode for model , , can be computed by . Let denote the Fisher information for given and denote the Hessian of the negative log joint density of and . The covariance of the Laplace approximation is given by

where

The Laplace approximation to the marginal likelihood is given by

where is the number of parameters in model . In practice, we can estimate the relevant quantities required for the Laplace approximation (i.e. the posterior mode, Hessian) using a numerical optimiser. The estimate of the Hessian matrix is determined using finite differencing. Using the Laplace approximation of the marginal likelihood, the posterior model probabilities can be approximated by

where again the inside the summation represents . We can plug this approximation into the model discrimination utility, given in (15), to obtain an approximated utility.

We can approximate the posterior distribution of the model parameters given the model and the data

via a normal distribution, such that

(18)

Since we are approximating the posterior of with the multivariate normal distribution displayed in (18), we can additionally consider approximating the prior distribution with a multivariate normal distribution. The benefit of approximating both prior and posterior in this way is that there is an explicit expression for the KLD between two multivariate normal distributions. Thus we can approximate the parameter estimation utility as:

where and are the mean and covariance matrix of the prior distribution of , respectively. This significantly simplifies the calculation of the parameter estimation utility. As usual, the expected utility is estimated by Monte Carlo integration over joint samples from and . By using the utility approximations described above, this can be carried out in a computationally efficient way. The only remaining objective is to maximise the expected utility over the design space.

ACE Algorithm

To maximise the expected utility in a high dimensional design space, we use the approximate coordinate exchange (ACE) algorithm (Overstall and Woods, 2017). Very briefly, this methodology is a cyclic descent algorithm which maximises the expected utility for each design point in turn in a series of optimisation steps. A Gaussian process emulator is used to approximate the expected utility at each optimisation step in order to reduce the number of expected utility calculations required. For a more detailed explanation of the algorithm and the theory behind it, see Overstall and Woods (2017).

Simulation Study Details

The purpose of this simulation study is to illustrate the benefit of using optimal sequential design over optimal static design for functional response experiments. For this simulation study, we have a similar set-up to the one described in Section 4 of the main text. However, instead of comparing the performance of different utility functions on our experimental goals, we are comparing the performance of sequential and static optimal designs for the different utility functions. To gain an insight into how optimal and non-optimal designs differ in efficiency, randomly generated designs have been included into the comparison.

For each of the 4 models in Table 1 of the main paper, we generate 10 “true parameter settings” from the posterior distribution based on the dataset from Papanikolaou et al. (2016) (SMC was used to sample from the posterior). Therefore, we have a total of 40 different true parameter/model configurations for this study. For each of these true data generating configurations, we repeat each optimal design algorithm 30 times. The total number of observations collected in each run of these algorithms is .

For this simulation study, we place a uniform prior on the model probabilities and place log-normal priors on our parameters. More specifically, the selected priors for the parameters , and are , and , respectively. In each functional response experiment, the design points in are selected from the set and the experimental responses can take the values for . The total time that a predator has access to the prey is fixed at hours.

The Bayesian D-posterior precision is a useful tool that enables us to quantify the precision of the true model parameters. By comparing this across different optimal design methodologies, we can identify which methodology is the most efficient for parameter estimation. In a similar way, comparing the model probabilities across different optimal design methodologies enables us to compare and assess the model discrimination power of each method. The subsequent section compares the sequential and static design methodologies through the simulation study results.

Results

For each true model, Figure 7 compares the distribution of the log Bayesian D-posterior precision for the designs generated from the random, sequential and static methodologies. In this particular case, the experimental aim for the static and sequential designs is precise parameter estimation. It is apparent from Figure 7 that the sequential experimental design methodology outperforms the other methods. In addition, we can see that the two optimal designs outperform the randomly generated designs. Therefore, it is clear that both optimal experimental design methods are beneficial for parameter estimation in functional response experiments.

(a) True Model 1
(b) True Model 2
(c) True Model 3
(d) True Model 4
Figure 7: Distribution of the log Bayesian D-posterior precision for random, sequential and static design methods for each true model. RG and PE signify randomly generated and parameter estimation, respectively. The goal of the optimal designs is precise parameter estimation.

The distributions of the posterior model probabilities for randomly generated, sequential and static designs for different true models are shown in Figure 8 for the case where the experimental aim is to discriminate between the models. Similar to the results of the parameter estimation experimental goal, the sequential design outperforms the other methods for the goal of model discrimination.

(a) True Model 1
(b) True Model 2
(c) True Model 3
(d) True Model 4
Figure 8: Distributions of the posterior model probabilities for randomly generated, sequential and static design methods. Distributions are displayed and compared for each true model. RG and MD signify randomly generated and model discrimination, respectively. The goal of the optimal designs is discriminating between models.

Figure 9 and Figure 10 compare the precision of parameters and model probabilities, respectively, between randomly generated, sequential and static designs when the experimental aim is dual-purpose. We can clearly see from the figures that the sequential design outperforms the static for both components of the dual-purpose experimental goal.

(a) True Model 1
(b) True Model 2
(c) True Model 3
(d) True Model 4
Figure 9: Distributions of the log Bayesian D-posterior precision for random, sequential and static design methods for each true model. RG and TE signify randomly generated and total entropy, respectively. The goal of the optimal designs is dual-purpose (parameter estimation and model discrimination).
(a) True Model 1
(b) True Model 2
(c) True Model 3
(d) True Model 4
Figure 10: Distributions of the posterior model probabilities for randomly generated, sequential and static design methods. Distributions are displayed and compared for each true model. RG and TE signify randomly generated and total entropy, respectively. The goal of the optimal designs is dual-purpose (parameter estimation and model discrimination).