Markov chain analysis of succession in a rocky subtidal community

We present a Markov chain model of succession in a rocky subtidal community based on a long-term (1986-1994) study of subtidal invertebrates (14 species) at Ammen Rock Pinnacle in the Gulf of Maine. The model describes successional processes (disturbance, colonization, species persistence, and replacement), the equilibrium (stationary) community, and the rate of convergence. We described successional dynamics by species turnover rates, recurrence times, and the entropy of the transition matrix. We used perturbation analysis to quantify the response of diversity to successional rates and species removals. The equilibrium community was dominated by an encrusting sponge (Hymedesmia) and a bryozoan (Crisia eburnea). The equilibrium structure explained 98% of the variance in observed species frequencies. Dominant species have low probabilities of disturbance and high rates of colonization and persistence. On average, species turn over every 3.4 years. Recurrence times varied among species (7-268 years); rare species had the longest recurrence times. The community converged to equilibrium quickly (9.5 years), as measured by Dobrushin's coefficient of ergodicity. The largest changes in evenness would result from removal of the dominant sponge Hymedesmia. Subdominant species appear to increase evenness by slowing the dominance of Hymedesmia. Comparison of the subtidal community with intertidal and coral reef communities revealed that disturbance rates are an order of magnitude higher in coral reef than in rocky intertidal and subtidal communities. Colonization rates and turnover times, however, are lowest and longest in coral reefs, highest and shortest in intertidal communities, and intermediate in subtidal communities.     

Computing a population budget from sequentially sacrificed,replicated cultures

When a population budget must be obtained from censuses based on replicated, sacrificed cultures, it is difficult to obtain estimates of transition probabilities and of the errors of such estimates, because there is no logical basis for pairing successive census counts. In a study of this nature estimating a population budget of immature stages of the housefly, the problem was solved by a randomization treatment of the original census results obtained at two densities. One hundred randomly generated census matrices over all census times for each density were smoothed to remove the effects of sampling error and a population budget constructed according to defined rules. Transition probabilities computed from the population budget were plotted on triangular coordinate paper and mean probabilities, 95% confidence regions for these means, and 95% equal frequency ellipses computed. All computations and the graphing of the results were carried out on a digital computer. The computer program, available from the authors, is written in FORTRAN IV and could be easily modified for similar studies.     

Calculating confidence intervals for prediction error in microarray classification using resampling

Cross-validation based point estimates of prediction accuracy are frequently reported in microarray class prediction problems. However these point estimates can be highly variable, particularly for small sample numbers, and it would be useful to provide confidence intervals of prediction accuracy. We performed an extensive study of existing confidence interval methods and compared their performance in terms of empirical coverage and width. We developed a bootstrap case cross-validation (BCCV) resampling scheme and defined several confidence interval methods using BCCV with and without bias-correction. The widely used approach of basing confidence intervals on an independent binomial assumption of the leave-one-out cross-validation errors results in serious under-coverage of the true prediction error. Two split-sample based methods previously proposed in the literature tend to give overly conservative confidence intervals. Using BCCV resampling, the percentile confidence interval method was also found to be overly conservative without bias-correction, while the bias corrected accelerated (BCa) interval method of Efron returns substantially anti-conservative confidence intervals. We propose a simple bias reduction on the BCCV percentile interval. The method provides mildly conservative inference under all circumstances studied and outperforms the other methods in microarray applications with small to moderate sample sizes.     

Estimating rates of local extinction and colonization in colonial species and an extension to the metapopulation and community levels

Coloniality has mainly been studied from an evolutionary perspective, but relatively few studies have developed methods for modelling colony dynamics. Changes in number of colonies over time provide a useful tool for predicting and evaluating the responses of colonial species to management and to environmental disturbance. Probabilistic Markov process models have been recently used to estimate colony site dynamics using presence–absence data when all colonies are detected in sampling efforts. Here, we define and develop two general approaches for the modelling and analysis of colony dynamics for sampling situations in which all colonies are, and are not, detected. For both approaches, we develop a general probabilistic model for the data and then constrain model parameters based on various hypotheses about colony dynamics. We use Akaike's Information Criterion (AIC) to assess the adequacy of the constrained models. The models are parameterised with conditional probabilities of local colony site extinction and colonization. Presence–absence data arising from Pollock's robust capture–recapture design provide the basis for obtaining unbiased estimates of extinction, colonization, and detection probabilities when not all colonies are detected. This second approach should be particularly useful in situations where detection probabilities are heterogeneous among colony sites. The general methodology is illustrated using presence–absence data on two species of herons. Estimates of the extinction and colonization rates showed interspecific differences and strong temporal and spatial variations. We were also able to test specific predictions about colony dynamics based on ideas about habitat change and metapopulation dynamics. We recommend estimators based on probabilistic modelling for future work on colony dynamics. We also believe that this methodological framework has wide application to problems in animal ecology concerning metapopulation and community dynamics.     

Searching for convergence in phylogenetic Markov chain Monte Carlo

Markov chain Monte Carlo (MCMC) is a methodology that is gaining widespread use in the phylogenetics community and is central to phylogenetic software packages such as MrBayes. An important issue for users of MCMC methods is how to select appropriate values for adjustable parameters such as the length of the Markov chain or chains, the sampling density, the proposal mechanism, and, if Metropolis-coupled MCMC is being used, the number of heated chains and their temperatures. Although some parameter settings have been examined in detail in the literature, others are frequently chosen with more regard to computational time or personal experience with other data sets. Such choices may lead to inadequate sampling of tree space or an inefficient use of computational resources. We performed a detailed study of convergence and mixing for 70 randomly selected, putatively orthologous protein sets with different sizes and taxonomic compositions. Replicated runs from multiple random starting points permit a more rigorous assessment of convergence, and we developed two novel statistics, delta and epsilon, for this purpose. Although likelihood values invariably stabilized quickly, adequate sampling of the posterior distribution of tree topologies took considerably longer. Our results suggest that multimodality is common for data sets with 30 or more taxa and that this results in slow convergence and mixing. However, we also found that the pragmatic approach of combining data from several short, replicated runs into a "metachain" to estimate bipartition posterior probabilities provided good approximations, and that such estimates were no worse in approximating a reference posterior distribution than those obtained using a single long run of the same length as the metachain. Precision appears to be best when heated Markov chains have low temperatures, whereas chains with high temperatures appear to sample trees with high posterior probabilities only rarely.     

Estimating sampling error of evolutionary statistics based on genetic covariance matrices using maximum likelihood

We explore the estimation of uncertainty in evolutionary parameters using a recently devised approach for resampling entire additive genetic variance–covariance matrices ( G ). Large‐sample theory shows that maximum‐likelihood estimates (including restricted maximum likelihood, REML) asymptotically have a multivariate normal distribution, with covariance matrix derived from the inverse of the information matrix, and mean equal to the estimated G . This suggests that sampling estimates of G from this distribution can be used to assess the variability of estimates of G , and of functions of G . We refer to this as the REML‐MVN method. This has been implemented in the mixed‐model program WOMBAT. Estimates of sampling variances from REML‐MVN were compared to those from the parametric bootstrap and from a Bayesian Markov chain Monte Carlo (MCMC) approach (implemented in the R package MCMCglmm). We apply each approach to evolvability statistics previously estimated for a large, 20‐dimensional data set for Drosophila wings. REML‐MVN and MCMC sampling variances are close to those estimated with the parametric bootstrap. Both slightly underestimate the error in the best‐estimated aspects of the G matrix. REML analysis supports the previous conclusion that the G matrix for this population is full rank. REML‐MVN is computationally very efficient, making it an attractive alternative to both data resampling and MCMC approaches to assessing confidence in parameters of evolutionary interest.     

Transition probabilities for general birth–death processes with applications in ecology, genetics, and evolution

A birth-death process is a continuous-time Markov chain that counts the number of particles in a system over time. In the general process with n current particles, a new particle is born with instantaneous rate λ(n) and a particle dies with instantaneous rate μ(n). Currently no robust and efficient method exists to evaluate the finite-time transition probabilities in a general birth-death process with arbitrary birth and death rates. In this paper, we first revisit the theory of continued fractions to obtain expressions for the Laplace transforms of these transition probabilities and make explicit an important derivation connecting transition probabilities and continued fractions. We then develop an efficient algorithm for computing these probabilities that analyzes the error associated with approximations in the method. We demonstrate that this error-controlled method agrees with known solutions and outperforms previous approaches to computing these probabilities. Finally, we apply our novel method to several important problems in ecology, evolution, and genetics.     

Markovian analysis of cervical cell images.

Markovian analysis is a method to measure optical texture based on gray-level transition probabilities in digitized images. Experiments are described that investigate that classification performance of parameters generated by Markovian analysis. Results using Markov texture parameters show that the selection of a Markov step size strongly affects classification error rates and the number of parameters required to achieve the maximum correct classification rates. Markov texture parameters are shown to achieve high rates of correct classification in discriminating images of normal from abnormal cervical cell nuclei.     

Bayesian phylogenetic inference using DNA sequences: a Markov Chain Monte Carlo Method

An improved Bayesian method is presented for estimating phylogenetic trees using DNA sequence data. The birth-death process with species sampling is used to specify the prior distribution of phylogenies and ancestral speciation times, and the posterior probabilities of phylogenies are used to estimate the maximum posterior probability (MAP) tree. Monte Carlo integration is used to integrate over the ancestral speciation times for particular trees. A Markov Chain Monte Carlo method is used to generate the set of trees with the highest posterior probabilities. Methods are described for an empirical Bayesian analysis, in which estimates of the speciation and extinction rates are used in calculating the posterior probabilities, and a hierarchical Bayesian analysis, in which these parameters are removed from the model by an additional integration. The Markov Chain Monte Carlo method avoids the requirement of our earlier method for calculating MAP trees to sum over all possible topologies (which limited the number of taxa in an analysis to about five). The methods are applied to analyze DNA sequences for nine species of primates, and the MAP tree, which is identical to a maximum-likelihood estimate of topology, has a probability of approximately 95%.      

Estimating transition probabilities from aggregate samples plus partial transition data

Longitudinal studies often collect only aggregate data, which allows only inefficient transition probability estimates. Barring enormous aggregate samples, improving the efficiency of transition probability estimates seems to be impossible without additional partial-transition data. This paper discusses several sampling plans that collect data of both types, as well as a methodology that combines them into efficient estimates of transition probabilities. The method handles both fixed and time-dependent categorical covariates and requires no assumptions (e.g., time homogeneity, Markov) about the population evolution.     

Modeling observation error and its effects in a random walk/extinction model

This paper examines the consequences of observation errors for the "random walk with drift", a model that incorporates density independence and is frequently used in population viability analysis. Exact expressions are given for biases in estimates of the mean, variance and growth parameters under very general models for the observation errors. For other quantities, such as the finite rate of increase, and probabilities about population size in the future we provide and evaluate approximate expressions. These expressions explain the biases induced by observation error without relying exclusively on simulations, and also suggest ways to correct for observation error. A secondary contribution is a careful discussion of observation error models, presented in terms of either log-abundance or abundance. This discussion recognizes that the bias and variance in observation errors may change over time, the result of changing sampling effort or dependence on the underlying population being sampled.     

Empirical models of transitions between coral reef states: effects of region, protection, and environmental change

There has been substantial recent change in coral reef communities. To date, most analyses have focussed on static patterns or changes in single variables such as coral cover. However, little is known about how community-level changes occur at large spatial scales. Here, we develop Markov models of annual changes in coral and macroalgal cover in the Caribbean and Great Barrier Reef (GBR) regions. We analyzed reef surveys from the Caribbean and GBR (1996-2006). We defined a set of reef states distinguished by coral and macroalgal cover, and obtained Bayesian estimates of the annual probabilities of transitions between these states. The Caribbean and GBR had different transition probabilities, and therefore different rates of change in reef condition. This could be due to differences in species composition, management or the nature and extent of disturbances between these regions. We then estimated equilibrium probability distributions for reef states, and coral and macroalgal cover under constant environmental conditions. In both regions, the current distributions are close to equilibrium. In the Caribbean, coral cover is much lower and macroalgal cover is higher at equilibrium than in the GBR. We found no evidence for differences in transition probabilities between the first and second halves of our survey period, or between Caribbean reefs inside and outside marine protected areas. However, our power to detect such differences may have been low. We also examined the effects of altering transition probabilities on the community state equilibrium, along a continuum from unfavourable (e.g., increased sea surface temperature) to favourable (e.g., improved management) conditions. Both regions showed similar qualitative responses, but different patterns of uncertainty. In the Caribbean, uncertainty was greatest about effects of favourable changes, while in the GBR, we are most uncertain about effects of unfavourable changes. Our approach could be extended to provide risk analysis for management decisions.     

Hidden Markov models for settings with interval-censored transition times and uncertain time origin: application to HIV genetic analyses

The simplicity and flexibility of Markov models make them appealing for investigations of the acquisition of HIV drug-resistance mutations, whose presence can define specific Markov states. Because the exact time of acquiring a mutation is not observed during clinical research studies on HIV infection, it is important that methods for fitting such models accommodate interval-censored transition times. Furthermore, many such studies include patients with extensive treatment experience prior to the onset of the studies. Therefore, the virus in these patients may have already acquired resistance mutations by study entry. Retrospective data regarding the time on treatment, which is often known or known with error, provide information about the acquisition rates before the start of a study. Finally, variability in the genetic sequences of circulating HIV creates uncertainty in the Markov states. This paper considers approaches to fitting Markov models to data with interval-censored transition times when the time origin and the Markov states are known with error. The methods were applied to AIDS Clinical Trial Group protocol 398, a randomized comparison of mono- versus dual-protease inhibitor use in heavily pretreated patients. We found that the estimated rates of acquiring certain classes of mutations depended on the presence of others, and that the precision of these estimates can be considerably improved by inclusion of retrospective data.     

A continuous time version and a generalization of a Markov-recapture model for trapping experiments

Wileyto et al. [E.P. Wileyto, W.J. Ewens, M.A. Mullen, Markov-recapture population estimates: a tool for improving interpretation of trapping experiments, Ecology 75 (1994) 1109] propose a four-state discrete time Markov process, which describes the structure of a marking-capture experiment as a method of population estimation. They propose this method primarily for estimation of closed insect populations. Their method provides a mark-recapture estimate from a single trap observation by allowing subjects to mark themselves. The estimate of the unknown population size is based on the assumption of a closed population and a simple Markov model in which the rates of marking, capture, and recapture are assumed to be equal. Using the one step transition probability matrix of their model, we illustrate how to go from an embedded discrete time Markov process to a continuous time Markov process assuming exponentially distributed holding times. We also compute the transition probabilities after time t for the continuous time case and compare the limiting behavior of the continuous and discrete time processes. Finally, we generalize their model by relaxing the assumption of equal per capita rates for marking, capture, and recapture. Other questions about how their results change when using a continuous time Markov process are examined.     

Incorporating Genotyping Error Into Non-Invasive DNA-Based Mark—Recapture Population Estimates

ABSTRACT Use of non-invasive sources of DNA, such as hair or scat, to obtain a genetic mark for population estimates is becoming commonplace. Unfortunately, with such marks, potentials for genotyping errors and for the shadow effect have resulted in use of many loci and amplification of each specimen many times at each locus, drastically increasing time and cost of obtaining a population estimate. We proposed a method, the Genotyping Uncertainty Added Variance Adjustment (GUAVA), which statistically adjusts for genotyping errors and the shadow effect, thereby allowing use of fewer loci and one amplification of each specimen per locus. Using allele frequencies and estimates of genotyping error rates, we determined, for each pair of specimens, the probability that the pair was obtained from the same individual, whether or not their observed genotypes match. Using these probabilities, we reconstructed possible capture history matrices and used this distribution to obtain a population estimate. With simulated data, we consistently found our estimates had lower bias and smaller variance than estimates based on single amplifications in which genotyping error was ignored and that were comparable to estimates based on data free of genotyping errors. We also demonstrated the method on a fecal DNA data set from a population of red wolves (Canis rufus). The GUAVA estimate based on only one amplification genotypes compares favorably to the estimate based on consensus genotypes. A program to conduct the analysis is available from the first author for UNIX or Windows platforms. Application of GUAVA may allow for increased accuracy in population estimates at reduced cost.     

Evaluation of several methods for estimating phylogenetic trees when substitution rates differ over nucleotide sites

Several maximum likelihood and distance matrix methods for estimating phylogenetic trees from homologous DNA sequences were compared when substitution rates at sites were assumed to follow a gamma distribution. Computer simulations were performed to estimate the probabilities that various tree estimation methods recover the true tree topology. The case of four species was considered, and a few combinations of parameters were examined. Attention was applied to discriminating among different sources of error in tree reconstruction, i.e., the inconsistency of the tree estimation method, the sampling error in the estimated tree due to limited sequence length, and the sampling error in the estimated probability due to the number of simulations being limited. Compared to the least squares method based on pairwise distance estimates, the joint likelihood analysis is found to be more robust when rate variation over sites is present but ignored and an assumption is thus violated. With limited data, the likelihood method has a much higher probability of recovering the true tree and is therefore more efficient than the least squares method. The concept of statistical consistency of a tree estimation method and its implications were explored, and it is suggested that, while the efficiency (or sampling error) of a tree estimation method is a very important property, statistical consistency of the method over a wide range of, if not all, parameter values is prerequisite.     

APIS: An auto‐adaptive parentage inference software that tolerates missing parents

In the context of parentage assignment using genomic markers, key issues are genotyping errors and an absence of parent genotypes because of sampling, traceability or genotyping problems. Most likelihood‐based parentage assignment software programs require a priori estimates of genotyping errors and the proportion of missing parents to set up meaningful assignment decision rules. We present here the R package APIS, which can assign offspring to their parents without any prior information other than the offspring and parental genotypes, and a user‐defined, acceptable error rate among assigned offspring. Assignment decision rules use the distributions of average Mendelian transmission probabilities, which enable estimates of the proportion of offspring with missing parental genotypes. APIS has been compared to other software (CERVUS, VITASSIGN), on a real European seabass (Dicentrarchus labrax) single nucleotide polymorphism data set. The type I error rate (false positives) was lower with APIS than with other software, especially when parental genotypes were missing, but the true positive rate was also lower, except when the theoretical exclusion power reached 0.99999. In general, APIS provided assignments that satisfied the user‐set acceptable error rate of 1% or 5%, even when tested on simulated data with high genotyping error rates (1% or 3%) and up to 50% missing sires. Because it uses the observed distribution of Mendelian transmission probabilities, APIS is best suited to assigning parentage when numerous offspring (>200) are genotyped. We have demonstrated that APIS is an easy‐to‐use and reliable software for parentage assignment, even when up to 50% of sires are missing.     

Importance of a pilot study for non-invasive genetic sampling: genotyping errors and population size estimation in red deer

Population size information is critical for managing endangered or harvested populations. Population size can now be estimated from non-invasive genetic sampling. However, pitfalls remain such as genotyping errors (allele dropout and false alleles at microsatellite loci). To evaluate the feasibility of non-invasive sampling (e.g., for population size estimation), a pilot study is required. Here, we present a pilot study consisting of (i) a genetic step to test loci amplification and to estimate allele frequencies and genotyping error rates when using faecal DNA, and (ii) a simulation step to quantify and minimise the effects of errors on estimates of population size. The pilot study was conducted on a population of red deer in a fenced natural area of 5440 ha, in France. Twelve microsatellite loci were tested for amplification and genotyping errors. The genotyping error rates for microsatellite loci were 0–0.83 (mean=0.2) for allele dropout rates and 0–0.14 (mean=0.02) for false allele rates, comparable to rates encountered in other non-invasive studies. Simulation results suggest we must conduct 6 PCR amplifications per sample (per locus) to achieve approximately 97% correct genotypes. The 3% error rate appears to have little influence on the accuracy and precision of population size estimation. This paper illustrates the importance of conducting a pilot study (including genotyping and simulations) when using non-invasive sampling to study threatened or managed populations.     

Evaluating pedigree data. II. Identifying the cause of error in families with inconsistencies

Pedigree data can be evaluated, and subsequently corrected, by analysis of the distribution of genetic markers, taking account of the possibility of mistyping . Using a model of pedigree error developed previously, we obtained the maximum likelihood estimates of error parameters in pedigree data from Tokelau. Posterior probabilities for the possible true relationships in each family are conditional on the putative relationships and the marker data are calculated using the parameter estimates. These probabilities are used as a basis for discriminating between pedigree error and genetic marker errors in families where inconsistencies have been observed. When applied to the Tokelau data and compared with the results of retyping inconsistent families, these statistical procedures are able to discriminate between pedigree and marker error, with approximately 90% accuracy, for families with two or more offspring. The large proportion of inconsistencies inferred to be due to marker error (61%) indicates the importance of discriminating between error sources when judging the reliability of putative relationship data. Application of our model of pedigree error has proved to be an efficient way of determining and subsequently correcting sources of error in extensive pedigree data collected in large surveys.