Bayesian Model Selection for Incomplete Data Using the Posterior Predictive Distribution

Summary We explore the use of a posterior predictive loss criterion for model selection for incomplete longitudinal data. We begin by identifying a property that most model selection criteria for incomplete data should consider. We then show that a straightforward extension of the Gelfand and Ghosh (1998, Biometrika, 85 , 1–11) criterion to incomplete data has two problems. First, it introduces an extra term (in addition to the goodness of fit and penalty terms) that compromises the criterion. Second, it does not satisfy the aforementioned property. We propose an alternative and explore its properties via simulations and on a real dataset and compare it to the deviance information criterion (DIC). In general, the DIC outperforms the posterior predictive criterion, but the latter criterion appears to work well overall and is very easy to compute unlike the DIC in certain classes of models for missing data.     

Comparison of Hierarchical Bayesian Models for Overdispersed Count Data using DIC and Bayes' Factors

Summary .  When replicate count data are overdispersed, it is common practice to incorporate this extra-Poisson variability by including latent parameters at the observation level. For example, the negative binomial and Poisson-lognormal (PLN) models are obtained by using gamma and lognormal latent parameters, respectively. Several recent publications have employed the deviance information criterion (DIC) to choose between these two models, with the deviance defined using the Poisson likelihood that is obtained from conditioning on these latent parameters. The results herein show that this use of DIC is inappropriate. Instead, DIC was seen to perform well if calculated using likelihood that was marginalized at the group level by integrating out the observation-level latent parameters. This group-level marginalization is explicit in the case of the negative binomial, but requires numerical integration for the PLN model. Similarly, DIC performed well to judge whether zero inflation was required when calculated using the group-marginalized form of the zero-inflated likelihood. In the context of comparing multilevel hierarchical models, the top-level DIC was obtained using likelihood that was further marginalized by additional integration over the group-level latent parameters, and the marginal densities of the models were calculated for the purpose of providing Bayes' factors. The computational viability and interpretability of these different measures is considered.     

Bayesian analysis for nonlinear regression model under skewed errors,with application in growth curves

We have considered a Bayesian approach for the nonlinear regression model by replacing the normal distribution on the error term by some skewed distributions, which account for both skewness and heavy tails or skewness alone. The type of data considered in this paper concerns repeated measurements taken in time on a set of individuals. Such multiple observations on the same individual generally produce serially correlated outcomes. Thus, additionally, our model does allow for a correlation between observations made from the same individual. We have illustrated the procedure using a data set to study the growth curves of a clinic measurement of a group of pregnant women from an obstetrics clinic in Santiago, Chile. Parameter estimation and prediction were carried out using appropriate posterior simulation schemes based in Markov Chain Monte Carlo methods. Besides the deviance information criterion (DIC) and the conditional predictive ordinate (CPO), we suggest the use of proper scoring rules based on the posterior predictive distribution for comparing models. For our data set, all these criteria chose the skew‐t model as the best model for the errors. These DIC and CPO criteria are also validated, for the model proposed here, through a simulation study. As a conclusion of this study, the DIC criterion is not trustful for this kind of complex model.     

A model of physiological and allometric factors in the self-thinning curve

A generalized self-thinning curve for plants is derived from the modified Von Bertallanfy equation. When an asymptotic relation between photosynthesis per unit of leaf area and stocking density is assumed, the self-thinning curve thus derived is also asymptotic on a log-log scale but is fitted quite well by a log-linear approximation. The model predicts that the slope of the log-linear approximation is a function of (a) photosynthetic response to density and (b) the relation between leaf area and total aboveground biomass. Intercept of the log-linear approximation is a function of these plus maximum attainable biomass, site productivity, density at which maximum photosynthesis is attained, and the nature of carbon loss within the plant community. Linkages between various parameters within the model act to reduce differences in slope and intercept for species with different life history's and physiological requirements.     

Generalized hierarchical multivariate CAR models for areal data

In the fields of medicine and public health, a common application of areal data models is the study of geographical patterns of disease. When we have several measurements recorded at each spatial location (for example, information on p>/= 2 diseases from the same population groups or regions), we need to consider multivariate areal data models in order to handle the dependence among the multivariate components as well as the spatial dependence between sites. In this article, we propose a flexible new class of generalized multivariate conditionally autoregressive (GMCAR) models for areal data, and show how it enriches the MCAR class. Our approach differs from earlier ones in that it directly specifies the joint distribution for a multivariate Markov random field (MRF) through the specification of simpler conditional and marginal models. This in turn leads to a significant reduction in the computational burden in hierarchical spatial random effect modeling, where posterior summaries are computed using Markov chain Monte Carlo (MCMC). We compare our approach with existing MCAR models in the literature via simulation, using average mean square error (AMSE) and a convenient hierarchical model selection criterion, the deviance information criterion (DIC; Spiegelhalter et al., 2002, Journal of the Royal Statistical Society, Series B64, 583-639). Finally, we offer a real-data application of our proposed GMCAR approach that models lung and esophagus cancer death rates during 1991-1998 in Minnesota counties.     

Bayesian design for minimizing prediction uncertainty in bivariate spatial responses with applications to air quality monitoring

Model-based geostatistical design involves the selection of locations to collect data to minimize an expected loss function over a set of all possible locations. The loss function is specified to reflect the aim of data collection, which, for geostatistical studies, could be to minimize the prediction uncertainty at unobserved locations. In this paper, we propose a new approach to design such studies via a loss function derived through considering the entropy about the model predictions and the parameters of the model. The approach includes a multivariate extension to generalized linear spatial models, and thus can be used to design experiments with more than one response. Unfortunately, evaluating our proposed loss function is computationally expensive so we provide an approximation such that our approach can be adopted to design realistically sized geostatistical studies. This is demonstrated through a simulated study and through designing an air quality monitoring program in Queensland, Australia. The results show that our designs remain highly efficient in achieving each experimental objective individually, providing an ideal compromise between the two objectives. Accordingly, we advocate that our approach could be adopted more generally in model-based geostatistical design.     

A New Approach for Handling Longitudinal Count Data with Zero‐Inflation and Overdispersion: Poisson Geometric Process Model

For time series of count data, correlated measurements, clustering as well as excessive zeros occur simultaneously in biomedical applications. Ignoring such effects might contribute to misleading treatment outcomes. A generalized mixture Poisson geometric process (GMPGP) model and a zero‐altered mixture Poisson geometric process (ZMPGP) model are developed from the geometric process model, which was originally developed for modelling positive continuous data and was extended to handle count data. These models are motivated by evaluating the trend development of new tumour counts for bladder cancer patients as well as by identifying useful covariates which affect the count level. The models are implemented using Bayesian method with Markov chain Monte Carlo (MCMC) algorithms and are assessed using deviance information criterion (DIC).     

Integrating variability in detection probabilities when designing wildlife surveys: a case study of amphibians from south-eastern Australia

Occupancy-based monitoring programs rely on survey data to infer presence or absence of the target species. However, species may occupy a site and go undetected, leading to erroneous inference of absence (‘false absence’). If detectability is influenced by the time of year or weather conditions, survey protocols can be adjusted to minimize the chance of false absences. In this study, detection probabilities for three amphibian species from south-eastern Australia were modelled using a Bayesian approach. For aural surveys, we compared basic models, which only included effects of survey date, duration and time of day on detection, to models including additional effects of weather. Model selection using deviance information criterion (DIC) suggested that the basic model was the most parsimonious for Crinia signifera, while models including relative humidity and water temperature were most supported for Limnodynastes dumerilii and L. tasmaniensis respectively. When predictive performance was assessed by cross validation, DIC results were largely matched for C. signifera and L. dumerilii, while models of detection for L. tasmaniensis were indistinguishable, AUC scores suggesting inadequate performance. We show how results such as these can be used to design surveys, developing protocols for individual surveys and estimating the number of surveys required under those protocols to achieve a threshold cumulative probability of detection. Conservation managers can use these models to maximize the efficiency of surveys. This will improve the accuracy of occupancy data, and reduce the risk of misdirected conservation actions resulting from false absences.     

A versatile ODE approximation to a network model for the spread of sexually transmitted diseases

 We develop a moment closure approximation (MCA) to a network model of sexually transmitted disease (STD) spread through a steady/casual partnership network. MCA has been used previously to approximate static, regular lattices, whereas application to dynamic, irregular networks is a new endeavour, and application to sociologically-motivated network models has not been attempted. Our goals are 1) to investigate issues relating to the application of moment closure approximations to dynamic and irregular networks, and 2) to understand the impact of concurrent casual partnerships on STD transmission through a population of predominantly steady monogamous partnerships. We are able to derive a moment closure approximation for a dynamic irregular network representing sexual partnership dynamics, however, we are forced to use a triple approximation due to the large error of the standard pair approximation. This example underscores the importance of doing error analysis for moment closure approximations. We also find that a small number of casual partnerships drastically increases the prevalence and rate of spread of the epidemic. Finally, although the approximation is derived for a specific network model, we can recover approximations to a broad range of network models simply by varying model parameters which control the structure of the dynamic network. Thus our moment closure approximation is very flexible in the kinds of network models it can approximate. Received: 26 August 2001 / Revised version: 15 March 2002 / Published online: 23 August 2002 C.T.B. was supported by the NSF. Key words or phrases: Moment closure approximation – Network model – Pair approximation – Sexually transmitted diseases – Steady/casual partnership network     

Drivers of growth for Atlantic cod (Gadus morhua L.) in Icelandic waters – A Bayesian approach to determine spatiotemporal variation and its causes

This study assesses spatiotemporal and sex-specific growth of Atlantic cod Gadus morhua in Icelandic waters. We use a Bayesian approach which lends itself to fitting and comparing nested models such as these. We then compare fitted parameters of these models to potential explanatory variables using a redundancy analysis (RDA) to look for drivers of growth in G. morhua. Results indicate that models that incorporate differences in growth among time, space and sex are the best-fitting models according to deviance information criterion (DIC). Results from RDA indicate that capelin Mallotus villosus recruitment and biomass is highly correlated with deviations in the von Bertalannfy growth parameter k and that L_∞ is correlated with G. morhua landings in the model that uses year to account for time-varying growth and estimated G. morhua recruitment in the model that uses cohort to account for time-varying growth.     

Is my study system good enough? A case study for identifying maternal effects

In this paper, we demonstrate how simulation studies can be used to answer questions about identifiability and consequences of omitting effects from a model. The methodology is presented through a case study where identifiability of genetic and/or individual (environmental) maternal effects is explored. Our study system is a wild house sparrow (Passer domesticus) population with known pedigree. We fit pedigree‐based (generalized) linear mixed models (animal models), with and without additive genetic and individual maternal effects, and use deviance information criterion (DIC) for choosing between these models. Pedigree and R‐code for simulations are available. For this study system, the simulation studies show that only large maternal effects can be identified. The genetic maternal effect (and similar for individual maternal effect) has to be at least half of the total genetic variance to be identified. The consequences of omitting a maternal effect when it is present are explored. Our results indicate that the total (genetic and individual) variance are accounted for. When an individual (environmental) maternal effect is omitted from the model, this only influences the estimated (direct) individual (environmental) variance. When a genetic maternal effect is omitted from the model, both (direct) genetic and (direct) individual variance estimates are overestimated.     

Modeling Harvest Intensity of Sooty Shearwater Chicks by Rakiura Māori in New Zealand

Abstract: Cultural evidence suggests that sooty shearwater (Puffinus griseus) chicks have been harvested by Rakiura Māori on islands in southern New Zealand since prehistoric times. Concerns exist that modern harvests may be impacting sooty shearwater abundance. We modeled human-related and ecological determinants of harvest (total no. of individuals harvested) of sooty shearwater chicks on 11 islands and examined the relationship between shearwater abundance and harvesting rates (chicks/hr) and harvester behavior throughout the harvesting season. Models best explaining variation in harvest between harvesting areas (manu), for both the early and late parts of the harvesting season, included harvester-days (included in all models with change in deviance information criteria [ΔDIC], ΔDIC < 8.36 and ΔDIC < 11.5, for the early and late periods, respectively). Other harvest determinants included shearwater density, size of the manu, and number of people helping harvesters (all included in the top 5 models within ΔDIC = 2.25 for the late period). Areas harvested by several families under a common-property harvesting system had higher harvest intensity for their size (24% points higher, 95% credible interval 11–36%) than those managed as an exclusive resource for one family. The slowest harvesters spent more time harvesting but on average only harvested 36% (95% credible interval 15–65%) and 34% (95% credible interval 12–63%) of the harvest taken by the fastest harvesters during the early and late periods, respectively. Our results highlight the possibility of elevated harvest intensity as the population of harvesters increases. However, our models suggested that a corresponding reduction in harvesting rate at low prey densities during the most productive period could potentially regulate harvest intensity. Future research will integrate these results into prospective shearwater demographic models to assess the utility of a range of harvesting strategies in ensuring harvest sustainability.     

Bayesian model selection for spatial capture–recapture models

A vast amount of ecological knowledge generated over the past two decades has hinged upon the ability of model selection methods to discriminate among various ecological hypotheses. The last decade has seen the rise of Bayesian hierarchical models in ecology. Consequently, commonly used tools, such as the AIC, become largely inapplicable and there appears to be no consensus about a particular model selection tool that can be universally applied. We focus on a specific class of competing Bayesian spatial capture–recapture (SCR) models and apply and evaluate some of the recommended Bayesian model selection tools: (1) Bayes Factor—using (a) Gelfand‐Dey and (b) harmonic mean methods, (2) Deviance Information Criterion (DIC), (3) Watanabe‐Akaike's Information Criterion (WAIC) and (4) posterior predictive loss criterion. In all, we evaluate 25 variants of model selection tools in our study. We evaluate these model selection tools from the standpoint of selecting the “true” model and parameter estimation. In all, we generate 120 simulated data sets using the true model and assess the frequency with which the true model is selected and how well the tool estimates N (population size), a parameter of much importance to ecologists. We find that when information content is low in the data, no particular model selection tool can be recommended to help realize, simultaneously, both the goals of model selection and parameter estimation. But, in general (when we consider both the objectives together), we recommend the use of our application of the Bayes Factor (Gelfand‐Dey with MAP approximation) for Bayesian SCR models. Our study highlights the point that although new model selection tools are emerging (e.g., WAIC) in the applied statistics literature, those tools based on sound theory even under approximation may still perform much better.     

Geographical gradients of species richness: a test of the water‐energy conjecture of Hawkins et al. (2003) using European data for five taxa

Aims We present an analysis of grid‐based species‐richness data for European plants, mammals, birds, amphibians and reptiles, designed to test the proposition of Hawkins et al. (2003a ) that the single best factor describing richness variation switches from the water regime to the energy regime in the mid‐latitudes and that the ‘breakpoint’ is related to the physiological character of the taxa. We go on to develop subregional models showing the extent to which regional model fits vary as a function of the extent of the study system, and compare the relative performance of ‘water’, ‘energy’ and ‘water–energy’ models of richness for southern, northern and pan‐European models. Location Western Europe. Methods We use atlas data comprising species range data for 187 species of mammals, 445 species of breeding birds, 58 amphibians, 91 reptiles and 2362 plant species, inserted into a c. 50 × 50 km grid cell system. We used 11 modelled climate variables, averaged for the period 1961–90. Statistical analyses were carried out using generalized additive models (GAMs), with splines simplified to a maximum of four degrees of freedom, and we tested for spatial autocorrelation using Moran's I values obtained at 10 different distance intervals. We selected favoured models on the grounds of deviance explained combined with a simple parsimony criterion, such that we selected either: (1) the best two‐variable energy, water or water–energy model, or (2) a four‐variable water–energy model, where the latter improved on the best two‐variable model by a minimum of 5% deviance explained. Results Threshold energy values, at which richness shows a transition from an increasing to a decreasing function of annual solar radiation, were identified for all taxa apart from reptiles. We found conditional support for the switch from dominance of water variables (southern models) to energy variables (northern models). Our favoured models switched between ‘water’ and ‘energy’ for mammals, and between ‘energy’ and ‘water–energy’ for birds, depending on whether we used data of pan‐European extent, southern or northern subsets. Deviance explained in our favoured models varied from 15% (birds, southern Europe) to 72% (amphibians, northern Europe), i.e. ranging from very poor to good fits with the data. Comparison with previous work indicates that our models are generally consistent with (if sometimes weaker than) previous findings. Main conclusions Our models are incomplete representations of factors influencing macro‐scale richness patterns across Europe, taking no explicit account of, for example, topographic variation, human influences or long‐term climatic variation. However, with the exception of birds, for which only the northern model attains over one‐third deviance explained, the models show that climate can account for meaningful proportions of the deviance. We find general support for considering water and energy regimes together in modelling species richness, and for the proposition that water is more limiting in southern Europe and energy in the north. Our analyses demonstrate the sensitivity of model outcomes to the geographical location and extent of the study system, illustrating that simple curve‐fitting exercises like these, particularly if based on regions with the complex history and geography characteristic of Europe, are unlikely to provide the basis for global, predictive models. However, such exercises may be of value in detecting which aspects of water and energy regimes may be of most importance in refining independently generated global models for regional application.     

The analysis of ring-recovery data using random effects

We show how random terms, describing both yearly variation and overdispersion, can easily be incorporated into models for mark-recovery data, through the use of Bayesian methods. For recovery data on lapwings, we show that the incorporation of the random terms greatly improves the goodness of fit. Omitting the random terms can lead to overestimation of the significance of weather on survival, and overoptimistic prediction intervals in simulations of future population behavior. Random effects models provide a natural way of modeling overdispersion-which is more satisfactory than the standard classical approach of scaling up all standard errors by a uniform inflation factor. We compare models by means of Bayesian p-values and the deviance information criterion (DIC).     

An integrated approach to identify distribution of Oncomelania hupensis, the intermediate host of Schistosoma japonicum, in a mountainous region in China

The aim of this study is to better understand ecological variability related to the distribution of Oncomelania hupensis, the snail intermediate host of Schistosoma japonicum, and predict the spatial distribution of O. hupensis at the local scale in order to develop a more effective control strategy for schistosomiasis in the hilly and mountainous regions of China. A two-pronged approach was applied in this study consisting of a landscape pattern analysis complemented with Bayesian spatial modelling. The parasitological data were collected by cross-sectional surveys carried out in 11 villages in 2006 and mapped based on global positioning system (GPS) coordinates. Environmental surrogates and landscape metrics were derived from remotely-sensed images and land-cover/land-use classification data. Bayesian non-spatial and spatial models were applied to investigate the variation of snail density in relation to environmental surrogates and landscape metrics at the local scale. A Bayesian spatial model, validated by the deviance information criterion (DIC), was found to be the best-fitting model. The mean shape index (MSI) and Shannon's evenness indexes (SEI) were significantly associated with snail density. These findings suggest that decreasing the heterogeneity of the landscape can reduce snail density. A prediction maps were generated by the Bayesian model together with environmental surrogates and landscape metrics. In conclusion, the risk areas of snail distribution at the local scale can be identified using an integrated approach with landscape pattern analysis supported by remote sensing and GIS technologies, as well as Bayesian modelling.     

European anchovy (<Emphasis Type="Italic">Engraulis encrasicolus</Emphasis>) landings and environmental conditions on the Catalan Coast (NW Mediterranean) during 2000–2005

Generalized additive models are proposed for a better understanding of the underlying mechanisms for anchovy variations in abundance. Environmental variables derived from satellite imagery (surface chlorophyll, sea surface temperature and wind-mixing index), river discharge (Rhône River and Ebre River) and anchovy landings (landings per unit of effort) as proxy for abundance were used, and three fishing zones were defined along the Catalan Coast. A time shift among wind index mixing, sea surface temperature and chlorophyll was observed for these variables to be significantly correlated with anchovy. Results pointed out to processes that appear to greatly influence species abundance and affect different life stages of anchovy (conditions preceding reproduction, larvae growth and survival and recruits growth). A high proportion of anchovy LPUE variability could be explained by environmental variables. Thus, some univariate models explained deviance are more than 50%, even up to around 70% of anchovy variability. In several cases the deviance explained by a given variable was even higher at the longer time-lags. Among all univariate and bivariate models fitted, the model that best explained anchovy LPUE variability, 79% of total deviance, was a model proposed for the central zone, based on the additive effect of surface chlorophyll and Rhône River discharge, considering time lags of 15 and 18 months, respectively, for each variable.