Modeling Tips
Avoiding Mistakes in Population Modeling
H. Resit Akçakaya
The following list includes some of the mistakes modelers make in developing ecological models at the population or metapopulation levels, for example for population viability analysis. It grew out of my experiences in reviewing and editing papers on modeling, answering questions from RAMAS users, as well as my own mistakes.
My goal is mostly to help users of RAMAS Metapop and RAMAS GIS avoid these common mistakes. Although the list focuses on models developed using these programs, most (if not all) of these issues are relevant for models developed in any program or using any modeling platform.
Note that whether a model feature is a mistake or not depends on the context (the species, other model components and features, etc.). Thus, for example, modeling only females is a mistake only under some circumstances. Such circumstances are discussed in the links below.
I welcome any feedback, comments, and suggestions, including mundane things like misspellings and broken links, and (more importantly) suggestions about additional types of mistakes, references, approaches, etc.
General Mistakes
Invalid model assumptions
Model too complex
Model too simple
Internal inconsistency (Not paying attention to error messages)
Modeling Demographic Structure
No demographic structure (scalar models for age-structured populations)
Bias in fecundity estimation
Fecundity vs. maternity
Not incorporating proportion breeding
Not incorporating sex ratios
Extra stage
Fecundity of first age class or stage
Bias in survival rate estimation
Uncertainty in survival rate estimation
Using survivorship instead of survival rate
Survival rates in the diagonal
Ignoring constraints
Too many (or too few) age classes or stages
Modeling Density Dependence
Using (the wrong type of) density dependence
Not using density dependence
Bias in Rmax estimation
Underestimating maximum growth rate (Rmax)
Overestimating maximum growth rate (Rmax)
Incorrectly modeling impact under density dependence
Toxicity vs. harvest under density dependence
Adding Variability (Stochasticity)
Ignoring variability
Not using demographic stochasticity
Environmental stochasticity distribution
Ignoring correlations
Too few replications
Random fluctuations vs. regular oscillations
Using catastrophes instead of environmental stochasticity
Overestimating catastrophe impact
Not incorporating delayed effects of catastrophes
Overestimating variation
Using the wrong standard deviation
Not correcting truncations
Phantom individuals
Duration (simulation time horizon) too long or too short
Confusing uncertainty and variability
Modeling Spatial Structure
Ignoring spatial structure
Too many (or too few) populations
Ignoring spatial correlation
Wrong dispersal rates
Symmetric dispersal rates
Too high (or too low) map resolution
Presenting Results
Implicit assumptions
Ignoring uncertainty
Emphasizing deterministic (not probabilistic) results
Using absolute (not relative) predictions
Estimating risk of extinction rather than decline
Impact Assessment
Failing to incorporate cumulative impacts
Selecting the wrong spatial scale (or a single scale)
Reference population includes impact effects
Uncertainty masking impact
Underestimating impacts of toxicity, habitat loss, habitat fragmentation, and harvest
Difference between experimental and model time steps
Wrong time horizon or threshold
General Mistakes
Invalid model assumptions
Every model makes assumptions, and many assumptions are not or cannot be validated (if they were validated, they would be parameters, not assumptions). However, there are several types of checks that a user can make to ensure that the model assumptions are consistent with each other, and with the available qualitative information about the species' life history.
It is important to remember that a program such as RAMAS GIS is not a model, but a tool that helps users build their own models. RAMAS GIS has been validated at the program level by checking all subunits and algorithms of the programs, by checking the lower-level algorithms for consistency (e.g., random distributions of vital rates), and by making sure that the program does what is described in the manual and help files. The validation at the model level is the user's responsibility. For tips on how to do this, see the following topics in the RAMAS Metapop help files (the sections given below provide short introductions):
Validation and Testing
RAMAS GIS and RAMAS Metapop have been validated at the program level by checking all subunits and algorithms of the subprograms, by making sure that the program does what is described in the manual, by checking the lower-level algorithms for consistency (e.g., random distributions of vital rates), etc. RAMAS GIS and RAMAS Metapop are not models, but tools which help users build models. The validation at the model level is your responsibility. To validate a model, start with a review of the types of results that the program gives, and make sure that they will address the question you are asking. Next, you should review the assumptions and limitations of the program (see the manual) and convince yourself that these fit the species you are modeling. If they don't fit (or you are not sure whether they do or not) then at least you are aware of these assumptions and limitations. The major component of the validation involves your input of the model parameters into the program. It is essential that you understand how these parameters are used by the program, and not rely on your intuitive understanding based on the names given to various parameters in the program. The same names or symbols may be used in entirely different contexts by different programs or textbooks, leading to several common mistakes. If you have large amounts of data, you can also validate your model by comparing the model predictions with observations (e.g., see Brook et al. 2000). However desirable this type of validation may be, it is difficult and often impossible for stochastic models since you need observations on abundances of several populations as well as their variation and risks of local extinction.
The manual provides several suggestions for testing or checking your model.
Program Limitations and Assumptions
Models in RAMAS are expressed in discrete time. The program does not use continuous-time models. We believe that most modeling needs in population dynamics can be met with (and perhaps are most naturally expressed as) discrete-time formulations.
Since RAMAS GIS and RAMAS Metapop models can only address a single species, they cannot explicitly represent competition, predation, mutualism or other interspecific interactions. These interactions can be modeled as constant, stationary (randomly fluctuating with constant mean), or deterministically varying (e.g., cyclic) influences on demographic parameters. In addition, it may be possible to model such interactions in a user-written code linked to the program.
Population dynamics within each population is modeled with a stage- or age-structured matrix model, which may also include sex structure. The model may contain non-linearities describing density dependence relations between matrix elements and the population abundance. The program also allows used-written code to be linked to a model. This provides a lot of flexibility that will let you model most species. However, some very complex structures (for example those requiring complicated rules representing social structure) may be difficult or impossible to incorporate unless they can be simplified to a matrix model.
Within-population correlations among vital rates are restricted to 1 (full), 0 (none) and -1. (Among populations, any correlation between -1 and +1 can be specified.)
A catastrophe can either be regional (occur at the same time for all the populations) or local (occur independently for each population).
If a model includes two different types of catastrophes, the correlation between the two catastrophes (in time) can be zero, maximum positive or maximum negative.
The program does not estimate any input parameters. The manual provides background information, simple examples, and references to the relevant literature to help the users with parameter estimation. The program includes sample files that are parameterized from literature data on several species.
The program makes certain assumptions about the way parameters are combined into the overall model. These are implicit in the equations used in the program, and include dispersal-distance and correlation-distance functions, density-dependent dispersal, density dependence functions, etc. For the exact methods and equations the program uses, see discussions in the manual (particularly the algorithm in Appendix I).
In addition to these limitations, any particular model built in RAMAS will have a set of assumptions, characterized by its parameters (or lack of them). As with any other program or model, how realistic the metapopulation model is largely depends on the amount and type of data. (See "Avoiding mistakes in population modeling"). This program will, hopefully, let the user utilize such data effectively in predicting extinction risks, and perhaps even guide additional field work by helping the user identify parameters to which the risk results seem to be most sensitive. However, no program or model is a substitute for empirical work.
Summary of Transparency Features
There are several features in RAMAS Metapop and RAMAS GIS that are specifically designed to make the programs transparent.
-
The methods used in the program are described in detail in the manual as well as in several peer-reviewed publications.
-
The source code of the program is summarized in the form of pseudo-code in a file distributed with the program, as well as in an appendix to the manual.
-
You can link your own code to the program, with a "user-defined function" option.
-
Several examples of "user-defined functions" are provided with the program, including their full source codes as well as the compiled versions.
-
The manual and the help file list the basic assumptions and technical limitations of the program.
-
The program allows you to print a list of "Model Summary and Assumptions" for your model.
-
The program displays your input parameters in the form of various graphs that allow you to examine the properties of your model.
-
The program indicates which input parameters it is not using (e.g., if catastrophe probability is zero, then the magnitude of the catastrophe impact is not relevant).
-
The manual and the help file include the exact definitions of the model parameters, how they are used by the program, tips on parameter estimation, and suggestions for various methods of testing and validating your model.
-
The program makes many checks, and displays warnings if it detects potential problems in your model.
-
The help file provides a list and detailed explanation of these warning messages, as well as suggestions about how to correct your model.
These help topics make a number of suggestions. Two specific examples are given below.
1. Print a summary of your model (select "Export Model Summary" from the File menu to save a file that you can view and print using your browser), and examine the features and parameters of the model.
2. Examine the structure of the program, which is summarized in the form of pseudo-code in a file distributed with the program (RAMAS Metapop Algorithm.pdf in the program folder). This describes the way parameters are combined into the overall model. It is important that you understand how the program uses each parameter, and not rely on your intuitive understanding based on the parameter's name. The same name or symbol may be used in entirely different contexts by different programs or textbooks.
Also, read about how to check for and correct internal inconsistencies of your model.
Model too complex
Many modelers are tempted to add a lot of complexity to realistically reflect the life history of the species. However, complexity means more parameters to estimate, for which data may not be available. Attempts to include more details than can be justified by the quality of the available data may result in decreased predictive power (Ginzburg and Jensen 2004). In many cases, simple assumptions may be preferable to adding new factors or parameters.
For example, if the number of individuals in older age classes is small, then pooling these data into a "composite age class" (e.g., a 3+ "age" class that includes all individuals that are 3 years old and older) would allow a more precise estimate of the survival rate of these individuals. Of course, this would mean making several assumptions (survival and fecundity do not change after age 3; there is no senescence), which should be compared to what is known about the biology of the species.
Another example is adding various factors (e.g., inbreeding depression, behavioral traits, Allee effects, etc.) that affect the dynamics of very small populations. Although it is preferable to model each of these factors when data are available, often adding these factors involves a lot of guesswork and assumptions. When data are not available, it may be better to simply ignore these factors, but focus on risk of decline to a threshold greater than zero (instead of risk of extinction). When the model includes multiple populations, this can be done in two different ways: (1) use an extinction threshold for the whole population (in Stochasticity dialog), (2) use a local threshold for each population. See the help files and the manual for setting these thresholds, and the difference between the two.
Another problem with complex models is that it is more difficult to interpret or understand their results, some of which may be counter-intuitive. In many cases, even when sufficient data exist to build a complex model, it is advisable to start from a very simple version of your model (e.g., scalar, deterministic, density-independent), and add one feature at a time, running the model and examining its results at each step. Alternatively (if you have already build a complex model), you can remove one feature of the model at a time. A process of addition or elimination may let you understand the factors that contribute to the results, and notice problems in your model.
Each feature is removed in a different way. For example, to remove all variability (all stochastic elements) and run a deterministic simulation, set Replications to 0 (in General information dialog). To only eliminate catastrophes, go to Catastrophes dialog and uncheck each box under "Affects". To remove random fluctuations in vital rates, go to Standard deviations matrix dialog and set all elements to zero by using "Autofill" with CV=0. To remove age/stage structure, go to Stages dialog and delete all but one of the stages, then go to Stage Matrix dialog and enter the population growth rate ("finite rate of increase") as the only "stage matrix" element, then go to Initial abundances dialog to enter the correct value. To remove density dependence, go to Density dependence dialog, select "All populations have the same density dependence type", and set this type to "Exponential". To remove all spatial features, go to Populations dialog and delete all but one of the populations.
Model too simple
Some modelers make the opposite mistake of building a model that is (1) too simple to use all the available data and incorporate all information about the life history of the species, or (2) too simple to properly address the question at hand, or (3) too simple to incorporate all the known threats to the species. The risk in building a very simple model is that it may lead to an underestimate of the risk of decline or extinction (because all the threats may not be modeled) or to an inefficient use of the available information. The question of the appropriate level of complexity (i.e., the trade-off between realism and functionality) depends on:
1. characteristics of the species under study (e.g., its ecology),
2. what you know of the species (the availability of data), and
3. what you want to know or predict (the questions addressed).
For example, if the question to be addressed relates to harvest (hunting, fishing, etc.) of the population, or population management by translocating individuals, the results will probably depend sensitively on which age classes are harvested or translocated. In such cases, a simple model that does not include the population's age structure may not be very useful. Similarly, if different age classes show different physiological sensitivities to a particular toxin, then an ecotoxicological assessment using a structured model would be more appropriate than one that uses a scalar or unstructured model (one that does not have age or stage structure).
Some modelers are under the impression that data from regular census of a population (i.e., "count data") can only be used to develop a scalar model. This impression perhaps arises because a population viability analysis using a scalar model is sometimes (confusingly) called a "count-based PVA." However, if the census is not restricted to only the total population, but also includes a count of individuals in each age class or stage, then such data can be used to develop an age or stage structured model.
An objective way to determine the appropriate level of complexity is to fit statistical models of different complexity to the available data, and compare their goodness of fit using the Akaike Information Criterion (AIC). If you are not sure of the appropriate level of complexity, you may want to use multiple models at different levels of complexity, instead of a single model, and consider the range of results you obtain as part of model uncertainty.
Internal inconsistency (Not paying attention to error messages)
A common mistake is to ignore signs that some aspects of the model may be inconsistent with each other. In the case of RAMAS GIS, these signs include a large number and variety of error and warning messages. Some of these messages are encountered while entering model parameters; others appear while running a simulation. An error message indicates that there is a potential mistake in your model. The warning given by these messages is one of the advantages of using a professional program like RAMAS. So, if you see such a message, read about it in the help files to see what it means, to determine if it needs a corrective action, and to learn what options are suggested for correction.
Another sign of inconsistency is that a parameter you had entered earlier becomes unavailable for editing (the parameter box will become gray, and you will not be able to change the parameter). This indicates that the program is ignoring the parameter, usually because of the value of some other parameter (e.g., if local catastrophe probability is zero, local multiplier that determines catastrophe impact is irrelevant and is ignored). Check each set of input parameters (all items under the Model menu) to see which parameters are not available for editing.
Another way to check for inconsistencies is to examine the various graphs and tables the program uses to summarize the input data. For this, go to the Populations dialog (under the Model menu) and click on the "Display" button, which lets you see a variety of input summaries for each population.
If you are programming your own model (or using the "user-defined density dependence" feature in RAMAS Metapop), internal model inconsistencies may be the most difficult mistakes to notice and correct. Two strategies may help: First, start from a very simple version of your model, and gradually add more complex features, running the model at each step and keeping an eye on any counter-intuitive results. Second, compare your results to those from a model which is developed by someone else but has similar features.
Modeling Demographic Structure
No demographic structure (scalar models for age-structured populations)
Scalar models (with no age or stage structure within populations) are often used in cases where the only available data are a time series of population size estimates (these methods are also known as "count-based models" and "diffusion approximation"). However, if the population being modeled has age structure, a scalar model of this population may overestimate the variability in the population size, and hence overestimate the risks faced by the population (Holmes 2004).
An analysis of scalar and structured models for a set of populations has indicated a precautionary bias (overestimated risks of decline) by scalar models, and a set of simulations has indicated that the bias increases as a function of the generation time of the species (Dunham et al. 2006). Correcting the bias (e.g., based on the species' generation time) seems difficult if not impossible, because the bias is not a simple function of generation time, and because any deviation of the initial age structure from the stable age structure adds uncertainty (Dunham et al. 2006).
Future developments in the analysis of time series data for building stochastic scalar models may address this issue. Until then, if a population is known to have age structure, but age-structured data are not available to build a matrix model, the results of scalar models of this population should be viewed as preliminary (and likely precautionary).
Bias in fecundity estimation
Fecundity should be based on total productivity of individuals over the time step of the model. Fecundity calculations based on a single nesting attempt underestimate fecundity. See Anders & Marshall (2005) for methods of estimating season-long productivity in landbirds.
Note that in addition to productivity or fertility (number of fledglings per pair per time step), the fecundity values in a matrix model must also incorporate survival rate, sex ratio, and proportion breeding.
A common mistake is to estimate maternity (see below) as juvenile:adult ratio, in cases where the probability of recording (seeing, counting, capturing) juveniles and adults are different. In such cases, the juvenile:adult ratio must be multiplied by the ratio of the relative capture (or sighting) probabilities of adults to juveniles. For an example, see Ryu et al (2016).
Fecundity vs. maternity
A common mistake is to use a measure of maternity (m(x), e.g., number of eggs per female; number of fledglings per nest, etc.) in the stage matrix, instead of fecundities, F(x). In a matrix model, the fecundities must incorporate two types of survival:
a. survival of breeders (from census to next breeding)
b. survival of newborns (from birth to next census)
Most matrix models assume either post-breeding or pre-breeding census, and thus either (a) or (b) is assumed to be 1.0 (but never both). In using maternities or life table data to estimate the stage matrix, the exact definition of m(x) is very important; check the source of the information in order to use the correct index. For more information, see Akçakaya (2000a), Caswell (2001) and Akçakaya et al. (1999).
Not incorporating proportion breeding
In a matrix model, the fecundities must incorporate the proportion of individuals breeding in an age class or stage. Fecundities have units of "per individual" or "per female" in that age class or stage (not "per breeder"). If, for example, two-year old breeders and three-year old breeders have the same fecundity, but only half of two-year old individuals and all of three-year old individuals breed, then fecundity of the two-year-old age class will be half of the fecundity of the three-year-old age class.
The fecundities must include the proportion of breeding individuals, regardless of the type of sex structure or the mating system. In RAMAS Metapop, when females and males are modeled in separate stages (i.e., when the model includes sex structure), the proportion of breeders is entered in the Stages dialog. However, the parameter in the Stages dialog is used only to determine the limiting sex. If less than 100% of individuals breed, then this must be reflected in the stage matrix, as described above (even if the model includes sex structure, and the proportion of breeders have been entered in the Stages dialog).
Not incorporating sex ratio
Sex ratio is the proportion of females in a population or in an age class. Sex ratio at birth is the proportion of daughters among the offspring. In a matrix model, the sex ratio must be incorporated into the estimation of fecundities.
If you are modeling only the female population, the fecundities should be in terms of "number of daughters per female", and the initial abundances should be the number of females in each age class.
If you are modeling all individuals, the fecundities should be in terms of "offspring per individual", and the initial abundances should be the number of all individuals in each age class.
If you are modeling females and males separately, fecundities are entered in two different elements of the stage matrix, as "number of daughters per female" and "number of sons per female" (assuming a monogamous or a polygynous mating system). For more information, see "Modeling sex structure" in the RAMAS Metapop help file.
Extra stage
Some matrix models include a stage that should not be included in the matrix. One example of this mistake is including a "seed" stage in a plant model, even though seeds do not need to "wait" for a whole time step (e.g., for 1 year) before they germinate. An example of a similar mistake in a model of an animal population might be an "egg" stage.
For an example of this common mistake, consider the following model. Flowering plants produce 50 seeds/plant per year, of which 20% germinate to become "seedling" in the next time step. There is no seed bank, so any seed that does not germinate by the next census dies. We also assume that the census is pre-reproductive. In this case, the first row of the correct matrix would have the following form (ignore the other rows).
In this case, seeds are not included in the model, because in a pre-reproductive census, seeds would not be counted. If this matrix were parameterized according to post-productive census, it would be correct to name the first stage"seed", but in that case, the fecundity value must include the survival rate to the flowering stage.
However, some models incorrectly include a "seed" stage:
This model is wrong (as a model of the same species as described above), because it introduces an extra time step between seed production and germination. It would also be wrong under post-productive census. For an example of this mistake, see Werner and Caswell (1977), and the discussion of the corrected model in Caswell (2001, p. 60-62, figs 4.3, 4.4).
As mentioned above, a similar mistake that is common in models of animal populations is including stages such as "egg" and "hatchling". This is a mistake if an egg does not take one time step to hatch. In general, the life cycle of the organism, in conjuction with the time step (projection interval) of the model determine the stages. Just because a phase in the life cycle of a species is important does not mean that that phase will always be represented with a stage in the matrix. In most cases, if the projection interval is longer than the time individuals spend in that phase of their life cycles, that phase will not be represented as a separate stage in the matrix model.
Case with seed bank: Now consider a similar model of a species with seed bank.Flowering plants produce 50 seeds/plant per year, of which 20% germinate to become "seedling" in the next time step. 30% of the seeds remain in the seed bank, of which only 5% germinate in the following year, and the rest die (i.e., no seed remains in the seed bank longer than a year). In this case, the first 2 rows of the correct matrix would have the following form:
Fecundity of first age class or stage
A common mistake is to set the fecundity of the first age class (the upper left corner element of the stage matrix) to zero, when it should be greater than zero. If a model has, for example, annual age classes, and some of the individuals that are born/fledged in one year can breed in the following year (when they are almost, but not quite, 1-year old), then this element should not be zero.
If at least some individuals can reproduce when they are, say, 11-months old (or earlier), then the upper left corner element of the stage matrix should be greater than zero (in a model with annual age classes). Similarly, if age of first reproduction is 2 years, then there should be only one zero in the first row, not two zeros.
Bias in survival rate estimation
Immigration and emigration, if not properly accounted for, can bias survival rate estimates. In order to observe and properly account for dispersing individuals in a mark-recapture study, the area covered for recoveries or resightings should be larger than the area within which individuals were marked (Anders & Marshall 2005). Otherwise, dispersing (emigrating) individuals would not be observed and would be assumed dead, biasing (underestimating) the survival rate. Other solutions include correcting based on estimated or independenly available population trend (Ryu et al. 2016), obtaining known-fate data from radio-tagged individuals (i.e., radiotelemetry data) (Pollock et al. 1989, 1995, 2004; Conroy et al. 1996; Powell et al. 2000; Nasution et al. 2001), and using a combination of band recoveries (dead recoveries) and live recaptures to estimate true survival (Barker & White 2002).
Uncertainty in survival rate estimation
Reducing uncertainty in survival rate estimation often means increasing the sample size, i.e., marking or censusing a larger number of individuals. However, in cases where individuals can be selected for marking (esp. in plant population studies), the distribution of marked individuals among stages can also determine the amount of sampling error or measurement uncertainty. In most studies, such sampling is random; i.e., on average a similar proportion of individuals in each stage are marked (for example, when sampling is based on fixed plots). Two alternatives to this design may allow more precise estimates of survival rates:
The first alternative is to determine the number of individuals to be sampled in each stage, based on how much the stages contribute to the growth rate of the population (Gross 2002). This method requires an "educated guess" of the stage matrix, based on prior studies or preliminary data.
The second alternative is to sample an equal number of individuals from each stage. A simulation study by Munzbergova and Ehrlen (2005) demonstrated that sampling an equal number of individuals per stage overall provides better model predictions than plot-based sampling. It is also better than the method proposed by Gross (2002), unless relatively good estimates of the stage matrix are available before the study.
Using either of these alternatives may be more difficult than a plot-based design because of practical reasons (e.g., ease of relocating individuals) and because they require a priori definition of stages (Munzbergova and Ehrlen 2005).
Using survivorship instead of survival rate
Life tables often include survivorship, l(x), instead of survival rate, S(x). Survivorship is the proportion of the original number of individuals in the cohort that are still alive at the beginning of age x (by definition, l(0) = 1.0). For a stage matrix, what is needed is survival rate, which is the probability of surviving from a given age to the next, whereas survivorship is the probability of surviving from birth to a given age. Calculate survival rate using S(x) = l(x+1) / l(x).
Survival rates in the diagonal
When entering a Leslie matrix (for an age-structured model), the survival rates should be in the sub-diagonal, not diagonal. The exception is that if the last age class is an composite age class, the last survival rate is in the lower right corner of the matrix.
Ignoring constraints
A model must have a number of restrictions or constraints on survival rates and the number of survivors:
1. A survival rate must be between 0 and 1.
2. The sum of all survival transitions from a given stage must be less than 1 in any time step.
3. The sum of the number of survivors from a given stage must be less than or equal to the number of individuals in that stage in the previous time step.
In RAMAS Metapop, the "Constraints Matrix" (accessed from the Stage matrix dialog) identifies the survival rates (as opposed to fecundities), and thus allows the program to impose these three types of restrictions. If this matrix is all zeros (thus assuming all elements are fecundities), or the "Ignore Constraints" box (in General information dialog) is checked, these three types of checks will not be done, and the model will give erroneous results. This is especially important if the model includes any stochasticity (variability) in survival rates.
If you are using programs other than RAMAS Metapop/GIS, make sure that the model does impose such constraints on survival rates, as well as constraining fecundities to be non-negative.
Too many (or too few) age classes or stages
The number of age classes or stages should be consistent with available data, with what is known about the life history of the species, and with the question being addressed (see Model too complex and Model too simple). If there are too many age classes or stages, the sample size for estimating each survival rate and fecundity will be small, increasing error variance (measurement error). If there are too few age classes or stages, the survival rate and fecundity for each stage may not be uniform with that class (because each stage will include a large proportion of the population and therefore a wide variety of individual demographic traits).
On the one hand, it is necessary to define a sufficiently large number of stages so that the demographic characteristics of individuals within a given stage are similar. On the other hand, it is necessary to have a sufficiently large number of individuals in each stage so that the transition probabilities can be calculated with reasonable accuracy (also see Uncertainty above). For a discussion of this trade-off, and methods of determining the appropriate number of stages, see Moloney (1986) and Vandermeer (1978).
For an age-structured model, you can reduce the number of age classes either by pooling data from older age classes and creating a composite age class (see Model too complex), or by defining multi-year age classes. For very long-lived species, 2-year, 5-year, or even 10-year age classes may be defined. If you have, say, 2-year age classes, you have to calculate survival rates, fecundities, etc. over 2 years; and each time step of the simulation will be 2 years (so, for a 50-year simulation, you'd set Duration=25).
Modeling only females
Most models of vertebrate populations include only females. In most cases, this is fine, as long as sex ratio is properly incorporated into fecundity estimates. However, there are some cases where it is necessary to model both males and females, for example by developing matrix models with different stages or age classes for males and females. If, for example, the purpose of building a model is to evaluate the consequences of different hunting regimes, and only males (and perhaps only males over a certain age) are hunted, then the model obviously needs to have both male and female age classes (e.g., see Sezen et al. 2004). Regardless of the model objective, if males have higher mortality than females (causing a skewed sex ratio), and the mating system is (or is close to) monogamous, then the number of breeding females (and thus, overall population productivity) may be limited by the availability of males; in such cases, the model should include the males in separate stages.
The stage matrix of such a model may look like the figure below.
In this example, there are 3 age classes of females, followed by 3 age classes of males. There are 2 sets of age-specific fecundities. One row of fecundities (highlighted yellow) represents female offspring (daughters per female), and the other row (highlighted blue) represents male offspring (sons per female). The upper left quadrant also includes female survival rates; this quadrant is also the 3x3 matrix of the female-only version of this model. The lower right quadrant includes male survival rates. To complete this model, you must also specify the mating system, and the degree of polygamy (see the program manual or the help file for details).
Ignoring genetics
Most population models do not incorporate genetic factors (such as inbreeding depression). In many cases, this is justified, because data on genetic effects are rarely available, and incorporating such effects usually makes the model too complex for the available data. However, genetic factors may be important, for some populations. Inbreeding depression, for example, may strongly affect populations of naturally outbreeding diploid species, if they are not rapidly declining, if they exhibit large variations in population size, and if the simulated time horizon is long (Brook et al. 2002).
If you do not have sufficient data to model inbreeding depression, consider estimating risk of decline (rather than total extinction). If you do have sufficient data, read about modeling inbreeding depression in the program manual or the help file (RAMAS Metapop/GIS version 4 or later).
Modeling Density Dependence
Using (the wrong type of) density dependence
Density dependence can have important effects on extinction risks (e.g., see Ginzburg et al. 1990), therefore it is necessary to use caution in selecting the type of density dependence, and specifying its parameters. For example, some types of density dependence assume that there is an equilibrium population abundance, K, and if the population size is below K, the population will have an average growth rate above 1 (will tend to grow). If this is not the case, in other words, if the population is subject to "systemic" pressure (Shaffer 1981) and declining deterministically, then assuming density-dependent regulation will cause an underestimation of the risks of decline and extinction.
Thus, if the population has a long-term declining trend, or the long-term average of survival rates and fecundities give a finite rate of population change (eigenvalue) of less than 1, then it may not be appropriate to use some types of density dependence (e.g., Scramble, Contest, Ricker, Beverton-Holt, logistic, etc.), depending on how exactly they are modeled and parameterized (however, there are some exceptions to this; see below). In such cases, a density-independent ("exponential") model, a Ceiling type of density dependence, or a user-defined density dependence model may be more appropriate.
Not using density dependence
Many species are threatened by habitat loss, which is most commonly (and perhaps most naturally) modeled as a decline in the carrying capacity of the population. Thus, models that do not include such a habitat-related parameter (e.g., most density-independent models) may not be appropriate for modeling threats such as habitat loss.
(Note that simulating habitat loss may require modifying Rmax or the stage matrix, in addition to the carrying capacity).
As discussed above, an observed population decline may indicate that the population is not under density-dependent regulation. However, there are some exceptions to this. In the following situations, using a density-dependent model may be appropriate, even though the population may be declining:
1. The population may be above the carrying capacity temporarily (for example due to recent environmental fluctuations or excess migration), and declining towards it. If the population was above the carrying capacity because the carrying capacity temporarily declined due to recent environmental fluctuations, then you should also use "Standard deviation of K" parameter to model stochasticity.
2. The population may be subject to Allee effects, and declining even though population is below the carrying capacity. In this case, you should specify a density-dependent model that includes Allee effects, or use a high extinction threshold.
3. The survival rates and fecundities (the stage matrix) may have been estimated during years at which the population happened to be declining due to environmental factors, and you have reason to believe that in the long run the population is stabilized by density-dependent factors.
4. The population may be declining because of habitat loss, or decreased habitat quality, and you have reason to believe that if habitat loss stopped, the population would be stabilized by density-dependent factors. In this case, you should use the "Temporal change in K" parameter to model habitat loss.
Conversely, if a population is observed to be increasing exponentially, this does not mean that a density-independent model is the most appropriate. The population may approach or reach its carrying capacity in the (future) simulated time horizon, even if currently it is far below the carrying capacity.
Not using density dependence may cause an underestimation of risks of decline or extinction (or an overestimation of risks of spread of an invasive species), if the population is under the influence of Allee effects ("positive density dependence"). Often there are not sufficient data at low population densities, so it is difficult to model Allee effects. In such cases, focusing on the risk of decline to a sufficiently high threshold abundance may be appropriate.
Bias in Rmax estimation
Fitting a density dependent model to a time series data (count of total population size over time) may introduce a bias. Suppose you have a time series of population size estimates, N(1), N(2), N(3), etc.; and for each time step t, you calculate the growth rate as R(t) = N(t+1) / N(t).
To estimate density dependence, you may want to perform a regression of R(t) [or r(t)=log(R(t))] on N(t), and use the y-intercept as the estimate of Rmax [or rmax], assuming it is a declining function.
However, note that N(t) appears both in the dependent and in the independent variable of the regression! This implies that if N(t) are measured with error (as they are in most cases), or are subject to other (stochastic) factors, then the estimate of the slope of the regression and, hence the estimate of Rmax, will be biased. In particular, if the data are from density-independent population fluctuations, you will often get an estimate of Rmax above 1.0, meaning that you will detect density dependence even though it does not exist.
For example, the left figure below shows a time series produced by a density-independent, age-structured model with environmental stochasticity in survival rates and fecundities. The figure on the right shows an attempt to fit a simple density dependence model to this data. The line shows linear regression of population growth rate, ln(R), on population size, N. The population growth rate is calculated as ln(R(t)) = ln( N(t+1) / N(t)). Even though the model that produced the data did not include density dependence, the negative slope of the regression line indicates density dependence, with an intercept of about 0.24, which corresponds to Rmax of about 1.3.
Detecting and parameterizing density dependence is a complicated problem. For a discussion of complexities inherent in detecting and estimating density dependence, see Ryu et al. (2016), Langton et al. (2002), Lande et al. (2002), Saether et al. (2002), as well as Hassell (1986), Hassell et al. (1989), Solow (1990), and Walters (1985).
Time series of population size (N), produced by a density independent, age-structured, stochastic model.
Fitting a density-dependent model to the time series on the left. The negative slope indicates density dependence, even though the time series was created by a density-independent model.
Underestimating maximum growth rate (Rmax)
Rmax is defined as the maximal growth rate in the absence of density effects, namely at low population sizes. Thus, the observed growth rates at higher population sizes, especially the population growth rate when the population is approaching its equilibrium (carrying capacity) often underestimates Rmax. For example, if the stage matrix is estimated from a population approaching carrying capacity, and the eigenvalue of the matrix is used as Rmax, then the risks of decline and extinction may be overestimated.
A similar mistake occurs when using Ceiling-type density dependence (which does not use Rmax). Ceiling-type density dependence assumes that the vital rates (the stage matrix) are not affected by the effects oif density until the population reaches the Ceiling level (K). If the vital rates are actually affected by density before the population reaches K, and the stage matrix is estimated from a population close to K, then this type of density dependence may overestimate the risks of extinction and decline.
Overestimating maximum growth rate (Rmax)
Although Rmax is defined as the maximal growth rate in the absence of density effects, merely observing a growth rate above 1.0 at low population densities does not justify using scramble- or contest-type density dependence (incl. logistic, Ricker, Beverton-Holt, theta logistic, etc.). You either have to fit an entire dataset as discussed above, or you need other evidence which shows that the population is regulated by these types of density dependence. This is because the observed growth rates are affected by factors other than density, such as stochasticity. Populations that are not regulated by scramble- or contest-type density dependence (for example those that are only subject to ceiling-type density dependence) will also frequently experience periods of positive growth, some of which will coincide with low population sizes. If you model such a population with a scramble- or contest-type density dependence, you may underestimate extinction risks, because of the stabilizing effect of these functions (see Ginzburg et al. 1990).
Similar cautions apply to the estimation of Rmax in cases where the use of scramble- or contest-type density dependence is justified. Simply using the maximum growth rate ever observed (or observed under laboratory conditions, etc.) may overestimate the strength of density dependence, because this observed value might have been the result of factors other than the low density of the population or the lack of competition. If you have census (count data), consider using the statistical methods referenced above.
Another way of overestimating Rmax is calculating its value for a long time step based on measurements at a short time step, especially if the long time step is longer than the generation time (or, even if it is shorter than generation time, but longer than period between reproduction events). For example, suppose that an insect species is observed to increase by two-fold in a generation (thus, Rmax is 2.0 per generation), and assume that there are 3 generations per year. Now, suppose you build a model with an annual time step. If the population can grow two-fold in one generation, it can grow 8-fold in three generations, if there is no density dependence. However, if you add density dependence to this model, it would be wrong to estimate Rmax as 8.0, because the effects of density would be felt in the population at a generation time scale, before the population grows by 8-fold. Thus, in this case it would be wrong to have a density dependent model with an annual time step. The correct model in this case would have a time step of 1/3 years and an Rmax of 2.0.
Incorrectly modeling impact under density dependence
When developing models for impact assessment, the interactions between density dependence and the simulated impact must be carefully considered; and the parameters of the model that need to be modified to simulate the impact must be carefully selected.
One common mistake is to modify only the stage matrix elements (survival rates or fecundity) under Scramble or Contest type of density dependence. If density dependence type is Scramble or Contest, the program modifies the stage matrix elements during a simulation, according to the population size at each time step and the parameters of the density dependence function (see the help file and manual for details). Because of this, when density dependence type is Scramble or Contest, changing only the stage matrix elements to model pollution (or any other impact that affects vital rates) is not adequate (it would not produce any substantial difference in dynamics). The solution is to change both stage matrix elements and the density dependence function.
Thus, modeling impacts on survival and fecundity under density dependence means changing either or both the maximum growth rate (Rmax) and carrying capacity (K). Which one should be modified to simulate impact depends on the interaction between the impact and the density dependence relationship. Changing only Rmaxor only K may not adequately describe what happens to the density dependence relationship under impact (see Figure 1).
For example, changing only Rmax assumes that there is no impact on the population if the population size (N) is close to carrying capacity (K), and that the growth rate under impact is actually higher if N>K (Figure 1A). Changing only K assumes that there is no impact on the population if the population size is small (see Figure 1B).
A. Impact (red): Rmax = 1.3; Baseline (blue): Rmax = 1.5
B. Impact (red): K = 60; Baseline (blue): K=100
Figure 1: The effect of changing only Rmax (A, left) and only K (B, right) to simulate impact under density dependence. In both graphs, the blue curve shows the baseline model, and the red curve shows the impact model. The baseline model has Rmax = 1.5, K = 100.
While the assumptions above may be valid in some situations, simulating impact in many cases requires modifying both the maximum growth rate (Rmax) and carrying capacity (K), as shown in Figure 2. The impact may reduce the growth rate of the population by the same amount regardless of population size, making the density dependence functions of baseline and impact models parallel (Figure 2A); or the effect may be stronger at high (Figure 2B) or low (Figure 2C) population sizes.
A. Impact (red): Rmax = 1.3, K = 55
B. Impact (red): Rmax = 1.3, K = 35
C. Impact (red): Rmax = 1.1, K = 35
Figure 2: The effect of changing both Rmax and K to simulate impact under density dependence.
In all three graphs, the blue curve shows the baseline model, and the red curve shows the impact model. The baseline model has Rmax=1.5; K=100.
Toxicity vs. harvest under density dependence
Another mistake is to model mortality due to toxicity as harvest. These two types of mortality are not equivalent; there is a fundamental difference between them under density dependence. When a population regulated by density dependence is harvested, the remaining individuals may be better off than before the impact (harvest) because there may be more resources per individual. This is called compensation. The situation under toxicity is entirely different: the remaining individuals may be few, and there may be more resources per individual, but they are still poisoned. Most likely, they have higher probability of mortality, lower reproduction, and slower growth; in other words, they are probably not better off than before the impact (pollution).
In most cases, harvest by itself does not change the density dependence relationship, and compensation effects are represented by the density dependence function. Toxicity, however, most likely shifts the density dependence curve downward, as in Figure 2 above.
Note that there are limits to compensation even with harvest. Depending on its rate and timing, harvest can drive a population to extinction, even if the population is regulated by density dependence. Harvest may also interact with other factors. For example, Allee effects and stochasticity may cause a population reduced by harvest to decline even further or increase its risk of extinction.
Adding Variability (Stochasticity)
Ignoring variability
Although deterministic measures such as predicted population size and population growth rate may be useful under certain circumstances, variability is such a pervasive feature of natural populations that ignoring it almost always leads to an incomplete picture of the state of the population and an incorrect prediction of its future.
Deterministic models make point (single-value) estimates of variables such as future population size, based on point estimates (e.g., average values) of parameters such as survival rates, fecundities, etc. In natural populations, these parameters show variation in time and space, and it is not the average values of these parameters that put the populations at risk; it is the lower values they take (or the tails of their distributions). Therefore, ignoring variability often gives a misleading, optimistic picture of a population's viability.
In RAMAS Metapop, natural variability in parameter values (environmental stochasticity) is modeled by random fluctuations in age or stage-specific fecundities and survivor rates, carrying capacities, and dispersal rates, as well as two types of catastrophes. In addition, demographic stochasticity and observation error can be simulated (see the Stochasticity dialog under the Model menu).
Not using demographic stochasticity
Demographic stochasticity is the sampling variation in the number of survivors and the number of offspring that occurs (even if survival rates and fecundities were constant) because a population is made up of a finite, integer number of individuals. In RAMAS Metapop, demographic stochasticity is modeled by sampling the number of survivors from a binomial distribution and the number of offspring from a Poisson distribution (Akçakaya 1991). Relative to other factors, demographic stochasticity becomes more important at small population sizes.
Demographic stochasticity option should be used for all models, unless you are modeling densities (such as number of animals per km-sq) instead of absolute numbers of individuals (however, in this case, note the fact that the program always models abundance as an integer). It is especially important to model demographic stochasticity when modeling impacts of habitat fragmentation. A model that ignores demographic stochasticity will underestimate the increased extinction risk when a population is fragmented into several smaller populations.
If the standard deviation estimates you are using incorporate the effects of demographic variability, it is more appropriate to estimate standard deviation without the contribution by demographic stochasticity (see below), than to exclude demographic stochasticity from your model. This is because if the population size decreases in the future, the component of observed variance due to demographic stochasticity in your data will underestimate the variance due to demographic stochasticity in these lower population sizes, thus your model will underestimate the risk of decline or extinction of the population.
Environmental stochasticity distribution
If the standard deviations are high, the results may be biased because of truncation. In this case, selecting a lognormal distribution instead of a normal distribution may be helpful. Lognormal distribution is recommended if (i) any survival rate or fecundity has a small average value with a large standard deviation, (ii) any survival rate has a high average value with a large standard deviation. For a more detailed discussion, see the help file and the manual.
Ignoring correlations
Correlations among vital rates (elements of the stage matrix) increases the variability of the population size, and hence increases the risks of decline or increase. Thus, when correlations are not known, assuming full correlation rather than independence gives results that are more conservative (precautionary). Note that the correlation discussed here is the correlation (across years) among the underlying vital rates, for example among the survival rates (or fecundities) of different age classes; it is not the correlation in the observed number or proportion of survivors (or offspring, in the case of fecundities). The distinction relates to demographic stochasticity: Even if the underlying survival rates have high correlation across years (e.g., between one-year old and two-year old individuals), the observed proportion of survivors may have a low correlation, because the observed proportion will include sampling variation (and demographic stochasticity), which by definition is independent between age classes. Thus, the observed correlation may underestimate the actual correlation if the sample size (or abundance in the observed age classes) is low. (RAMAS Metapop assumes full correlation between survival rates of different stages, and between fecundities of different stages. By default, the program sets full correlation between survival rates, fecundities and carrying capacities, which can be changed by the user).
A related issue is the correlation among vital rates of different populations in a metapopulation model, which may have a large effect on simulation results.
Too few replications
Models with environmental stochasticity should have a sufficient number of replications to adequately simulate the tails of the distributions of random variables, and to estimate risks with sufficient precision. Use the maximum number of replications (10000) unless you are running a test simulation or making a demonstration. With 1000 replications, the risk curves have 95% confidence interval of about ±0.03 (see the figure and discussion in Chapter 11 of the manual on confidence intervals for risk curves).
However, one should also be careful not to be too confident because there are a large number of replications. The confidence limits are purely statistical; they do not tell anything about the reliability of your model. If the model has large uncertainties, the risk results will not be reliable, regardless of the number replications.
Random fluctuations vs. regular oscillations
In some cases, the observed variability in a population is part of a regular oscillation, or other fluctuation that has periods longer than the time step of the model. Examples include multi-annual predator-prey cycles, and long-term oscillations in demographic parameters caused by global cycles such as ENSO. In many cases, annual fluctuations would be superimposed on such cycles in demographic rates (see the figure below for an example). When observed variability in a demographic rate includes a long-term (multi-annual) cycle, using the total observed variance as the variance of a short-term (annual) fluctuations in a model could give biased results. Note that the cycles and fluctuations we are interested in here are in vital rates (survival and fecundity); they are not cycles or fluctuations in population size or density. This is an important distinction (see below).
The correct way of incorporating an observed pattern of random fluctuations superimposed on longer-term cycles is to model the long-term cycle as a temporal trend in vital rate (using .SCH and .FCH files in the Populations dialog), and to use the remaining (residual) variance to model short-term (annual) fluctuations (environmental stochasticity as modeled in the "Standard deviations" dialog). For an example of this approach used to model predator-induced cyclicity, see Akçakaya et al. (2004b). See the program manual and help files for details about temporal trends files.
This issue is related to the temporal autocorrelation in survival and fecundity, i.e., whether deviations in a vital rate are correlated in time (e.g., "bad" years more likely to be followed by "bad" years). Population sizes are often temporally autocorrelated, because population size at any time step is (at least in part) a function of the population size in the previous time step(s). However, this does not necessarily mean that vital rates are temporally autocorrelated. In an exponential (density-independent) model, temporally uncorrelated vital rates ("white noise") will result in a "random walk" pattern of temporally autocorrelated population sizes, associated with spectral reddening. Most natural populations show dynamics "halfway" between white noise and random walk (in population size, not vital rates), a pattern that can be explained by white-noise (uncorrelated) environmental fluctuations and a combination of weak to no density dependence, age structure, and/or measurement (observation) error (see Akçakaya et al. 2003a). For these reasons, in RAMAS Metapop, temporal autocorrelation is not explicitly modeled. However, autocorrelated environmental fluctuations can be added as discussed above, by using temporal trend files in addition to environmental stochasticity in the model.
Using catastrophes instead of environmental stochasticity
Catastrophes are infrequent events that cause large reductions (or increases) in population parameters. (Events that increase vital rates are called "bonanzas", but both catastrophes and bonanzas can be modeled as Catastrophes in RAMAS Metapop, using different "multipliers".) Some models include catastrophes that would be better modeled as part of normal environmental variability. Many factors that are modeled as catastrophes (such as El Nino, fire, drought) form a continuum, and can be modeled either as catastrophes or as part of normal environmental variability (perhaps overlaid on top of a temporal trend), depending on temporal and spatial scales. The best (perhaps the only) way to determine whether to model such a factor as a catastrophe or as part of annual fluctuations is to check the distribution of vital rates (or other population parameter affected by the factor). If the distribution is bimodal (with some years having distinctly different values), then adding the factor as a catastrophe is justified (for example, see Figure 3 in Akçakaya et al. 2003b).
Overestimating catastrophe impact
If a model includes a catastrophe (e.g., a hurricane that lowers annual survival rate), then the estimates of the mean and standard deviation of the affected demographic rates (e.g., survival rate) should exclude the rates observed during catastrophe years (e.g., see Akçakaya & Atwood 1997). Otherwise, the model would include the effects of the catastrophe twice, overestimating its impact.
Not incorporating delayed effects of catastrophes
Some catastrophes have delayed effects. For example, the effects of a fire that reduces the carrying capacity by destroying part of the available habitat will last for several years until the habitat is fully restored. During this time the carrying capacity may be gradually increasing to its pre-fire level. Similarly, the effects of a toxic spill may last for several years, during which average survival rates may be gradually increasing to their normal level, even as they fluctuate from year to year. In RAMAS Metapop, such effects can be modeled by a combination of catastrophes and temporal trends in average population parameters (such as carrying capacities or survival rate). In modeling such effects, it is important to correctly specify "Time since last catastrophe" and "Reset vital rates" parameters. See the help file or the manual for more details.
Overestimating variation
The standard deviation parameters to model environmental stochasticity should be based on observed temporal variance. Variance components due to sampling and demographic stochasticity must be subtracted from total observed variance. Otherwise, standard deviations may be overestimated, which may cause overestimated risks, as well as truncation (and, consequently, bias) in vital rates.
If survival rates are based on a mark-recapture analysis, see Gould & Nichols (1998), White et al. (2002), or the help file of Program MARK on how to remove demographic/sampling variance.
If survival or fecundity estimates are based on age-structured census data, subtract the expected demographic variance from the observed variance (Akçakaya 2002). The example used in this paper is worked out in this Excel file.
If the data are from a census of the total (or part of the) population, use methods discussed by Holmes (2001, 2004).
If you have repeated or parallel estimates of the same population, the covariance of the parallel measurements or the within-year variance estimates can be used to subtract observation error from total observed variance (see Dunning et al. 2002; and Morris & Doak 2002, Chapter 5).
Using the wrong standard deviation
The standard deviation parameters to model environmental stochasticity should be based on the temporal variation in these parameters. When such data are lacking, some modelers use estimates of spatial variation (or even measurement error or sampling variation) from a single time step. There is no reason to expect that spatial variance (or sampling variation) should be similar to temporal variation.
Another mistake is to base the estimates on the standard error of mean, rather than on standard deviation of the series of estimates of the vital rate.
When "Pool variance for survivals" option is selected in RAMAS Metapop (in the "Advanced stochasticity settings"), the standard deviations must be estimated in a specific way (see the help file).
If you are estimating fecundities as a product (e.g., maternity multiplied by zero-year-old survival rate), remember that variance of the product of two random numbers is a function of their means, variances and covariance (see Akçakaya & Raphael 1998 for an example).
Not correcting truncations
Truncation of sampled values may introduce a bias, so that the realized means may be different from the average value used as the model parameter. (Note any error messages that the program displays while running a simulation.) The most common causes of truncation (and suggestions for correcting them) are listed below.
1. The standard deviations and means you entered are not based on data. If you guessed the values of stage matrix and standard deviation matrix, you may have overestimated them. The means and standard deviations should be based on observed survival rates.
2. The standard deviations are large because the distribution of the survival rates is bimodal. This may happen, for example, if males and females in the same age class or stage have different average survival rates. In this case, it may be better to include separate stages for males and females, or to model only females.
3. Bimodality may also result from catastrophes, if the survival rates in years with catastrophes are very different than those without catastrophes, and they are combined in a single distribution. In this case, it may be better to add catastrophes to the model explicitly, and separate their effects from normal year-to-year variability.
4. The standard deviations are large because they include variability due to sampling and measurement error and demographic stochasticity, and spatial variability. In this case, the standard deviations estimates should exclude these components (see above). Also, the "Pooled variance for survivals" option (in "Advanced Stochasticity Settings") will reduce truncations.
5. The distribution for environmental variation is Normal when there are survival rates close to 0 or 1. Use lognormal distribution instead.
6. The population you are modeling does not fit the assumption of the program that all survival rates within a population are perfectly correlated. The "Negative correlation for largest survival" option (in "Advanced Stochasticity Settings") provides an alternative assumption.
Phantom individuals
As discussed elsewhere, constraints must be imposed when sampling vital rates, especially survival rates. (RAMAS Metapop does this automatically, as long as the "Constraints Matrix" is properly specified.)
Even when constraints are imposed, demographic stochasticity may make the number of individuals surviving from a given stage larger than the number in the stage in the previous time step, thus creating "phantom" individuals. This happens if there are two or more survival rates (to different stages) from a given stage, the total survival rate is close to 1.0, and the number of individuals is small. This is automatically corrected by RAMAS Metapop, but if you are using another program, make sure a similar correction is made.
Duration (simulation time horizon) too long or too short
The effect of uncertainties in the model structure and parameters on model results get compounded in time. In other words, the range of predicted outcomes expands with time, so that for long time horizons (simulation durations), the results may become too uncertain to be of any use. Long time horizons also make it difficult to justify model assumptions (e.g., that the vital rates will fluctuate around the same average values with the same variability as was observed).
In contrast a simulation time horizon that is too short may yield results that are not relevant or interesting. For example, when simulation time horizon is less than the generation time of the species, the risk of extinction or substantial decline may be very low, but this low risk may not be relevant to the management question at hand.
Also, as discussed elsewhere, impacts may be assessed as very small when simulation time horizon is too long or too short (see also Akçakaya & Sjögren-Gulve 2000).
The appropriate duration depends on the biology of the species (especially its generation time), the amount of data (especially the length of the time period from which data are available), and the question addressed (especially whether an absolute or relative prediction is required).
Confusing uncertainty and variability
Confusion about the terminology related to uncertainty and variability contributes to some modeling mistakes. As modeled in RAMAS Metapop, natural variability results from temporal and spatial environmental fluctuations, catastrophes, and demographic stochasticity. Natural variability can be modeled and translated into risk (probability) estimates using a stochastic model.
Model or parameter uncertainty results from measurement error, lack of knowledge, and model misspecification and biases. It determines model accuracy, and its effects on the uncertainty of model results increases with time. In principle, this kind of uncertainty can be reduced by collecting additional data. This type of uncertainty should be incorporated in the form of ranges (lower and upper bounds) of each model parameter, anmd results used in a sensitivity analysis (see below). For more information, see the "Measurement error and Uncertainty" topic in the RAMAS Metapop help file or manual.
Note that the sensitivity analysis module of RAMAS GIS (see Chapter 13 of the manual) allows analyzing sensitivity to one parameter at a time. A more comprehensive method, called global sensitivity analysis (GSA), considers all parameters simultaneously. A new GSA method, implemented as an R package, allows global sensitivity analyses using RAMAS (Aiello-Lammens and Akçakaya 2016). The R package, including sample data and a tutorial, is freely available at https://github.com/mlammens/demgsa.
Modeling Spatial Structure
Ignoring spatial structure
Spatial structure refers to the spatial distribution of the population, so that instead of one mixed (panmictic) population, there are several subpopulations that are relatively isolated from each other (i.e., the dispersal among the subpopulations is more restricted or limited than within a subpopulation). There are several ways models ignore spatial structure.
First, if a population is composed of multiple subpopulations, then modeling it as a single population often gives the wrong results in terms of risk of decline and extinction (Burgman et al.1993; Akçakaya 2000b). In addition, certain types of threats or impacts (such as habitat fragmentation and dispersal barriers) and certain types of conservation measures (reserve design, translocations, habitat corridors, etc.) can only be evaluated by taking account the interactions between neighboring populations. Even measures that target individual populations may need to be evaluated in a metapopulation context, because the presence of other populations may change the relative effectiveness of alternative options (e.g., see Regan & Auld 2004). Determining harvest levels (for example, fishing quotas) without taking spatial structure into account may lead to overharvest (Smedbol & Stephenson 2004).
Second, a model can include multiple populations, but ignore their locations, which determine the rate of dispersal among the populations, and the correlation or similarity of the environmental conditions they experience. Both of these factors have important effects on population viability, so it is important to correctly estimate these factors (which often requires taking the locations of the populations into account).
Third, a model can include multiple populations, their locations, and interactions among the populations, but ignore the spatial variation in within-population characteristics such as survival rate, fecundity, carrying capacity, maximum growth rate, etc. Such variation may be a function of local habitat factors (such as habitat quality, habitat area, edge:core ratio, etc.) or other local factors (such as harvest rate, predation rate, etc.). Differences in these types of within-population characteristics may have important effects on metapopulation viability, e.g., through source-sink dynamics. Such differences may also result from human impacts (such as pollution or harvest), in which case ignoring them in models developed for impact assessment may cause underestimation of impacts.
Too many (or too few) populations
When a species exists in a metapopulation (e.g., in a fragmented habitat), determining the number of discrete subpopulations of the metapopulation in an arbitrary manner (e.g., by the number of sampling locations, etc.) may lead to the wrong number of populations.
Considering the difficulty of defining a species, a much more fundamental concept, it is perhaps not surprising that the definition and delineation of a population presents problems. A biological population can be defined as a group of interbreeding (i.e., panmictic) individuals. Assuming that the distribution of a species is more-or-less continuous across parts of the landscape, the question of delineating a population can be rephrased as: how far apart must two individuals be in order to be considered to be in different populations? This depends on the movement distance, home range, or some other measure related to the possibility of interbreeding. This approach, combined with modelling and prediction of suitable habitat, is used in RAMAS GIS to delineate populations (Akçakaya et al. 1995; Akçakaya 2000b).
The main parameter used to address the above question in RAMAS GIS is the Neighborhood distance, which is used in the "Spatial Data" subprogram of RAMAS GIS to find patches in the habitat map. The value is specified in units of cell lengths, and should be consistent with the biology of the species you are modeling. This parameter represents the spatial scale at which the population can be assumed to be panmictic. If you are modeling an animal, you might think of it as the foraging distance of the species modeled. Technically, suitable cells that are separated by a distance of less than or equal to the neighborhood distance are regarded to be in the same patch. The unit of distance is one cell, and the distances are measured from the centers of the cells. Thus, a Neighborhood distance of 1.0 means that two cells that are touching only in one corner are in two different patches (because the distance between their centers is about 1.4 cells). In many cases, the resolution of the maps are high enough that the minimum value of 1.0 is not appropriate.
For more information and examples, see the User Manual, Akçakaya & Raphael (1998) and Akçakaya & Atwood (1997).
Ignoring spatial correlation
Spatial correlation refers to the similarity (synchrony) of environmental fluctuations in different parts of the landscape and, in the case of a metapopulation, in different populations. In many metapopulations, fluctuations in the environment are at least partially correlated (Liebhold et al. 2004). In these cases, models that ignore spatial correlations (assuming independence) may be misleadingly optimistic in their estimation of risks of extinction and decline (e.g. Harrison & Quinn 1989; Akçakaya & Ginzburg 1991; LaHaye et al. 1994).
This is especially important when modeling impacts of habitat fragmentation. In a model that ignores spatial correlations, when a population is fragmented into several smaller populations, the model assumes that new fragments have independent dynamics, even though they were part of the same population before fragmentation. Although the new fragments may have partially independent dynamics (e.g., because some threats such diseases, fires, etc. may become less likely to spread to the whole population), making the extreme assumption of no correlation is unrealistic and may severely underestimate the true impacts of fragmentation, even making the fragmented populations appear as being more viable than the single population they formed before fragmentation.
For modeling spatial correlation in RAMAS Metapop, select "Correlations" from the Model menu, and then press F1 for help.
Wrong dispersal rates
There are several ways dispersal rates can be incorrectly estimated. Please read the help file and the manual for details of how dispersal is modeled, and how dispersal parameters should be estimated.
One important point to keep in mind is that in RAMAS Metapop dispersal is modeled in terms of the proportion of individuals moving (successfully dispersing) from one population to another. Thus, the program does not make a distinction among dispersers that die during dispersal, those which die before they leave the source patch, and those that die after they reach the target patch (but before the next census). In other words, it assumes that dispersal mortality is incorporated into the vital rates (e.g., specified in a stage matrix). This assumption does not mean that the dispersers necessarily have the same mortality as residents; it means that mortality during dispersal is accounted for by both the vital rates and the dispersal rates. This point is a pragmatic assumption because, in most field studies, one cannot really measure the proportion of individuals leaving a patch (unless all individuals are radio-tagged and continuously monitored) but can only count those that arrive in another patch. The distinction (as to where exactly dispersers die) is usually not only difficult to quantify but is also mostly irrelevant to the estimation of extinction risk. For a simple example that illustrates this point, see the program manual or Akçakaya (2000b, page 47-48).
Symmetric dispersal rates
In some cases, assuming that the dispersal rate is the same in both directions between two populations (i.e., assuming a symmetric dispersal matrix) may be problematic. This is especially true if one population is much larger than the other. If dispersal rate is the same between a large and a small population in both directions, the number of dispersers from the large to the small population would be much larger than the number in the other direction. The large number of migrants from a large to a small population will overshoot the small population's ceiling or carrying capacity (and thus not contribute much to its persistence), whereas the small number of migrants from the small population to large population will not compensate for the number that leaves the large population (Akçakaya and Baur 1996; Akçakaya & Raphael 1998). Thus, dispersal from the large to small population will drain the large population and cause N>K in the small population, and as a result, the small population will act as a pseudo-sink.
Although this may be realistic scenario in some metapopulations, in most models, it is just an artifact of assuming symmetric dispersal between all pairs of populations.
In RAMAS GIS, this problem can be corrected in several ways:
Modify dispersal rates as proportional to the carrying capacities of the source and the target populations (see Akçakaya & Raphael 1998; page 885). There is not an automatic way of doing this, and it may be complicated when carrying capacities have a temporal trend. It can be done by carefully modifying the *.DCH files (for dispersal rates) based on the *.KCH files (for K).
Use center-to-edge distances in the Spatial Data program. When the patch area is proportional to population size (so that there is a big asymmetry in area among close-by populations), center-to-edge distances make the large-to-small patch distance longer than the small-to-large distance, resulting in asymmetrical dispersal rates, with larger dispersal rate from small to large populations than from large to small populations. This would be realistic since most of the large patch lies quite far away from the small patch, whereas the entire small patch is about the same distance from the edge of the large patch.
Change the "Density-dependent dispersal as a function of target population K" parameter to a large number (using the same or similar value for all populations). This parameter is in the "Dens. Dep." tab of the Populations dialog, and is the threshold K under which dispersal rates (into this population) decrease linearly as a function of the carrying capacity of this population. (It is available only if the "Dispersal depends on Target pop. K" parameter in the Dispersal dialog is checked.) Using a large value for this parameter would resolve the problem, because it would reduce immigration populations with small K much more than immigration into larger populations.
Too high (or too low) map resolution
When the spatial structure of the metapopulation is based on habitat maps (using the "Spatial Data" subprogram of RAMAS GIS), the map resolution (cell length or pixel size) should be consistent with the biology of the species being modeled. If the resolution is too high (cell size is too small), then the maps will be unnecessarily large and some program operations will take a very long time or may fail due to memory shortage. If, on the other hand, the resolution is too low (cell size is too large), the resulting patch structure will not reflect the spatial structure of the species' habitat and populations very precisely.
When determining the appropriate resolution, it is important to keep in mind that RAMAS is for population-level modeling, not for individual-based modeling. Thus, a patch is defined as an area that supports one subpopulation of a metapopulation, not as a territory or home range of an individual or a pair. In many cases, a cell area that is close to or slightly smaller than average territory or home range size may provide sufficient resolution.
Another guide is the Neighborhood distance parameter. You must specify this parameter in the program in units of cells (cell lengths), based on the biology of the species. If the biology of the species you are modeling suggests a neighborhood distance that is larger than about 4 or 5 cells, then the resolution of your maps may be too high (you may need to increase cell length, and decrease the number of rows and columns). If it suggests a neighborhood distance that is smaller than 1 cell, then the resolution of your maps may be too low.
For tips on how to change the resolution of maps, see the help topic titled "About Map Size and Resolution."
Presenting Results
Implicit assumptions
When presenting the results of a modeling effort, it is important to explicitly list the assumptions made. Although modelers know the assumptions of their models, and many assumptions are implicit in the description of a model, those using the model results or evaluating the model may not realize what these assumptions are.
Ignoring uncertainty
Uncertainty is a prevalent feature of ecological data, and includes natural variability (which can be modeled and translated into risk estimates using a stochastic model), as well as model or parameter uncertainty (which results from measurement error, lack of knowledge, and model misspecification and biases). If data for a particular parameter of the model are unavailable or uncertain, parameter uncertainty can be incorporated by using ranges (lower and upper bounds, instead of point estimates). Uncertainties in structure of the model can be incorporated by building multiple models (e.g., with different types of density dependence).
There are various methods of propagating such uncertainties in calculations and simulations (Regan et al. 2002; Tucker and Ferson 2003). For incorporating model and parameter uncertainty in population models, one of the simplest methods is to build best-case and worst-case models. A best-case (or optimistic) model includes a combination of the lower bounds of parameters that have a negative effect of viability (such as variation in survival rate), and upper bounds of those that have a positive effect (such as average survival rate). A worst-case or pessimistic model includes the reverse bounds. Combining the results of these two models gives a range of estimates of extinction risk and other assessment endpoints. This allows the users of the PVA results (managers, conservationists) to understand the effect of uncertain input, and to make decisions with full knowledge of the uncertainties.
The uncertainties can also be used in a sensitivity analysis. When assessing impacts, it is important to prevent uncertainties making it difficult to detect impact. For more information, see the "Measurement error and Uncertainty" topic in the RAMAS Metapop help file or manual.
Emphasizing deterministic (not probabilistic) results
Although deterministic measures such as predicted population size and population growth rate may be useful under certain circumstances, variability is such a pervasive feature of natural populations that ignoring variability almost always leads to an incomplete picture of the state of the population and an incorrect prediction of its future. When a model includes variability, deterministic results such as the future population size will be unreliable, and must be accompanied by measures of their variability (such as variance, percentiles, etc.). They may also be irrelevant. If a population is subject to variability and has a substantial risk of going extinct in the near future, the long-term estimated growth rate, or the point estimates of a future population size are not relevant. In this case, probabilistic results such as risk of decline may be more relevant.
Estimating risk of extinction rather than decline
Very small populations may be subject to Allee effects, inbreeding depression, social dysfunction, or other effects that are often difficult to model because of lack of information and uncertainties associated with these factors. Thus, predictions of models for very small population sizes are not reliable. Because of this, the results are more reliable if risk of decline (rather than total extinction) is used. Thus, results should be expressed as the probability that the population size falls to or below a critical population level (e.g., for social dysfunction or other severe effects), say, 20, 50, or 100 individuals. Estimating risks of decline rather than extinction may also allow the results to be more sensitive to effects of simulated impact.
When the model includes multiple populations, there are two different types of thresholds you can set: (1) an extinction threshold for the whole population (in Stochasticity dialog), (2) a local threshold for each population. See the help files and the manual for setting these thresholds, and the difference between the two.
Using absolute (not relative) predictions
Results are usually more reliable if they are relative (e.g., "Which option gives higher viability?"), rather than absolute (e.g., "What is the risk of extinction?"). As discussed elsewhere, relative results may not be as sensitive to uncertainties in the data, even in cases where the uncertainties in the data result in uncertainties in absolute results (e.g.., see Akçakaya & Raphael 1998).
Analyzing sensitivity to single variables
Sensitivity analysis identifies important parameters; those that, if known with a higher precision, would decrease the uncertainty in model results to the largest extent. The importance of a parameter in determining viability depends on both the range of plausible values (its current uncertainty), and practical limitations (such as cost considerations).
Another use of sensitivity analysis is determining the most effective management action. This is often done by evaluating the sensitivity of the model results to each parameter. However, most management actions cause changes in more than one parameter. For example, an effort to increase the survival of newborns affects both the first survival rate and the fecundity in a matrix model based on pre-reproductive census. It most likely affects the survival of other age classes as well. In such cases, it is better to evaluate “whole-model” sensitivity (with respect to management actions) instead of parameter-by-parameter sensitivities (see Akçakaya et al. 1999 and Akçakaya 2000a).
Impact Assessment
Failing to incorporate cumulative impacts
When different impacts affecting a population or metapopulation are assessed singly, their effects may be underestimated. For example, if a fish population is affected by both fishing and pollution, the population may be assessed to be viable under each of these threats, but it may be threatened when both factors are assessed together in one model. Although singly assessing threats is a valid method of analyzing their effects, if such analyses are not made in the context of other existing threats, they may give a biased picture of the status of the population.
One of the advantages of population modeling is the ease with which effects of multiple factors can be incorporated to assess their population-level consequences. This is possible (although more complicated) even when there is interaction among the factors or threats such that the presence of one factor changes how much the population is affected by another factor. For example, lower survival rates, carrying capacity or Rmax caused by a toxic substance may make the population more susceptible to overharvest. In such a situation, if data are available to adequately model the population dynamics under different levels of pollution (esp. when there is density dependence), then the interaction of this factor with harvest will emerge as one of the results of the model.
Selecting the wrong spatial scale (or a single scale)
When a threat factor impacts a certain area, the identification of the assessment population has important consequences for the relevance and significance of the results of impact assessment. For example, suppose that an oil spill affects a small area that is only part of a population of an affected species. The assessment can focus (1) only on the part of the population affected by the spill, (2) on the whole (biological) population, part of which is affected, or (3) on the regional metapopulation (one population of which is partly affected).
The first option will result in the largest magnitude of predicted impact (e.g., in terms of proportional reduction in the population), but the relevance of the results will depend on the context. For example, whether the species is a sensitive species used as an indicator of threats to the ecosystem, or whether it is a threatened or an an economically important species, may determine how the results at different spatial scales are interpreted. If the species is a common and widespread species, used as an indicator for human health purposes, then any effect would be relevant, even if the biological population is not affected. For example, deformations or tumors in common but sensitive aquatic species in a small area may signal a relevant impact even if there are no population-level effects. If, on the other hand, the impact is specific to the species under study, and the species is of interest by itself (not necessarily an indicator species), then larger scales may also be relevant, although the spatial scale may also be a political or sociological question (for "charismatic" species, such considerations may imply local assessments). Larger scales may be relevant especially if recovery actions (in addition to impacts) are being considered: successful recovery of an impacted population may depend on its interaction with other populations.
In many cases, it may be best to make assessments at multiple spatial scales, with at least one scale for providing the larger regional context, and a smaller scale for assessing the more local implications.
Reference population includes impact effects
Impact assessment often involves predicting population-level effects at various levels of a threatening factor. In its simplest formulation, this means comparing two scenarios: baseline (or, reference) and impact. If the model parameters for the baseline scenario are estimated from data that includes some level of the threatening factor, the impact may be overestimated.
Uncertainty masking impact
When a model includes large parameter or structural uncertainties, the resulting uncertainty in model outcome may, if not properly handled, mask the impact. For a very simple example, suppose that the risk of a 50% decline in a population is predicted as ranging from 0.1 to 0.4 under the baseline (no impact) scenario and 0.2 to 0.5 under the impact scenario. The large overlap in these risk estimates may be mistakenly interpreted to mean that significant impact cannot be demonstrated. However, this interpretation may be wrong, depending on the source of uncertainty. Suppose that uncertainty in two factors (say, average adult survival, S, and temporal variation in fecundity, VF) contribute most to the uncertainty in results, such that when these factors are held constant, there is very little uncertainty in the results. Suppose that for each extreme value of these factors, the results look as follows:
Thus, even though the results are highly uncertain, and there is a large overlap between the risks under the baseline and impact scenarios, it is clear that this particular impact is causing an additional risk of 0.1, under all combinations of the uncertain parameters. Thus, it is possible to be rather certain of the magnitude of the impact, even when it is not possible to be certain of the risk under either scenario (for a full example, see Akçakaya & Raphael 1998).
Underestimating impacts of toxicity, habitat loss, habitat fragmentation, and harvest
Several of the mistakes discussed in this page may result in underestimating impacts in general. For example, impacts may be underestimated if simulation time horizon is too long or too short, if uncertainty is not properly incorporated, if impacts are estimated from short-term experiments, if cumulative effects are ignored, or if impacts are assessed at the wrong spatial scale. In addition, there are several other issues discussed elsewhere that may lead to an underestimation of the specific impacts of toxicity, habitat loss, habitat fragmentation, and harvest.
Toxicity/Pollution
Incorrectly modeling impact under density dependence
Modeling mortality due to toxicity as harvest
Overestimating Rmax; Bias in Rmax
Ignoring spatial structure
Habitat Fragmentation
Ignoring spatial correlation
Ignoring spatial structure
Not using demographic stochasticity
Habitat Loss
Not using density dependence
Not incorporating delayed effects of catastrophes
Harvest
Using (the wrong type of) density dependence
Overestimating Rmax ; Bias in Rmax
Ignoring spatial structure
Difference between experimental and model time steps
When impacts on particular parameters (such as effect of pollution on survival rate) are estimated from short-term experiments, these data cannot be directly used in a population model with longer time steps and time horizon. This type of data comes, for example, from dose-response experiments. The experiments are usually run for short periods of time (e.g., 2 weeks), compared to a time step of the model (e.g., 1 year). Using the "response" from an experiment to model population-level impacts requires extrapolating, which depends not only on the difference between the two time scales, but also on the expected trends in the concentration of the pollutant during the longer time step (e.g., the expected mortality under the "average" concentration during a year, if the model’s time step is 1 year, or the cumulative mortality as a result of 1 year of exposure to a changing concentration of the pollutant).
Wrong time horizon or threshold
Impacts may be underestimated if simulation time horizon is too long or too short. In a very short time horizon, both baseline (non-impact; reference) and impact scenarios may yield very low risks, concealing the difference between the scenarios. For example, the risk of extinction of the spotted owl in the next 2-4 years may be very low, even under the most severe impacts. This is because the generation time of this species is about 7 years, and annual adult survival rate is relatively high (70% to 90%), so that even if mortality rate increases substantially and there is no reproduction at all, the population is unlikely to go extinct in such a short time. If the risk of extinction is low under both baseline and impact scenarios, the difference between them will also be small.
In a very long time horizon, both baseline and impact scenarios may yield very high risks (both close to 1) or very uncertain risks, making the difference between the two scenarios very small or very uncertain. In most situations, there will be a time horizon at which the difference between the scenarios is maximized, a time horizon at which the model results are most sensitive to impact. To find this time horizon, run both baseline and impact scenarios for a very long time (say, 500 years), and compare the "Time to extinction" results. For example, in the figure below, the difference between the scenarios is small at short and long time horizons, and maximum at about 200 time steps.
If the difference between the scenarios is too small at short time horizons and too uncertain at long time horizons, one solution is to focus on risk of decline (rather than extinction) at shorter time horizons. Risk of decline is the probability that a population will decline by a given percentage (or, to a given threshold). As in the case of selecting time horizons, too low or too high thresholds may be undesirable. The risk will be very small (and uncertain) for very low thresholds (for example, a threshold of zero, or extinction), so the difference between the scenarios may be very small.
This risk will be very high (close to 1 for both scenarios) for very high thresholds, again making the difference between the scenarios small. The reason is that even without an impact, the population is likely to hit a very high threshold (or decline by a small percentage) just due to natural fluctuations.
As in the case of selecting a time horizon, selecting an intermediate threshold will often maximize the difference. For example, in the figure below, the difference between the scenarios is small at both small and large thresholds, and maximum at about 30 individuals.
References
Aiello-Lammens, M.A. and H.R. Akçakaya. 2016. Using global sensitivity analysis of demographic models for ecological impact assessment. Conservation Biology (in press). [Methods and data available at https://github.com/mlammens/demgsa.]
Akçakaya, H. R. 1991. A method for simulating demographic stochasticity. Ecological Modelling 54:133-136.
Akçakaya, H. R. 2000a. Population viability analyses with demographically and spatially structured models. Ecological Bulletins 48:23-38.
Akçakaya, H. R. 2000b. Viability analyses with habitat-based metapopulation models. Population Ecology 42:45-53. (The original publication is available at http://www.springerlink.com)
Akçakaya, H. R. 2002. Estimating the variance of survival rates and fecundities. Animal Conservation 5:333-336.
Akçakaya, H. R. and J. L. Atwood. 1997. A habitat-based metapopulation model of the California Gnatcatcher. Conservation Biology 11:422-434.
Akçakaya, H. R. and L. R. Ginzburg. 1991. Ecological risk analysis for single and multiple populations. Pages 73-87 in: Species Conservation: A Population-Biological Approach. A. Seitz and V. Loeschcke, eds. Birkhauser Verlag, Basel.
Akçakaya, H.R. and M.G. Raphael. 1998. Assessing human impact despite uncertainty: viability of the northern spotted owl metapopulation in the northwestern USA.Biodiversity and Conservation 7:875-894.
Akçakaya H.R. and P. Sjögren-Gulve. 2000. Population viability analysis in conservation planning: an overview. Ecological Bulletins 48:9-21.
Akçakaya, H. R., M. A. McCarthy and J. Pearce. 1995. Linking landscape data with population viability analysis: management options for the helmeted honeyeater. Biological Conservation 73:169-176.
Akçakaya, H. R., M. Burgman and L. Ginzburg. 1999. Applied Population Ecology: principles and computer exercises using RAMAS EcoLab 2.0. Second edition. Sinauer, Sunderland, MA. 285 pp.
Akçakaya, H.R., J.M. Halley, and P. Inchausti. 2003a. Population-level mechanisms for reddened spectra in ecological time series. Journal of Animal Ecology 72: 698-702.
Akçakaya, H.R., J.L. Atwood, D. Breininger, C.T. Collins, and B. Duncan. 2003b. Metapopulation dynamics of the California least tern. Journal of Wildlife Management67:829-842.
Akçakaya, H.R., M.A. Burgman, O. Kindvall, C. Wood, P. Sjögren-Gulve, J. Hatfield, and M.A. McCarthy (editors). 2004a. Species Conservation and Management: Case Studies. Oxford University Press, New York.
Akçakaya, H.R., V.C. Radeloff, D.J. Mladenoff, and H.S. He. 2004b. Integrating landscape and metapopulation modeling approaches: viability of the sharp-tailed grouse in a dynamic landscape. Conservation Biology 18: 526-537.
Anders, A. D., and M. R. Marshall. 2005. Increasing the accuracy of productivity and survival estimates in assessing landbird population status. Conservation Biology 19:66-74.
Barker, R. J., and G. C. White. 2002. Joint analysis of live and dead encounters of marked animals. Proceedings of the 2nd International Wildlife Management Congress, Godollo, Hungary.
Brook, B.W., J. J. O'Grady, A. P. Chapman, M. A. Burgman, H. R. Akçakaya, R. Frankham. 2000. Predictive accuracy of population viability analysis in conservation biology.Nature 404:385-387.
Brook, B. W., D. W. Tonkyn, J. J. O'Grady, and R. Frankham. 2002. Contribution of inbreeding to extinction risk in threatened species. Conservation Ecology 6(1): 16. [online] URL: http://www.consecol.org/vol6/iss1/art16/
Burgman, M., S. Ferson and H. R. Akçakaya. 1993. Risk Assessment in Conservation Biology. Chapman & Hall, London. Population and Community Biology Series, 314 pp.
Caswell, H. 2001. Matrix Population Models: construction, analysis and interpretation. Sinauer, Sunderland, Massachusetts.
Conroy, M. J., J. E. Anderson, S. L. Rathbun, and D. G. Krementz. 1996. Statistical inference on patch-specific survival and movement rates from marked animals.Environmental and Ecological Statistics 3:99-116.
Dunham, A.E., H.R. Akçakaya, T. S. Bridges. 2006. Using scalar models for precautionary assessments of threatened species. Conservation Biology 20: 1499-1506.
Dunning, D., Ross, Q., Munch, S., & Ginzburg, L.R. (2002) Measurement error affects risk estimates for recruitment to the Hudson River stock of striped bass. The Scientific World Journal 2(S1), 238-253.
Ginzburg, L.R., S.F. Ferson and H. R. Akçakaya. 1990. Reconstructibility of density dependence and the conservative assessment of extinction risk. Conservation Biology4:63-70.
Ginzburg, L.R. and C.X.J. Jensen. 2004. Rules of thumb for judging ecological theories. Trends in Ecology and Evolution 19(3): 121-126.
Gould, W.R. and J.D. Nichols. 1998. Estimation of temporal variability of survival in animal populations. Ecology 79: 2531-2538.
Gross, K. 2002. Efficient data collection for estimating growth rates of structured populations. Ecology 83:1762-1767.
Harrison, S. and J.F. Quinn. 1989. Correlated environments and the persistence of metapopulations. Oikos 56:293-298.
Hassell, M.P. 1986. Detecting density dependence. Trends in Ecology and Evolution 1: 90-93.
Hassell, M.P., J. Latto and R.M. May. 1989. Seeing the wood for the trees: detecting density dependence from existing life table studies. Journal of Animal Ecology 58:883-892.
Holmes, E. E. 2001. Estimating risks in declining populations with poor data. Proceedings of the National Academy of Sciences of the United States of America 98:5072-5077.
Holmes, E. E. 2004. Beyond theory to application and evaluation: Diffusion approximations for population viability analysis. Ecological Applications 14:1272-1293.
LaHaye, W. S., R. J. Gutierrez and H. R. Akçakaya. 1994. Spotted owl metapopulation dynamics in southern California. Journal of Animal Ecology 63:775-785.
Lande et al. 2002. Estimating density dependence in time-series of age-structured populations. Phil. Trans. of the Royal Soc. B 357:1179-1184.
Langton et al. 2002. The estimation of density dependence using census data from several sites. Oecologia133:466-473.
Liebhold, A., W. D. Koenig, and O. N. Bjornstad. 2004. Spatial synchrony in population dynamics. Annual Review of Ecology Evolution and Systematics 35:467-490.
Moloney, K.A. 1986. A generalized algorithm for determining category size. Oecologia 69: 176-180.
Morris, W.F. and D.F. Doak 2002. Quantitative Conservation Biology. Sinauer, Sunderland, Massachusetts.
Munzbergova Z., and J. Ehrlen. 2005. How best to collect demographic data for population viability analysis models. Journal of Applied Ecology 42:1115-1120.
Nasution, M. D., C. Brownie, K. H. Pollock, and R. E. Bennetts. 2001. Estimating survival from joint analysis of resighting and radiotelemetry capture-recapture data for wild animals. Journal of Agricultural Biological and Environmental Statistics 6:461-478.
Pollock, K. H., C. M. Bunck, S. R. Winterstein, and C. L. Chen. 1995. A capture-recapture survival analysis model for radio-tagged animals. Journal of Applied Statistics 22:661-672.
Pollock, K. H., H. H. Jiang, and J. E. Hightower. 2004. Combining telemetry and fisheries tagging models to estimate fishing and natural mortality rates. Transactions of the American Fisheries Society 133:639-648.
Pollock, K. H., S. R. Winterstein, C. M. Bunck, and P. D. Curtis. 1989. Survival Analysis in Telemetry Studies - the Staggered Entry Design. Journal of Wildlife Management53:7-15.
Powell, L. A., M. J. Conroy, J. E. Hines, J. D. Nichols, and D. G. Krementz. 2000. Simultaneous use of mark-recapture and radiotelemetry to estimate survival, movement, and capture rates. Journal of Wildlife Management 64:302-313.
Ryu, H.Y., K.T. Shoemaker, E. Kneip, A. Pidgeon, P. Heglund, B. Bateman, W. Thogmartin, and H.R. Akçakaya. 2016. Developing population models with data from marked individuals. Biological Conservation 197:190-199. [Methods and data available at github.com/Akcakaya/MAPS-to-Models.]
Regan, H. M. and T. D. Auld. 2004. Australian shrub Grevillea caleyi: recovery through management of fire and predation. Pages 23-35 in Akçakaya, H.R., M.A. Burgman, O. Kindvall, C. Wood, P. Sjögren-Gulve, J. Hatfield, and M.A. McCarthy, editors. Species Conservation and Management: Case Studies. Oxford University Press.
Regan, H. M., M. Colyvan, and M. A. Burgman. 2002. A taxonomy and treatment of uncertainty for ecology and conservation biology. Ecological Applications 12: 618–628.
Saether et al. 2002. Stochastic population dynamics of an introduced Swiss population of the ibex. Ecology 83:3457-3465
Sezen, Z., H. R. Akçakaya, and C. C. Bilgin. 2004. Turkish Mouflon (Ovis gmelinii anatolica) in Central Anatolia: population viability under scenarios of harvesting for trophy. Pages 459-468 in Akçakaya, H.R., M.A. Burgman, O. Kindvall, C. Wood, P. Sjögren-Gulve, J. Hatfield, and M.A. McCarthy, editors. Species Conservation and Management: Case Studies. Oxford University Press.
Smedbol, R. K., and R. L. Stephenson. 2004. Atlantic Herring (Clupea harengus) in the Northwest Atlantic Ocean: dynamics of nested population components under several harvest regimes. Pages 245-255 in Akçakaya, H.R., M.A. Burgman, O. Kindvall, C. Wood, P. Sjögren-Gulve, J. Hatfield, and M.A. McCarthy, editors. Species Conservation and Management: Case Studies. Oxford University Press.
Shaffer, M.L. 1981. Minimum population sizes for species conservation. BioScience 31:131-134.
Solow, A.R. 1990. Testing for density dependence. A cautionary note. Oecologia 83:47-49.
Tucker, W. T. and S. Ferson. 2003. Probability bounds analysis in environmental risk assessments. Applied Biomathematics.
Vandermeer, J.H. 1978. Choosing category size in a stage projection matrix. Oecologia 32: 79-84.
Walters, C.J. 1985. Bias in the estimation of functional relationships from time series data. Canadian Journal of Fisheries and Aquatic Sciences 42:147-149.
Werner, P.A. and Caswell, H. 1977. Population growth rates and age versus stage-distribution models for teasel (Dipsacus sylvestris Huds.) Ecology 58:1103-1111
White, G.C., A.B. Franklin, and T.M. Shenk. 2002. Estimating parameters of PVA models from data on marked animals. Pages 169-190 in Beissinger, S.R. and D.R. Mccullough (eds). Population Viability Analysis. University of Chicago Press, Chicago.
Zens, M.S. and D.R. Peart. 2003. Dealing with death data: individual hazards, mortality and bias. Trends in Ecology and Evolution 18 (7): 366-373