- Email: [email protected]
Bayesian methods in plant conservation biology
Biological Conservation 113 (2003) 379–387 www.elsevier.com/locate/biocon
Bayesian methods in plant conservation biology Juan M. Marı´n*, Raquel Montes Diez, David Rı´os Insua Statistics and Decision Sciences Group, Department of Experimental Sciences and Engineering, Universidad Rey Juan Carlos, Tulipa´n s/n. 28933 Mo´stoles, Madrid, Spain
Abstract We provide an introduction to Bayesian methods in conservation biology, illustrating inferences, prediction and decision making issues. After presenting the basic framework with a recovery plan evaluation problem, we illustrate more complex issues related to forecasting trends of structured populations using matrix population models and, finally, describe relevant topics in spatial and logistic regression problems. Computational and other implementation difficulties are also discussed. # 2003 Published by Elsevier Science Ltd. Keywords: Conservation biology; Bayesian analysis; Matrix population models; Recovery plans; Viability; Mixtures
1. Introduction The Bayesian framework provides a unified and coherent approach to solving inference and decision making problems under conditions of uncertainty. As such, we would expect Bayesian methods to have a great impact on conservation biology. Conservation scientists gather and analyse data to improve resource management in environmental problems. These data are usually affected by uncertainty and therefore require a statistical analysis. In environmental and ecological problems, the end use of such analysis is often to reach a decision, hence requiring uncertainty to be formally accounted for in the decision making process. Bayesian analysis seems, therefore, the appropriate approach to bioconservation statistical analysis, as it provides a coherent framework that integrates multiple sources of uncertain information (data, expert opinion) within decision making. French and Rı´os Insua (2000) provide a complete description. Many other advantages support the use of Bayesian methods in biological conservation studies. For example, prior information is frequently available, possibly from previous related studies from which our analysis might benefit, if appropriately modelled. Also, we sometimes need to combine several studies, see Osenberg et al. (1999) for a review on
* Corresponding author. Tel.: +34-916647475; fax: 916647434. E-mail address: [email protected] (J.M. Marı´n).
+34-
metanalysis issues in ecology, and the Bayesian approach provides a natural way of tackling those problems through hierarchical analysis. Spatial issues are also becoming increasingly important in ecology and recent Bayesian models provide useful tools to handle them. However powerful we might think the Bayesian approach is, we should admit that most statistical analyses in this field are still performed within the classical paradigm, although the situation seems to be gradually changing. As in other areas, controversy has permeated the eventual application of Bayesian methods in ecology—see Dennis (1996) for an example. In spite of its foundational strength, we believe that the best way to introduce and further advance these methods in conservation biology is through practice. This is precisely the goal of this paper, to provide a motivated introduction to Bayesian methods. Other reviews in the field, such as Reckhow (1990), Ellison (1996), or Wade (2000), have remained at a more conceptual level. We shall emphasise issues related to plant conservation biology. After introducing the basic Bayesian framework, we first illustrate the ideas with a simple example relating to the assessment of alternative actions in recovery plans. We then describe an alternative approach to matrix population models based on dynamic linear models. In Section 5, we illustrate a logistic model to forecast germination, comparing Bayesian and standard methods. Finally, we introduce a novel model to handle spatial issues, when spatial clustering of individuals takes place.
0006-3207/03/$ - see front matter # 2003 Published by Elsevier Science Ltd. doi:10.1016/S0006-3207(03)00124-1
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
380
2. An overview of Bayesian methods The Bayesian framework for inference and decision making problems can be easily described. Indeed, we feel that one of its strengths is, precisely, the ease with which basic concepts are put into place. Frequently one of the goals in conservation biology is to learn about one or more parameters (e.g. plant mortality rates, expected lifetimes or dispersion) which describe an ecological phenomenon of interest. Shall designate such parameters. To learn about them, we shall observe the phenomenon, collect data (e.g. counts of plants at various ages or their locations) and form the likelihood pðxjÞ which describes how data x depends (probabilistically) on the parameters , that is the distribution of data x when the parameter values are fixed at . The first novel ingredient in the Bayesian approach is to take into account another source of information about the parameters: very frequently, the conservation biologist will have access to expert information, based on past experience or related studies, which we shall model through the prior distribution pðÞ. We now need a way to combine both sources of information. This is done through Bayes’ formula. It provides the posterior distribution pðxjÞ, which is the distribution of the parameter given the data x: pðjxÞ ¼ Ð
pðÞpðxjÞ : pðÞpðxjÞd
Note that the posterior will be proportional ( /) to the numerator, which we shall represent pðjxÞ / pðÞpðxjÞ: The posterior distribution summarises all available information about the parameters and may be used to solve all inference problems of interest in conservation biology, such as point and interval estimation, hypothesis testing or prediction. For example, for point estimation, we could summarise the posterior through, for example, its mean, that is ð EðjxÞ ¼ pðjxÞd: To predict further y values of the phenomenon, we can use the predictive distribution. If we knew the value of , we would use the distribution pðyjÞ. Since there is uncertainty about , modelled through the posterior pðjxÞ, we integrate it out to get the predictive ð pðyjxÞ ¼ pðyjÞpðjxÞd: The ultimate aim of statistical research in conservation biology is usually to support decision making, for instance, to decide between several alternative actions in a recovery plan. For each action a and each future result
y, we would associate a consequence cða; yÞ. For example, if action a succeeds (y, the future result in this case), we would have incurred certain expenses, created some jobs and recovered a certain piece of land for conservation purposes. Such consequence will be evaluated with its utility uða; yÞ and we will choose the action maximising the expected utility ð uða; yÞpðyjxÞdy: Our aim now is to illustrate these ideas with several applications in plant conservation biology, showing their potentiality and their computational and modelling subtleties.
3. A simple example: evaluating recovery plans To illustrate the basic concepts, we provide examples concerning the evaluation of actions in recovery plans. We assume that a certain protocol establishes action A for a certain plant. Our aim is to evaluate its success. The parameter in this case is the probability 1 of a successful intervention through action A. Let us suppose that 9 out of 12 cases ended in success. If X1 is the number of successes in 12 trials, we have a binomial model Binð12; 1 Þ and, accordingly, we have that the probability of obtaining nine successful interventions in 12 trials given 1 is: 12 9 1 ð 1 1 Þ 3 : PrðX1 ¼ 9j1 Þ ¼ 9 As we are novel to the problem, let us view all values of 1 as being equally likely and use the uniform prior pð1 Þ ¼ 1; 1 2 ½0; 1 . Hence, applying Bayes’ formula, we get the posterior pð1 jX1 ¼ 9Þ / pð1 Þ PrðX1 ¼ 9j1 Þ / 1 19 ð1 1 Þ3 ; 1 2 ½0; 1 which corresponds to a beta distribution with parameters ¼ 10 and ¼ 4 that we represent as Beð10; 4Þ. Remember that a random variable has a Beð; Þ disÞ 1 ð1 Þ1 , 8 2 tribution, ; > 0, if pð Þ ¼ GGððþ ÞGðÞ ½0; 1 its mean being þ. Given the sequential nature of Bayes’ formula, this posterior would serve as prior for further studies. For example, suppose that in a subsequent study, three out of five implementations of our action A were successful. The new posterior would be pð1 jX1 ¼ 9; X2 ¼ 3Þ / pð1 jX1 ¼ 9ÞPrðX2 ¼ 3j1 Þ / 19 ð11 Þ3 13 ð11 Þ2 / 112 ð11 Þ5 ; 1 2 ½0; 1 corresponding to a Beð13; 6Þ distribution.
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
As mentioned, we may sometimes have expert information available about parameters. For example, let us suppose that we believe that higher values of 1 are more likely and we expect that around 80% of the times that this action is applied in a similar environment it will be successful. We could model this with a Beð4; 1Þ prior distribution. Note that, as required, its mean is 0:8 and its mode is 1. With such prior, the posteriors would be Beð13; 4Þ, after the first sample, and Beð16; 6Þ, after the second one. Fig. 1 provides the six distributions involved in this example. Those in the left column refer to the evolution of beliefs about 1 , the probability of successful implementation of action A stemming from the uniform prior, whereas those on the right stem from the Beð4; 1Þ prior. The left column shows a remarkable change in beliefs, from total uncertainty to a distribution fairly concentrated in the region [0.5,0.9]. Also note that the posterior becomes increasingly concentrated as more data become available. Similar comments hold for the right column, though we model a lot of prior information, and, therefore, the evolution of beliefs is not that remarkable. Finally, note that the bottom figures are fairly similar, showing the usual fact that, as data increase, prior relevance mitigates. The posterior distribution contains all relevant information for further inferences and predictions. The
381
following comments refer to the case that stems from uniform prior beliefs and we have observed nine successes and three failures, i.e. the current beliefs about 1 follow a Beð10; 4Þ distribution. Let us suppose we need a point estimate for the probability of success of action A. We could use the posterior mean, which is 10=ð10 þ 4Þ ¼ 0:72. Should we need an interval estimate, we could provide a region which covers 1 with posterior probability 0:9; for example, we have that pð1 2 ½0:505; 0:887 jX1 ¼ 9Þ 0:9; where we have chosen the interval so that pð1 4 0:505jX1 ¼ 9Þ 0:05 and pð1 5 0:887jX1 ¼ 9Þ 0:05. From an applied point of view, we may be more interested in computing predictive probabilities. For example, suppose that we are about to implement action A in seven more locations and we are interested in computing the probability of succeeding five or more times. We have, for each k; k ¼ 0; . . . ; 7, that ð PrðY ¼ kjX1 ¼ 9Þ ¼ PrðY ¼ kj1 Þpð1 jX1 ¼ 9Þd1 ð Gð14Þ 7 k 9 ð1 1 Þ3 d1 ¼ ð1 1 Þ7-k k 1 Gð10ÞGð4Þ 1 ð Gð14Þ 7 9þk ð1 1 Þ10k d1 ¼ k Gð10ÞGð4Þ 1 7 13!ð9 þ kÞ!ð10 kÞ! ¼ k 9!3!20!
Fig. 1. Distributions involved in the recovery plan evaluation problem. 1 : probability of a successful intervention through action A.
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
382
Therefore, PrðY 5 5jX1 ¼ 9Þ ¼
7 X PrðY ¼ kjX1 ¼ 9Þ; 0:6641: k¼5
The comparison of alternative actions follows a similar path. Suppose, for example, that another action B has been used in a similar context, resulting in 10 applications with Z ¼ 6 successes. We compare 1 , the probability of success with action A, with 2 , the probability of success with action B. The model we consider is Xj1 Binð12; 1 Þ Zj2 Binð10; 2 Þ 1 ; 2 Unif ½0; 1 ; independent: For comparison purposes, we compute r ¼ Prð1 5 2 jX1 ¼ 9; Z ¼ 6Þ. To do so, we see that, a posteriori, 1 jX1 ¼ 9 Beð10; 4Þ, 2 jZ ¼ 6 Beð7; 5Þ and that they are independent. The posterior distribution of 1 2 cannot be obtained in closed form, but we may easily simulate it. For example, we may draw 1000 values from the posterior distributions of 1 and 2 , obtain their differences 1 2 , count the number of such differences which are non-negative and divide by 1000 to approximate r. In our example, we undertook such a procedure to obtain r 0:772. This result suggests that action A seems better than action B, as we have a better chance of success with the former. However, to formally make a decision, we should consider the utility of various consequences. We shall use a simplified utility function here; for further details see French and Rı´os Insua (2000). Suppose that action A implies a considerable cost, but action B implies virtually no cost. Let us also suppose that the utility function is as described in Table 1 where, for example, if we implement action A successfully, the utility obtained is 0.8. Then, conditional on 1 and 2 , the expected utilities of the actions are, respectively, 0:81 þ 0ð1 1 Þ ¼ 0:81 and 2 þ 0:2ð1 2 Þ ¼ 0:2 þ 0:82 , whereas their posterior expected utilities are, respectively, 0:8Eð1 jX1 ¼ 9Þ ¼ 0:8
10 4 ¼ ; 14 7
and Table 1 Utilities for actions A and B in a recovery plan
Action A Action B
7 2 ¼ : 12 3 Therefore, we would actually recommend plan B, since its expected utility is higher.
0:2 þ 0:8Eð2 jZ ¼ 6Þ ¼ 0:2 þ 0:8
Succeeds
Does not succeed
0.8 1
0 0.2
4. A complex example: forecasting population viability One of the most important issues in conservation biology is the evaluation of a population’s stability in terms of its size at each stage of its biological cycle. The final goal is usually to predict whether a population is decreasing, increasing or is stabilised, hence assessing its viability. To do this, a demographic analysis and a methodology to study population growth trajectories and identify life history stages are required. Many authors, for example, Schemske et al. (1994) and Caswell (2001), have advocated the use of matrix population models in the context of conservation biology. Individuals are classified into several categories that reflect natural stages of the life cycle. Fig. 2 shows a hypothetical graph of a plant population model with five states: seed, seedling, juvenile, subadult and adult, with the corresponding valid transitions, pij probabilities and Fj per capita fertilities at each stage. For example, p12 describes the probability of a seed becoming a seedling after one year, and F5 is the number of seeds produced by an adult per year. As usual in conservation biology, we do not show a death stage either in the life cycle or in the transition matrix in Table 2 where there are implicit transitions to the death stage with the corresponding probability. For example, the probability of seed mortality would be 1 p11 p12 . A Markov chain model is implicitly assumed, transition probabilities are estimated and the evolution of the population size can be studied, for instance, to compare different management strategies. Thus, this procedure can be used to estimate seeding needs or the expected number of adults in a period of time. To characterise the process, we shall consider the probabilities given in the associated transition matrix in Table 2. The parameters defining our matrix population model process are, therefore, y=(p11,p12,p22,p33,p34,p35, p44,p45,p53,p54,p55,F5). We assume they are, a priori, independent. Inference procedures for such parameters follow. For the p parameters corresponding to each column, that is, the fates of a particular life stage after a time period, we shall assume a Dirichlet-multinomial model (see O’Hagan, 1994). This model, in which more than two classes are used, is a straightforward generalisation of the b-binomial model applied in Section 3. Let us suppose that there are three classes and let X ¼ ðX1 ; X2 ; X3 Þ designate the number of occurrences in classes 1, 2 and 3 after n trials, with X1 þ X2 þ X3 ¼ n. Then, we assume Xj Multðn; 1 ; 2 Þ, where i is the
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
383
Fig. 2. Transitions from different states of the life cycle of a plant. Table 2 Associated population transition matrix
Seed Seedling Juvenile Subadult Adult
Seed
Seedling
Juvenile
Subadult
Adult
p11 p12 0 0 0
0 p22 p23 0 0
0 0 p33 p34 p35
0 0 0 p44 p45
F5 0 p53 p54 p55
probability of an observation in class i, i ¼ 1; 2 and 3 ¼ 1 1 2 is the probability of an observation in class 3. A priori, Dirð1 ; 2 ; 3 Þ, a Dirichlet distribution with parameters 1 ; 2 and 3 . Remember that the Dirichlet distribution is a multivariate generalisation of the Beta distribution. For example, the Dirð1 ; 2 ; 3 Þ 1 1 2 1 3 1 density where 0 < i < 1 P is pð1 ; 2 ; 3 Þ / 1 2 3 and i i ¼ 1. A posteriori, we see that jð x1 ; x2 ; x3 Þ Dirð1 þ x1 ; 2 þ x2 ; 3 þ x3 Þ;
which 800 germinate, 100 die and 100 remain dormant. From the 800 seedlings, 600 pass to juveniles, 180 die and the rest remain in that stage, and so on until we complete a whole cycle. For the first column, we would have that, a posteriori, ð1 p11 p12 ; p11 ; p12 Þ Dirð100 þ 1; 100 þ 1; 800 þ 1Þ; with, respectively, posterior means 100þ1 100þ1 800þ1 ¼ 0:1007; ¼ 0:1007; ¼ 0:7986 1000þ3 1000þ3 1000þ3 and very small posterior variances, namely,
9:0197 105 ; 9:0197 105 ; 1:6019 104 : The whole process is completed resulting in very concentrated distributions around the posterior means, with matrix
with posterior means i þ xi ; i ¼ 1; 2; 3: Eði jxÞ ¼ P j þ n j
If there is little information from experts, we use noninformative priors. We start with 1000 seeds, out of
0
p11 B p12 B P¼B B0 @0 0
0 p22 p23 0 0
0 0 p33 p34 p35
0 0 0 p44 p45
1 0 0 C C p53 C C p54 A p55
384
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
The above transition matrix P is conditional on the parameter values, that is, the entries are actually pij j ¼ PrðXnþ1 ¼ jjXn ¼ i; Þ. From a Bayesian point of view, the transition matrix that we should consider is the corresponding marginal distribution. However, as our posteriors are very concentrated (see Marı´n et al., 2002b), we may approximate pij jdata by Pr Xnþ1 ¼ jjXn ¼ i; ^ ¼ pij j^ , where ^ are our posterior point estimates. F5 may be similarly estimated by counting the number of seeds produced by several adults. We could then use the methods in Caswell (2001) to undertake forecasting tasks. Other cases with more uncertainty in the point estimates require simulation, as we shall illustrate in a forthcoming paper. A shortcoming of such models is that the involved parameters are fixed through time. We describe here a Dynamic Linear Model that permits dynamic variation on such parameters. Our aim is to forecast the population structure several periods (years) ahead, by using the available population structure history. a Let us suppose that nst ; nslt ; njt ; nsa t and nt are, respectively, the number of seeds, seedlings, juveniles, subadults and adults at a given year t. We model the underlying population evolution in terms of an (autoregressive) multivariate dynamic linear model, with observation equations nstþ1 ¼ pt11 nst þ F5t nat þ "ts nsltþ1 ¼ pt12 nst þ pt22 nslt þ "tsl njtþ1 ¼ pt23 nst þ pt33 njt þ pt53 nat þ "tj t j t sa t a t nsa tþ1 ¼ p34 nt þ p44 nt þ p54 nt þ "sa t a t natþ1 ¼ pt35 njt þ pt45 nsa t þ p55 nt þ "a
where ptij ; i; j ¼ 1; . . . ; 5 are the probabilities of changing (surviving, growing,. . .) from stage i to stage j and F t;5 is the fertility per capita parameter, within period t. For example, the first equation describes that the number of seeds ns;tþ1 in year t þ 1 depends on the number of seeds pt11 nst , remaining at such stage from the previous year, and seeds F5t nat coming from adults in the previous year. Similarly, the last equation suggests that the number of adults natþ1 in year t þ 1 depends on the numbers of j adults nat , subadults nsa t and juveniles nt in the previous year. The error terms are assumed to be jointly independent and normal with zero mean and unknown variance. Note that if we drop these terms, we recover a standard matrix population model with time dependent parameters. We now assume that the parameters will vary randomly with time according to the system equations ptþ1 ¼ ptij þ tij ij F5tþ1 ¼ F5t þ t5
where the pij ’s satisfy the constraints ptij 5 0 and P5 t j¼1 pij ¼ 1, 8i, and the terms are also independent and normally distributed. This model, in fact, generalises classical matrix population models by considering random evolving parameters. As the parameters of the model may evolve with time, the model includes a system of equations to upgrade their values. Following the Bayesian paradigm, we would assign a prior distribution to parameters pij and F5 at time 0. The posterior distribution of parameters is updated at each stage, also giving the predictive distributions of the variables of interest several stages ahead. This class of models has been used in many other fields within engineering and applied sciences. Its main advantage is its flexibility to adapt to many situations, where the calculations are automatically performed by using the Kalman Filter procedure, a set of embedded linear equations that resume the main stages of estimations of the model. A full description of this technique is outside the scope of this paper, but may be seen, for example, in West and Harrison (1997).
5. Logistic regression Logistic regression is a popular technique in life sciences. It can be regarded as an extension of standard regression models. Its main aim is to predict the probability of success within a given trial in relation to a group of predictive variables. For example, in a simple logistic (binomial) regression, we may have independent observations from I populations. We assume that each one is distributed according to a binomial law: yi Binðni ; i Þ, for i ¼ 1; . . . ; I, where Eðyi Þ ¼ ni i , and that the explanatory variable x;i is available. Odds are defined as i =ð1 i Þ. The logistic regression specifies a linear structure for the logarithm of the odds (log odds): i logitði Þ log ¼ 0 þ 1 xi ; i ¼ 1; . . . ; I: 1-i Note that by varying xi , the term on the right of equation can take any real value, but i is always restricted between 0 and 1. Standard methods for analysing binomial regression data rely on asymptotic inference, requiring a great amount of data to provide accurate estimates. This problem is straightforwardly worked out by Bayesian methodology, allowing us to include prior information, use diagnostic studies, undertake model selection among different models and perform sensitivity analysis. This has been facilitated by the recent introduction of Markov chain Monte Carlo methods, which in our model permit computations for prediction and making inferences on regression coefficients.
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
To fix ideas, suppose we have data yi ; x0i ; i ¼ 1; . . . n, available where xi s are known vectors of k covariates and the yi s are independent binomial random variables with ni trials. The success probability for any single trial y with covariate x is Fðx0 Þ, i.e. 0 Fðx Þ Prðy ¼ 1jx; Þ, where is an unknown k-dimensional vector of regression coefficients, and F is Fðx0 Þ ¼
expðx0 Þ : 1 þ expðx0 Þ
The likelihood function for the complete data Y ¼ ðy1 ; . . . yn Þ is
385
Table 3 Proportion of germinated seeds of two varieties of Orobanche aegyptiaca in presence of bean and cucumber root extracts. (Crowder, 1978) O. aegyptiaca 75 Bean
O. aegyptiaca 73 Cucumber
Bean
Cucumber
y
n
y/n
y
n
y/n
y
n
y/n
y
n
y/n
10 23 23 26 17 –
39 62 81 51 39 –
0.26 0.37 0.28 0.51 0.44 –
5 53 55 32 46 10
6 74 72 51 79 13
0.83 0.72 0.76 0.63 0.58 0.77
8 10 8 23 0 –
16 30 28 45 4 –
0.50 0.33 0.29 0.51 0.0 –
3 22 15 32 3 –
12 41 30 51 7 –
0.25 0.54 0.50 0.63 0.43 –
n Y n y ni 0 yi 1 F x0i i i F xi LðjYÞ yi
y: number of germinated seeds in each plate; n: total number of seeds in each plate.
For a given prior distribution on , namely ðÞ, to obtain posterior and predictive distributions we need to compute the posterior distribution of :
yi Binðni ; i Þ
i¼1
LðjYÞ ðÞ : LðjYÞ ðÞd As many important quantities are obtained through cumbersome and intractable integrals involving ðjYÞ, approximations must be undertaken. Monte Carlo methods provide discrete approximations to the posterior distribution. Given a function hðÞ, the posterior expectation E½hðÞjY may be approximated by
ðjYÞ ¼ Ð
ð t
1X E½hðÞjY ¼ hðÞ ðjYÞd h i ; t i¼1 t for a sample i i¼1 of size t, from the posterior distribution ðjYÞ. The Strong Law of Large Numbers ensures that the error in the approximation converges to zero as the simulation sample size t increases (see French and Rı´os Insua, 2000 for details). A number of strategies have been devised to obtain samples from the posterior, as e.g. Gibbs and Metropolis Hastings samplers. Interestingly, the whole process may be automated for relatively complex models, as in BUGS, a free software for Bayesian analysis using the Gibbs sampler (see Spiegelhalter et al., 1999). We illustrate the previous concepts with an example taken from Crowder (1978). It concerns the proportion of seeds that germinated on each of 21 plates, arranged according to a 2 2 factorial layout by seed and type of root extract. The data are shown in Table 3, where yi and ni are the number of germinations and the total number of seeds on the ith plate i ¼ 1; . . . ; n, respectively. The model is essentially a random effects logistic regression model, allowing for over-dispersion. If i is the probability of germination on the ith plate, we assume
logitði Þ ¼ 0 þ 1 x1i þ 2 x2i þ 12 x1i x2i þ "i
"i N 0; 2 where x1i , x2i are the seed type and root extract of the ith plate, and an interaction term 12 x1i x2i is included. 0 , 1 , 2 , 12 , 2 are given independent noninformative priors. Table 4 shows the different estimations for the logistic regression model obtained with maximum likelihood and the Gibbs sampler results using BUGS. Note, however, that in this case classical methods based on asymptotic properties are not appropriate because of the small sample sizes; thus confidence intervals are unrealistic. The results based on posterior inference show larger intervals, therefore acknowledging more uncertainty.
6. Spatial issues Spatial issues are of extreme importance in conservation biology as we try to relate biological phenomena to location-related variables (soil, climate, presence of human populations. . .) (e.g. Rubio and Escudero, 2000). They have also been of recent interest in Bayesian statistics (see Besag and Green, 1993, for general references and Gregoire et al., 1997, for references related to ecology). Here we outline a model which, we feel, may be relevant in conservation biology. The underlying theme is modelling with mixtures which, again, we feel is an important issue in this field for classification and discrimination purposes, as we may model population heterogeneity in a natural way. Diebolt and Robert (1994) provide a general framework for mixture modelling. For further details and extensions, see Marı´n et al. (2002a). Our goal is to model the spatial structure of a plant population. Our basic data are the locations of these plants. We shall consider an ecologically distributed
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
386
As another example, let us suppose we have 48 mixed measurements on the time of seed germination of a given species collected in a population (Table 5). We wish to know if the collected seeds come from plants occurring in two different microhabitats. Thus, we fit a mixture of two normal distributions with the same variance, so that each observation yi is assumed to be drawn from one (but we do not know which) of the two microhabitats (1 or 2). We denote Ti ¼ 1; 2 as the true microhabitat of the ith observation, (i ¼ 1; . . . ; 48), where microhabitat j (j ¼ 1; 2) has a normal distribution with mean j and variance 2 . We assume that an unknown fraction p of observations belongs to microhabitat 1 and, hence, the fraction ð1 pÞ belongs to microhabitat 2. The model can be written as
y i N T i ; 2
Table 4 Different estimations for the logistic regression model. Maximum likelihood and BUGS (Gibbs sampler results using BUGS software) Variable
Maximum likelihood SE
BUGS SE
constant (0 ) seed (1 ) extract (2 ) interaction (12 ) scale ( )
0.5480.167 0.0970.278 1.3370.237 0.8110.385 0.2360.110
0.5420.178 0.0280.340 1.3680.253 0.7920.426 0.2920.152
environment, for example regulated by wind. Under usual circumstances, normal winds will facilitate seed dispersal and establishment of individuals in nearby locations. Note that we shall have uncertainty about the number of clusters in the population. At the level of observed data, we induce clustering by assuming a mixture of normal models with an unknown number M of terms. Therefore, we assume that plant characteristics yi are distributed as
where Ti , i ¼ 1; . . . ; 48, follows a discrete distribution 1 with probability p Ti ¼ 2 with probability ð1 pÞ Note that this formulation may be easily generalised to additional components of the mixture. For the sake of computational efficiency we make a reparameterisation assuming that 2 ¼ 1 þ ; > 0: In this example 1 , , 2 and p, are given independent noninformative priors, including a uniform prior for p on ð0; 1Þ. Then, we apply a Markov chain Monte Carlo method, with the Gibbs sampler, using BUGS. We obtain samples (of size 10 000) from the posterior distributions of the parameters of interest. Then, we consider the mean and median of the samples as a summary of the distributions. A measure of precision of these values, the standard deviation and 0.95 posterior intervals are shown in Table 6. We obtain that 60% of the seeds would come from microhabitat 1. We can also see that the seeds from microhabitat 1 have a slightly shorter mean time of germination than the seeds from microhabitat 2 (Table 6). The difference in the mean time of germination might involve differences in seed physiology originated by contrasting microhabitats.
M X
qj N j j ; Sj ; i ¼ 1; . . . ; N: j¼1
By this, we understand that plants are clustered (normally) within one of M groups, each cluster centre being
j and each cluster dispersion matrix being j , j ¼ 1; . . . ; M. The proportion of plantsPwithin the jth group would be qj , j ¼ 1; . . . ; M, with M j¼1 qj ¼ 1. The model is completed with priors on the parameters, including the size M of the mixture, and the weights q ¼ ðq1 ; . . . ; qM Þ, as described in Marı´n et al. (2002a). We have already mentioned the computational difficulties associated with performing posterior inference. Section 3 illustrated a direct Monte Carlo approach, whereas in Section 4 we described how to undertake an MCMC based approach using BUGS. In other cases, however, such automated approaches are not feasible and we must perform specific analysis for the problem at hand.
Table 5 Simulated times of germination of seeds collected in a population with two potential microhabitats 529.0 535.3 537.5 542.8 548.9 552.8
530.0 535.4 538.3 543.0 549.0 552.9
532.0 535.9 538.5 543.5 549.4 553.2
533.1 536.1 538.6 543.8 549.9
533.4 536.3 539.4 543.9 550.6
533.6 536.4 539.6 545.3 551.2
533.7 536.6 540.4 546.2 551.4
534.1 537.0 540.8 548.8 551.5
534.8 537.4 542.0 548.7 551.6
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
387
Table 6 Sample estimators for the parameters obtained applying a Markov chain Montecarlo method with the Gibbs sampler Parameter
Mean
S.D.
2.5%
Median
97.5%
p
1
2
0.5982 536.7 548.9 3.768
0.08526 0.9265 1.256 0.6137
0.426 535.0 546.3 2.935
0.6016 536.7 548.9 3.668
0.7554 538.6 551.2 5.25
7. Concluding remarks The final emphasis in conservation biology should be decision-making issues concerning interventions. The intrinsically uncertain nature of conservation biology makes Bayesian analysis a natural approach within this field whenever statistical problems are considered. To some extent we have conveyed the relevance of Bayesian methods through various examples and models, covering areas of wide interest including space issues and population structure forecasting. We have outlined some typical problems and possible solutions under the Bayesian paradigm but now much further work has to be done. We hope that conservation scientists will emphasise Bayesian ideas in their future studies.
Acknowledgements Research supported by grants from CICYT, CAM and URJC. We are grateful to discussions with Adria´n Escudero, Jose´ M. Iriondo, the participants at the Threatened Plant Conservation Biology Seminar held at Universidad Polite´cnica de Madrid and the referees.
References Besag, J., Green, P., 1993. Spatial Statistics and Bayesian computation (with discussion). Journal Royal Statistics Society B 55, 25–37. Caswell, H., 2001. Matrix Population Models. Sinauer, Sunderland. Crowder, M.J., 1978. Beta–binomial Anova for proportions. Applied Statistics 27, 34–37.
Dennis, B., 1996. Should ecologists become Bayesian. Ecological Applications 6, 1095–1103. Diebolt, J., Robert, C., 1994. Estimation of finite mixture distributions through Bayesian sampling. Journal Royal Statistics Society B 56, 163–175. Ellison, A., 1996. An introduction to Bayesian inference for ecological research and environmental decision making. Ecological Applications 6, 1036–1046. French, S., Rı´os Insua, D., 2000. Statistical Decision Theory. Arnold, London. Gregoire, T.G., Billinger, D.R., Diggle, P.J., Russek-Cohen, E., Warren, W.G., Wolfinger, R.D. (Eds.), 1997. Modelling Longitudinal and Spatially Correlated Data: Methods, Applications, and Future Directions. Lecture Notes in Statistics 122. Springer-Verlag, New York. Marı´n, J.M., Muller, P., Rı´os Insua, D., 2002a. Spatial clustering with multi-level mixture models. Technical Report, U. Rey Juan Carlos. Marı´n, J.M., Pla, L.M., Rı´os Insua, D., 2002b. Inference for some stochastic processes related with sow farm management. Technical Report, U. Rey Juan Carlos. O’Hagan, A., 1994. Bayesian Inference. Arnold, London. Osenberg, C., Sarnelle, O., Goldberg, D., 1999. Meta-analysis in Ecology. Ecology 80, 1103–1104. Reckhow, K., 1990. Bayesian inference in non-replicated ecological studies. Ecology 71, 2053–2059. Rubio, A., Escudero, A., 2000. Small-scale spatial soil–plant relationship in semi-arid gypsum environment. Plant and Soil 220, 139–150. Schemske, D., Husband, B., Ruckelshaus, M., Goodwillie, C., Parker, J., Bishop, J., 1994. Evaluating approaches to the conservation of rare and endangered plants. Ecology 75, 2053–2059. Spiegelhalter, D.J., Thomas,A. and Best, N.G., 1999. WinBUGS Version 1.2 User Manual MRC Biostatistics Unit. http://www.mrcbsu.cam.ac.uk/bugs/. Wade, P., 2000. Bayesian methods in conservation biology. Conservation Biology 14, 1308–1316. West, M., Harrison, J., 1997. Bayesian Forecasting and Dynamic Models. Springer, New York.
Bayesian methods in plant conservation biology Juan M. Marı´n*, Raquel Montes Diez, David Rı´os Insua Statistics and Decision Sciences Group, Department of Experimental Sciences and Engineering, Universidad Rey Juan Carlos, Tulipa´n s/n. 28933 Mo´stoles, Madrid, Spain
Abstract We provide an introduction to Bayesian methods in conservation biology, illustrating inferences, prediction and decision making issues. After presenting the basic framework with a recovery plan evaluation problem, we illustrate more complex issues related to forecasting trends of structured populations using matrix population models and, finally, describe relevant topics in spatial and logistic regression problems. Computational and other implementation difficulties are also discussed. # 2003 Published by Elsevier Science Ltd. Keywords: Conservation biology; Bayesian analysis; Matrix population models; Recovery plans; Viability; Mixtures
1. Introduction The Bayesian framework provides a unified and coherent approach to solving inference and decision making problems under conditions of uncertainty. As such, we would expect Bayesian methods to have a great impact on conservation biology. Conservation scientists gather and analyse data to improve resource management in environmental problems. These data are usually affected by uncertainty and therefore require a statistical analysis. In environmental and ecological problems, the end use of such analysis is often to reach a decision, hence requiring uncertainty to be formally accounted for in the decision making process. Bayesian analysis seems, therefore, the appropriate approach to bioconservation statistical analysis, as it provides a coherent framework that integrates multiple sources of uncertain information (data, expert opinion) within decision making. French and Rı´os Insua (2000) provide a complete description. Many other advantages support the use of Bayesian methods in biological conservation studies. For example, prior information is frequently available, possibly from previous related studies from which our analysis might benefit, if appropriately modelled. Also, we sometimes need to combine several studies, see Osenberg et al. (1999) for a review on
* Corresponding author. Tel.: +34-916647475; fax: 916647434. E-mail address: [email protected] (J.M. Marı´n).
+34-
metanalysis issues in ecology, and the Bayesian approach provides a natural way of tackling those problems through hierarchical analysis. Spatial issues are also becoming increasingly important in ecology and recent Bayesian models provide useful tools to handle them. However powerful we might think the Bayesian approach is, we should admit that most statistical analyses in this field are still performed within the classical paradigm, although the situation seems to be gradually changing. As in other areas, controversy has permeated the eventual application of Bayesian methods in ecology—see Dennis (1996) for an example. In spite of its foundational strength, we believe that the best way to introduce and further advance these methods in conservation biology is through practice. This is precisely the goal of this paper, to provide a motivated introduction to Bayesian methods. Other reviews in the field, such as Reckhow (1990), Ellison (1996), or Wade (2000), have remained at a more conceptual level. We shall emphasise issues related to plant conservation biology. After introducing the basic Bayesian framework, we first illustrate the ideas with a simple example relating to the assessment of alternative actions in recovery plans. We then describe an alternative approach to matrix population models based on dynamic linear models. In Section 5, we illustrate a logistic model to forecast germination, comparing Bayesian and standard methods. Finally, we introduce a novel model to handle spatial issues, when spatial clustering of individuals takes place.
0006-3207/03/$ - see front matter # 2003 Published by Elsevier Science Ltd. doi:10.1016/S0006-3207(03)00124-1
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
380
2. An overview of Bayesian methods The Bayesian framework for inference and decision making problems can be easily described. Indeed, we feel that one of its strengths is, precisely, the ease with which basic concepts are put into place. Frequently one of the goals in conservation biology is to learn about one or more parameters (e.g. plant mortality rates, expected lifetimes or dispersion) which describe an ecological phenomenon of interest. Shall designate such parameters. To learn about them, we shall observe the phenomenon, collect data (e.g. counts of plants at various ages or their locations) and form the likelihood pðxjÞ which describes how data x depends (probabilistically) on the parameters , that is the distribution of data x when the parameter values are fixed at . The first novel ingredient in the Bayesian approach is to take into account another source of information about the parameters: very frequently, the conservation biologist will have access to expert information, based on past experience or related studies, which we shall model through the prior distribution pðÞ. We now need a way to combine both sources of information. This is done through Bayes’ formula. It provides the posterior distribution pðxjÞ, which is the distribution of the parameter given the data x: pðjxÞ ¼ Ð
pðÞpðxjÞ : pðÞpðxjÞd
Note that the posterior will be proportional ( /) to the numerator, which we shall represent pðjxÞ / pðÞpðxjÞ: The posterior distribution summarises all available information about the parameters and may be used to solve all inference problems of interest in conservation biology, such as point and interval estimation, hypothesis testing or prediction. For example, for point estimation, we could summarise the posterior through, for example, its mean, that is ð EðjxÞ ¼ pðjxÞd: To predict further y values of the phenomenon, we can use the predictive distribution. If we knew the value of , we would use the distribution pðyjÞ. Since there is uncertainty about , modelled through the posterior pðjxÞ, we integrate it out to get the predictive ð pðyjxÞ ¼ pðyjÞpðjxÞd: The ultimate aim of statistical research in conservation biology is usually to support decision making, for instance, to decide between several alternative actions in a recovery plan. For each action a and each future result
y, we would associate a consequence cða; yÞ. For example, if action a succeeds (y, the future result in this case), we would have incurred certain expenses, created some jobs and recovered a certain piece of land for conservation purposes. Such consequence will be evaluated with its utility uða; yÞ and we will choose the action maximising the expected utility ð uða; yÞpðyjxÞdy: Our aim now is to illustrate these ideas with several applications in plant conservation biology, showing their potentiality and their computational and modelling subtleties.
3. A simple example: evaluating recovery plans To illustrate the basic concepts, we provide examples concerning the evaluation of actions in recovery plans. We assume that a certain protocol establishes action A for a certain plant. Our aim is to evaluate its success. The parameter in this case is the probability 1 of a successful intervention through action A. Let us suppose that 9 out of 12 cases ended in success. If X1 is the number of successes in 12 trials, we have a binomial model Binð12; 1 Þ and, accordingly, we have that the probability of obtaining nine successful interventions in 12 trials given 1 is: 12 9 1 ð 1 1 Þ 3 : PrðX1 ¼ 9j1 Þ ¼ 9 As we are novel to the problem, let us view all values of 1 as being equally likely and use the uniform prior pð1 Þ ¼ 1; 1 2 ½0; 1 . Hence, applying Bayes’ formula, we get the posterior pð1 jX1 ¼ 9Þ / pð1 Þ PrðX1 ¼ 9j1 Þ / 1 19 ð1 1 Þ3 ; 1 2 ½0; 1 which corresponds to a beta distribution with parameters ¼ 10 and ¼ 4 that we represent as Beð10; 4Þ. Remember that a random variable has a Beð; Þ disÞ 1 ð1 Þ1 , 8 2 tribution, ; > 0, if pð Þ ¼ GGððþ ÞGðÞ ½0; 1 its mean being þ. Given the sequential nature of Bayes’ formula, this posterior would serve as prior for further studies. For example, suppose that in a subsequent study, three out of five implementations of our action A were successful. The new posterior would be pð1 jX1 ¼ 9; X2 ¼ 3Þ / pð1 jX1 ¼ 9ÞPrðX2 ¼ 3j1 Þ / 19 ð11 Þ3 13 ð11 Þ2 / 112 ð11 Þ5 ; 1 2 ½0; 1 corresponding to a Beð13; 6Þ distribution.
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
As mentioned, we may sometimes have expert information available about parameters. For example, let us suppose that we believe that higher values of 1 are more likely and we expect that around 80% of the times that this action is applied in a similar environment it will be successful. We could model this with a Beð4; 1Þ prior distribution. Note that, as required, its mean is 0:8 and its mode is 1. With such prior, the posteriors would be Beð13; 4Þ, after the first sample, and Beð16; 6Þ, after the second one. Fig. 1 provides the six distributions involved in this example. Those in the left column refer to the evolution of beliefs about 1 , the probability of successful implementation of action A stemming from the uniform prior, whereas those on the right stem from the Beð4; 1Þ prior. The left column shows a remarkable change in beliefs, from total uncertainty to a distribution fairly concentrated in the region [0.5,0.9]. Also note that the posterior becomes increasingly concentrated as more data become available. Similar comments hold for the right column, though we model a lot of prior information, and, therefore, the evolution of beliefs is not that remarkable. Finally, note that the bottom figures are fairly similar, showing the usual fact that, as data increase, prior relevance mitigates. The posterior distribution contains all relevant information for further inferences and predictions. The
381
following comments refer to the case that stems from uniform prior beliefs and we have observed nine successes and three failures, i.e. the current beliefs about 1 follow a Beð10; 4Þ distribution. Let us suppose we need a point estimate for the probability of success of action A. We could use the posterior mean, which is 10=ð10 þ 4Þ ¼ 0:72. Should we need an interval estimate, we could provide a region which covers 1 with posterior probability 0:9; for example, we have that pð1 2 ½0:505; 0:887 jX1 ¼ 9Þ 0:9; where we have chosen the interval so that pð1 4 0:505jX1 ¼ 9Þ 0:05 and pð1 5 0:887jX1 ¼ 9Þ 0:05. From an applied point of view, we may be more interested in computing predictive probabilities. For example, suppose that we are about to implement action A in seven more locations and we are interested in computing the probability of succeeding five or more times. We have, for each k; k ¼ 0; . . . ; 7, that ð PrðY ¼ kjX1 ¼ 9Þ ¼ PrðY ¼ kj1 Þpð1 jX1 ¼ 9Þd1 ð Gð14Þ 7 k 9 ð1 1 Þ3 d1 ¼ ð1 1 Þ7-k k 1 Gð10ÞGð4Þ 1 ð Gð14Þ 7 9þk ð1 1 Þ10k d1 ¼ k Gð10ÞGð4Þ 1 7 13!ð9 þ kÞ!ð10 kÞ! ¼ k 9!3!20!
Fig. 1. Distributions involved in the recovery plan evaluation problem. 1 : probability of a successful intervention through action A.
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
382
Therefore, PrðY 5 5jX1 ¼ 9Þ ¼
7 X PrðY ¼ kjX1 ¼ 9Þ; 0:6641: k¼5
The comparison of alternative actions follows a similar path. Suppose, for example, that another action B has been used in a similar context, resulting in 10 applications with Z ¼ 6 successes. We compare 1 , the probability of success with action A, with 2 , the probability of success with action B. The model we consider is Xj1 Binð12; 1 Þ Zj2 Binð10; 2 Þ 1 ; 2 Unif ½0; 1 ; independent: For comparison purposes, we compute r ¼ Prð1 5 2 jX1 ¼ 9; Z ¼ 6Þ. To do so, we see that, a posteriori, 1 jX1 ¼ 9 Beð10; 4Þ, 2 jZ ¼ 6 Beð7; 5Þ and that they are independent. The posterior distribution of 1 2 cannot be obtained in closed form, but we may easily simulate it. For example, we may draw 1000 values from the posterior distributions of 1 and 2 , obtain their differences 1 2 , count the number of such differences which are non-negative and divide by 1000 to approximate r. In our example, we undertook such a procedure to obtain r 0:772. This result suggests that action A seems better than action B, as we have a better chance of success with the former. However, to formally make a decision, we should consider the utility of various consequences. We shall use a simplified utility function here; for further details see French and Rı´os Insua (2000). Suppose that action A implies a considerable cost, but action B implies virtually no cost. Let us also suppose that the utility function is as described in Table 1 where, for example, if we implement action A successfully, the utility obtained is 0.8. Then, conditional on 1 and 2 , the expected utilities of the actions are, respectively, 0:81 þ 0ð1 1 Þ ¼ 0:81 and 2 þ 0:2ð1 2 Þ ¼ 0:2 þ 0:82 , whereas their posterior expected utilities are, respectively, 0:8Eð1 jX1 ¼ 9Þ ¼ 0:8
10 4 ¼ ; 14 7
and Table 1 Utilities for actions A and B in a recovery plan
Action A Action B
7 2 ¼ : 12 3 Therefore, we would actually recommend plan B, since its expected utility is higher.
0:2 þ 0:8Eð2 jZ ¼ 6Þ ¼ 0:2 þ 0:8
Succeeds
Does not succeed
0.8 1
0 0.2
4. A complex example: forecasting population viability One of the most important issues in conservation biology is the evaluation of a population’s stability in terms of its size at each stage of its biological cycle. The final goal is usually to predict whether a population is decreasing, increasing or is stabilised, hence assessing its viability. To do this, a demographic analysis and a methodology to study population growth trajectories and identify life history stages are required. Many authors, for example, Schemske et al. (1994) and Caswell (2001), have advocated the use of matrix population models in the context of conservation biology. Individuals are classified into several categories that reflect natural stages of the life cycle. Fig. 2 shows a hypothetical graph of a plant population model with five states: seed, seedling, juvenile, subadult and adult, with the corresponding valid transitions, pij probabilities and Fj per capita fertilities at each stage. For example, p12 describes the probability of a seed becoming a seedling after one year, and F5 is the number of seeds produced by an adult per year. As usual in conservation biology, we do not show a death stage either in the life cycle or in the transition matrix in Table 2 where there are implicit transitions to the death stage with the corresponding probability. For example, the probability of seed mortality would be 1 p11 p12 . A Markov chain model is implicitly assumed, transition probabilities are estimated and the evolution of the population size can be studied, for instance, to compare different management strategies. Thus, this procedure can be used to estimate seeding needs or the expected number of adults in a period of time. To characterise the process, we shall consider the probabilities given in the associated transition matrix in Table 2. The parameters defining our matrix population model process are, therefore, y=(p11,p12,p22,p33,p34,p35, p44,p45,p53,p54,p55,F5). We assume they are, a priori, independent. Inference procedures for such parameters follow. For the p parameters corresponding to each column, that is, the fates of a particular life stage after a time period, we shall assume a Dirichlet-multinomial model (see O’Hagan, 1994). This model, in which more than two classes are used, is a straightforward generalisation of the b-binomial model applied in Section 3. Let us suppose that there are three classes and let X ¼ ðX1 ; X2 ; X3 Þ designate the number of occurrences in classes 1, 2 and 3 after n trials, with X1 þ X2 þ X3 ¼ n. Then, we assume Xj Multðn; 1 ; 2 Þ, where i is the
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
383
Fig. 2. Transitions from different states of the life cycle of a plant. Table 2 Associated population transition matrix
Seed Seedling Juvenile Subadult Adult
Seed
Seedling
Juvenile
Subadult
Adult
p11 p12 0 0 0
0 p22 p23 0 0
0 0 p33 p34 p35
0 0 0 p44 p45
F5 0 p53 p54 p55
probability of an observation in class i, i ¼ 1; 2 and 3 ¼ 1 1 2 is the probability of an observation in class 3. A priori, Dirð1 ; 2 ; 3 Þ, a Dirichlet distribution with parameters 1 ; 2 and 3 . Remember that the Dirichlet distribution is a multivariate generalisation of the Beta distribution. For example, the Dirð1 ; 2 ; 3 Þ 1 1 2 1 3 1 density where 0 < i < 1 P is pð1 ; 2 ; 3 Þ / 1 2 3 and i i ¼ 1. A posteriori, we see that jð x1 ; x2 ; x3 Þ Dirð1 þ x1 ; 2 þ x2 ; 3 þ x3 Þ;
which 800 germinate, 100 die and 100 remain dormant. From the 800 seedlings, 600 pass to juveniles, 180 die and the rest remain in that stage, and so on until we complete a whole cycle. For the first column, we would have that, a posteriori, ð1 p11 p12 ; p11 ; p12 Þ Dirð100 þ 1; 100 þ 1; 800 þ 1Þ; with, respectively, posterior means 100þ1 100þ1 800þ1 ¼ 0:1007; ¼ 0:1007; ¼ 0:7986 1000þ3 1000þ3 1000þ3 and very small posterior variances, namely,
9:0197 105 ; 9:0197 105 ; 1:6019 104 : The whole process is completed resulting in very concentrated distributions around the posterior means, with matrix
with posterior means i þ xi ; i ¼ 1; 2; 3: Eði jxÞ ¼ P j þ n j
If there is little information from experts, we use noninformative priors. We start with 1000 seeds, out of
0
p11 B p12 B P¼B B0 @0 0
0 p22 p23 0 0
0 0 p33 p34 p35
0 0 0 p44 p45
1 0 0 C C p53 C C p54 A p55
384
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
The above transition matrix P is conditional on the parameter values, that is, the entries are actually pij j ¼ PrðXnþ1 ¼ jjXn ¼ i; Þ. From a Bayesian point of view, the transition matrix that we should consider is the corresponding marginal distribution. However, as our posteriors are very concentrated (see Marı´n et al., 2002b), we may approximate pij jdata by Pr Xnþ1 ¼ jjXn ¼ i; ^ ¼ pij j^ , where ^ are our posterior point estimates. F5 may be similarly estimated by counting the number of seeds produced by several adults. We could then use the methods in Caswell (2001) to undertake forecasting tasks. Other cases with more uncertainty in the point estimates require simulation, as we shall illustrate in a forthcoming paper. A shortcoming of such models is that the involved parameters are fixed through time. We describe here a Dynamic Linear Model that permits dynamic variation on such parameters. Our aim is to forecast the population structure several periods (years) ahead, by using the available population structure history. a Let us suppose that nst ; nslt ; njt ; nsa t and nt are, respectively, the number of seeds, seedlings, juveniles, subadults and adults at a given year t. We model the underlying population evolution in terms of an (autoregressive) multivariate dynamic linear model, with observation equations nstþ1 ¼ pt11 nst þ F5t nat þ "ts nsltþ1 ¼ pt12 nst þ pt22 nslt þ "tsl njtþ1 ¼ pt23 nst þ pt33 njt þ pt53 nat þ "tj t j t sa t a t nsa tþ1 ¼ p34 nt þ p44 nt þ p54 nt þ "sa t a t natþ1 ¼ pt35 njt þ pt45 nsa t þ p55 nt þ "a
where ptij ; i; j ¼ 1; . . . ; 5 are the probabilities of changing (surviving, growing,. . .) from stage i to stage j and F t;5 is the fertility per capita parameter, within period t. For example, the first equation describes that the number of seeds ns;tþ1 in year t þ 1 depends on the number of seeds pt11 nst , remaining at such stage from the previous year, and seeds F5t nat coming from adults in the previous year. Similarly, the last equation suggests that the number of adults natþ1 in year t þ 1 depends on the numbers of j adults nat , subadults nsa t and juveniles nt in the previous year. The error terms are assumed to be jointly independent and normal with zero mean and unknown variance. Note that if we drop these terms, we recover a standard matrix population model with time dependent parameters. We now assume that the parameters will vary randomly with time according to the system equations ptþ1 ¼ ptij þ tij ij F5tþ1 ¼ F5t þ t5
where the pij ’s satisfy the constraints ptij 5 0 and P5 t j¼1 pij ¼ 1, 8i, and the terms are also independent and normally distributed. This model, in fact, generalises classical matrix population models by considering random evolving parameters. As the parameters of the model may evolve with time, the model includes a system of equations to upgrade their values. Following the Bayesian paradigm, we would assign a prior distribution to parameters pij and F5 at time 0. The posterior distribution of parameters is updated at each stage, also giving the predictive distributions of the variables of interest several stages ahead. This class of models has been used in many other fields within engineering and applied sciences. Its main advantage is its flexibility to adapt to many situations, where the calculations are automatically performed by using the Kalman Filter procedure, a set of embedded linear equations that resume the main stages of estimations of the model. A full description of this technique is outside the scope of this paper, but may be seen, for example, in West and Harrison (1997).
5. Logistic regression Logistic regression is a popular technique in life sciences. It can be regarded as an extension of standard regression models. Its main aim is to predict the probability of success within a given trial in relation to a group of predictive variables. For example, in a simple logistic (binomial) regression, we may have independent observations from I populations. We assume that each one is distributed according to a binomial law: yi Binðni ; i Þ, for i ¼ 1; . . . ; I, where Eðyi Þ ¼ ni i , and that the explanatory variable x;i is available. Odds are defined as i =ð1 i Þ. The logistic regression specifies a linear structure for the logarithm of the odds (log odds): i logitði Þ log ¼ 0 þ 1 xi ; i ¼ 1; . . . ; I: 1-i Note that by varying xi , the term on the right of equation can take any real value, but i is always restricted between 0 and 1. Standard methods for analysing binomial regression data rely on asymptotic inference, requiring a great amount of data to provide accurate estimates. This problem is straightforwardly worked out by Bayesian methodology, allowing us to include prior information, use diagnostic studies, undertake model selection among different models and perform sensitivity analysis. This has been facilitated by the recent introduction of Markov chain Monte Carlo methods, which in our model permit computations for prediction and making inferences on regression coefficients.
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
To fix ideas, suppose we have data yi ; x0i ; i ¼ 1; . . . n, available where xi s are known vectors of k covariates and the yi s are independent binomial random variables with ni trials. The success probability for any single trial y with covariate x is Fðx0 Þ, i.e. 0 Fðx Þ Prðy ¼ 1jx; Þ, where is an unknown k-dimensional vector of regression coefficients, and F is Fðx0 Þ ¼
expðx0 Þ : 1 þ expðx0 Þ
The likelihood function for the complete data Y ¼ ðy1 ; . . . yn Þ is
385
Table 3 Proportion of germinated seeds of two varieties of Orobanche aegyptiaca in presence of bean and cucumber root extracts. (Crowder, 1978) O. aegyptiaca 75 Bean
O. aegyptiaca 73 Cucumber
Bean
Cucumber
y
n
y/n
y
n
y/n
y
n
y/n
y
n
y/n
10 23 23 26 17 –
39 62 81 51 39 –
0.26 0.37 0.28 0.51 0.44 –
5 53 55 32 46 10
6 74 72 51 79 13
0.83 0.72 0.76 0.63 0.58 0.77
8 10 8 23 0 –
16 30 28 45 4 –
0.50 0.33 0.29 0.51 0.0 –
3 22 15 32 3 –
12 41 30 51 7 –
0.25 0.54 0.50 0.63 0.43 –
n Y n y ni 0 yi 1 F x0i i i F xi LðjYÞ yi
y: number of germinated seeds in each plate; n: total number of seeds in each plate.
For a given prior distribution on , namely ðÞ, to obtain posterior and predictive distributions we need to compute the posterior distribution of :
yi Binðni ; i Þ
i¼1
LðjYÞ ðÞ : LðjYÞ ðÞd As many important quantities are obtained through cumbersome and intractable integrals involving ðjYÞ, approximations must be undertaken. Monte Carlo methods provide discrete approximations to the posterior distribution. Given a function hðÞ, the posterior expectation E½hðÞjY may be approximated by
ðjYÞ ¼ Ð
ð t
1X E½hðÞjY ¼ hðÞ ðjYÞd h i ; t i¼1 t for a sample i i¼1 of size t, from the posterior distribution ðjYÞ. The Strong Law of Large Numbers ensures that the error in the approximation converges to zero as the simulation sample size t increases (see French and Rı´os Insua, 2000 for details). A number of strategies have been devised to obtain samples from the posterior, as e.g. Gibbs and Metropolis Hastings samplers. Interestingly, the whole process may be automated for relatively complex models, as in BUGS, a free software for Bayesian analysis using the Gibbs sampler (see Spiegelhalter et al., 1999). We illustrate the previous concepts with an example taken from Crowder (1978). It concerns the proportion of seeds that germinated on each of 21 plates, arranged according to a 2 2 factorial layout by seed and type of root extract. The data are shown in Table 3, where yi and ni are the number of germinations and the total number of seeds on the ith plate i ¼ 1; . . . ; n, respectively. The model is essentially a random effects logistic regression model, allowing for over-dispersion. If i is the probability of germination on the ith plate, we assume
logitði Þ ¼ 0 þ 1 x1i þ 2 x2i þ 12 x1i x2i þ "i
"i N 0; 2 where x1i , x2i are the seed type and root extract of the ith plate, and an interaction term 12 x1i x2i is included. 0 , 1 , 2 , 12 , 2 are given independent noninformative priors. Table 4 shows the different estimations for the logistic regression model obtained with maximum likelihood and the Gibbs sampler results using BUGS. Note, however, that in this case classical methods based on asymptotic properties are not appropriate because of the small sample sizes; thus confidence intervals are unrealistic. The results based on posterior inference show larger intervals, therefore acknowledging more uncertainty.
6. Spatial issues Spatial issues are of extreme importance in conservation biology as we try to relate biological phenomena to location-related variables (soil, climate, presence of human populations. . .) (e.g. Rubio and Escudero, 2000). They have also been of recent interest in Bayesian statistics (see Besag and Green, 1993, for general references and Gregoire et al., 1997, for references related to ecology). Here we outline a model which, we feel, may be relevant in conservation biology. The underlying theme is modelling with mixtures which, again, we feel is an important issue in this field for classification and discrimination purposes, as we may model population heterogeneity in a natural way. Diebolt and Robert (1994) provide a general framework for mixture modelling. For further details and extensions, see Marı´n et al. (2002a). Our goal is to model the spatial structure of a plant population. Our basic data are the locations of these plants. We shall consider an ecologically distributed
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
386
As another example, let us suppose we have 48 mixed measurements on the time of seed germination of a given species collected in a population (Table 5). We wish to know if the collected seeds come from plants occurring in two different microhabitats. Thus, we fit a mixture of two normal distributions with the same variance, so that each observation yi is assumed to be drawn from one (but we do not know which) of the two microhabitats (1 or 2). We denote Ti ¼ 1; 2 as the true microhabitat of the ith observation, (i ¼ 1; . . . ; 48), where microhabitat j (j ¼ 1; 2) has a normal distribution with mean j and variance 2 . We assume that an unknown fraction p of observations belongs to microhabitat 1 and, hence, the fraction ð1 pÞ belongs to microhabitat 2. The model can be written as
y i N T i ; 2
Table 4 Different estimations for the logistic regression model. Maximum likelihood and BUGS (Gibbs sampler results using BUGS software) Variable
Maximum likelihood SE
BUGS SE
constant (0 ) seed (1 ) extract (2 ) interaction (12 ) scale ( )
0.5480.167 0.0970.278 1.3370.237 0.8110.385 0.2360.110
0.5420.178 0.0280.340 1.3680.253 0.7920.426 0.2920.152
environment, for example regulated by wind. Under usual circumstances, normal winds will facilitate seed dispersal and establishment of individuals in nearby locations. Note that we shall have uncertainty about the number of clusters in the population. At the level of observed data, we induce clustering by assuming a mixture of normal models with an unknown number M of terms. Therefore, we assume that plant characteristics yi are distributed as
where Ti , i ¼ 1; . . . ; 48, follows a discrete distribution 1 with probability p Ti ¼ 2 with probability ð1 pÞ Note that this formulation may be easily generalised to additional components of the mixture. For the sake of computational efficiency we make a reparameterisation assuming that 2 ¼ 1 þ ; > 0: In this example 1 , , 2 and p, are given independent noninformative priors, including a uniform prior for p on ð0; 1Þ. Then, we apply a Markov chain Monte Carlo method, with the Gibbs sampler, using BUGS. We obtain samples (of size 10 000) from the posterior distributions of the parameters of interest. Then, we consider the mean and median of the samples as a summary of the distributions. A measure of precision of these values, the standard deviation and 0.95 posterior intervals are shown in Table 6. We obtain that 60% of the seeds would come from microhabitat 1. We can also see that the seeds from microhabitat 1 have a slightly shorter mean time of germination than the seeds from microhabitat 2 (Table 6). The difference in the mean time of germination might involve differences in seed physiology originated by contrasting microhabitats.
M X
qj N j j ; Sj ; i ¼ 1; . . . ; N: j¼1
By this, we understand that plants are clustered (normally) within one of M groups, each cluster centre being
j and each cluster dispersion matrix being j , j ¼ 1; . . . ; M. The proportion of plantsPwithin the jth group would be qj , j ¼ 1; . . . ; M, with M j¼1 qj ¼ 1. The model is completed with priors on the parameters, including the size M of the mixture, and the weights q ¼ ðq1 ; . . . ; qM Þ, as described in Marı´n et al. (2002a). We have already mentioned the computational difficulties associated with performing posterior inference. Section 3 illustrated a direct Monte Carlo approach, whereas in Section 4 we described how to undertake an MCMC based approach using BUGS. In other cases, however, such automated approaches are not feasible and we must perform specific analysis for the problem at hand.
Table 5 Simulated times of germination of seeds collected in a population with two potential microhabitats 529.0 535.3 537.5 542.8 548.9 552.8
530.0 535.4 538.3 543.0 549.0 552.9
532.0 535.9 538.5 543.5 549.4 553.2
533.1 536.1 538.6 543.8 549.9
533.4 536.3 539.4 543.9 550.6
533.6 536.4 539.6 545.3 551.2
533.7 536.6 540.4 546.2 551.4
534.1 537.0 540.8 548.8 551.5
534.8 537.4 542.0 548.7 551.6
J.M. Marı´n et al. / Biological Conservation 113 (2003) 379–387
387
Table 6 Sample estimators for the parameters obtained applying a Markov chain Montecarlo method with the Gibbs sampler Parameter
Mean
S.D.
2.5%
Median
97.5%
p
1
2
0.5982 536.7 548.9 3.768
0.08526 0.9265 1.256 0.6137
0.426 535.0 546.3 2.935
0.6016 536.7 548.9 3.668
0.7554 538.6 551.2 5.25
7. Concluding remarks The final emphasis in conservation biology should be decision-making issues concerning interventions. The intrinsically uncertain nature of conservation biology makes Bayesian analysis a natural approach within this field whenever statistical problems are considered. To some extent we have conveyed the relevance of Bayesian methods through various examples and models, covering areas of wide interest including space issues and population structure forecasting. We have outlined some typical problems and possible solutions under the Bayesian paradigm but now much further work has to be done. We hope that conservation scientists will emphasise Bayesian ideas in their future studies.
Acknowledgements Research supported by grants from CICYT, CAM and URJC. We are grateful to discussions with Adria´n Escudero, Jose´ M. Iriondo, the participants at the Threatened Plant Conservation Biology Seminar held at Universidad Polite´cnica de Madrid and the referees.
References Besag, J., Green, P., 1993. Spatial Statistics and Bayesian computation (with discussion). Journal Royal Statistics Society B 55, 25–37. Caswell, H., 2001. Matrix Population Models. Sinauer, Sunderland. Crowder, M.J., 1978. Beta–binomial Anova for proportions. Applied Statistics 27, 34–37.
Dennis, B., 1996. Should ecologists become Bayesian. Ecological Applications 6, 1095–1103. Diebolt, J., Robert, C., 1994. Estimation of finite mixture distributions through Bayesian sampling. Journal Royal Statistics Society B 56, 163–175. Ellison, A., 1996. An introduction to Bayesian inference for ecological research and environmental decision making. Ecological Applications 6, 1036–1046. French, S., Rı´os Insua, D., 2000. Statistical Decision Theory. Arnold, London. Gregoire, T.G., Billinger, D.R., Diggle, P.J., Russek-Cohen, E., Warren, W.G., Wolfinger, R.D. (Eds.), 1997. Modelling Longitudinal and Spatially Correlated Data: Methods, Applications, and Future Directions. Lecture Notes in Statistics 122. Springer-Verlag, New York. Marı´n, J.M., Muller, P., Rı´os Insua, D., 2002a. Spatial clustering with multi-level mixture models. Technical Report, U. Rey Juan Carlos. Marı´n, J.M., Pla, L.M., Rı´os Insua, D., 2002b. Inference for some stochastic processes related with sow farm management. Technical Report, U. Rey Juan Carlos. O’Hagan, A., 1994. Bayesian Inference. Arnold, London. Osenberg, C., Sarnelle, O., Goldberg, D., 1999. Meta-analysis in Ecology. Ecology 80, 1103–1104. Reckhow, K., 1990. Bayesian inference in non-replicated ecological studies. Ecology 71, 2053–2059. Rubio, A., Escudero, A., 2000. Small-scale spatial soil–plant relationship in semi-arid gypsum environment. Plant and Soil 220, 139–150. Schemske, D., Husband, B., Ruckelshaus, M., Goodwillie, C., Parker, J., Bishop, J., 1994. Evaluating approaches to the conservation of rare and endangered plants. Ecology 75, 2053–2059. Spiegelhalter, D.J., Thomas,A. and Best, N.G., 1999. WinBUGS Version 1.2 User Manual MRC Biostatistics Unit. http://www.mrcbsu.cam.ac.uk/bugs/. Wade, P., 2000. Bayesian methods in conservation biology. Conservation Biology 14, 1308–1316. West, M., Harrison, J., 1997. Bayesian Forecasting and Dynamic Models. Springer, New York.