Fast and flexible inference for joint models of multivariate longitudinal and survival data using integrated nested Laplace approximations.

Modeling longitudinal and survival data jointly offers many advantages such as addressing measurement error and missing data in the longitudinal processes, understanding and quantifying the association between the longitudinal markers and the survival events, and predicting the risk of events based on the longitudinal markers. A joint model involves multiple submodels (one for each longitudinal/survival outcome) usually linked together through correlated or shared random effects. Their estimation is computationally expensive (particularly due to a multidimensional integration of the likelihood over the random effects distribution) so that inference methods become rapidly intractable, and restricts applications of joint models to a small number of longitudinal markers and/or random effects. We introduce a Bayesian approximation based on the integrated nested Laplace approximation algorithm implemented in the R package R-INLA to alleviate the computational burden and allow the estimation of multivariate joint models with fewer restrictions. Our simulation studies show that R-INLA substantially reduces the computation time and the variability of the parameter estimates compared with alternative estimation strategies. We further apply the methodology to analyze five longitudinal markers (3 continuous, 1 count, 1 binary, and 16 random effects) and competing risks of death and transplantation in a clinical trial on primary biliary cholangitis. R-INLA provides a fast and reliable inference technique for applying joint models to the complex multivariate data encountered in health research.

Bayesian survival analysis with INLA.

This tutorial shows how various Bayesian survival models can be fitted using the integrated nested Laplace approximation in a clear, legible, and comprehensible manner using the INLA and INLAjoint R-packages. Such models include accelerated failure time, proportional hazards, mixture cure, competing risks, multi-state, frailty, and joint models of longitudinal and survival data, originally presented in the article "Bayesian survival analysis with BUGS." In addition, we illustrate the implementation of a new joint model for a longitudinal semicontinuous marker, recurrent events, and a terminal event. Our proposal aims to provide the reader with syntax examples for implementing survival models using a fast and accurate approximate Bayesian inferential approach.

Bayesian estimation of two-part joint models for a longitudinal semicontinuous biomarker and a terminal event with INLA: Interests for cancer clinical trial evaluation.

Two-part joint models for a longitudinal semicontinuous biomarker and a terminal event have been recently introduced based on frequentist estimation. The biomarker distribution is decomposed into a probability of positive value and the expected value among positive values. Shared random effects can represent the association structure between the biomarker and the terminal event. The computational burden increases compared to standard joint models with a single regression model for the biomarker. In this context, the frequentist estimation implemented in the R package frailtypack can be challenging for complex models (i.e., a large number of parameters and dimension of the random effects). As an alternative, we propose a Bayesian estimation of two-part joint models based on the Integrated Nested Laplace Approximation (INLA) algorithm to alleviate the computational burden and fit more complex models. Our simulation studies confirm that INLA provides accurate approximation of posterior estimates and to reduced computation time and variability of estimates compared to frailtypack in the situations considered. We contrast the Bayesian and frequentist approaches in the analysis of two randomized cancer clinical trials (GERCOR and PRIME studies), where INLA has a reduced variability for the association between the biomarker and the risk of event. Moreover, the Bayesian approach was able to characterize subgroups of patients associated with different responses to treatment in the PRIME study. Our study suggests that the Bayesian approach using the INLA algorithm enables to fit complex joint models that might be of interest in a wide range of clinical applications.

Sensitivity and identification quantification by a relative latent model complexity perturbation in Bayesian meta-analysis.

In recent years, Bayesian meta-analysis expressed by a normal-normal hierarchical model (NNHM) has been widely used for combining evidence from multiple studies. Data provided for the NNHM are frequently based on a small number of studies and on uncertain within-study standard deviation values. Despite the widespread use of Bayesian NNHM, it has always been unclear to what extent the posterior inference is impacted by the heterogeneity prior (sensitivity S ) and by the uncertainty in the within-study standard deviation values (identification I ). Thus, to answer this question, we developed a unified method to simultaneously quantify both sensitivity and identification ( S - I ) for all model parameters in a Bayesian NNHM, based on derivatives of the Bhattacharyya coefficient with respect to relative latent model complexity (RLMC) perturbations. Three case studies exemplify the applicability of the method proposed: historical data for a conventional therapy, data from which one large study is first included and then excluded, and two subgroup meta-analyses specified by their randomization status. We analyzed six scenarios, crossing three RLMC targets with two heterogeneity priors (half-normal, half-Cauchy). The results show that S - I explicitly reveals which parameters are affected by the heterogeneity prior and by the uncertainty in the within-study standard deviation values. In addition, we compare the impact of both heterogeneity priors and quantify how S - I values are affected by omitting one large study and by the randomization status. Finally, the range of applicability of S - I is extended to Bayesian NtHM. A dedicated R package facilitates automatic S - I quantification in applied Bayesian meta-analyses.

A Bayesian spatiotemporal statistical analysis of out-of-hospital cardiac arrests.

We propose a Bayesian spatiotemporal statistical model for predicting out-of-hospital cardiac arrests (OHCAs). Risk maps for Ticino, adjusted for demographic covariates, are built for explaining and forecasting the spatial distribution of OHCAs and their temporal dynamics. The occurrence intensity of the OHCA event in each area of interest, and the cardiac risk-based clustering of municipalities are efficiently estimated, through a statistical model that decomposes OHCA intensity into overall intensity, demographic fixed effects, spatially structured and unstructured random effects, time polynomial dependence, and spatiotemporal random effect. In the studied geography, time evolution and dependence on demographic features are robust over different categories of OHCAs, but with variability in their spatial and spatiotemporal structure. Two main OHCA incidence-based clusters of municipalities are identified.

Estimating multilevel regional variation in excess mortality of cancer patients using integrated nested Laplace approximation.

Models of excess mortality with random effects were used to estimate regional variation in relative or net survival of cancer patients. Statistical inference for these models based on the Markov chain Monte Carlo (MCMC) methods is computationally intensive and, therefore, not feasible for routine analyses of cancer register data. This study assessed the performance of the integrated nested Laplace approximation (INLA) in monitoring regional variation in cancer survival. Poisson regression model of excess mortality including both spatially correlated and unstructured random effects was fitted to the data of patients diagnosed with ovarian and breast cancer in Finland during 1955-2014 with follow up from 1960 through 2014 by using the period approach with five-year calendar time windows. We estimated standard deviations associated with variation (i) between hospital districts and (ii) between municipalities within hospital districts. Posterior estimates based on the INLA approach were compared to those based on the MCMC simulation. The estimates of the variation parameters were similar between the two approaches. Variation within hospital districts dominated in the total variation between municipalities. In 2000-2014, the proportion of the average variation within hospital districts was 68% (95% posterior interval: 35%-93%) and 82% (60%-98%) out of the total variation in ovarian and breast cancer, respectively. In the estimation of regional variation, the INLA approach was accurate, fast, and easy to implement by using the R-INLA package.

Bayesian bivariate meta-analysis of diagnostic test studies with interpretable priors.

In a bivariate meta-analysis, the number of diagnostic studies involved is often very low so that frequentist methods may result in problems. Using Bayesian inference is particularly attractive as informative priors that add a small amount of information can stabilise the analysis without overwhelming the data. However, Bayesian analysis is often computationally demanding and the selection of the prior for the covariance matrix of the bivariate structure is crucial with little data. The integrated nested Laplace approximations method provides an efficient solution to the computational issues by avoiding any sampling, but the important question of priors remain. We explore the penalised complexity (PC) prior framework for specifying informative priors for the variance parameters and the correlation parameter. PC priors facilitate model interpretation and hyperparameter specification as expert knowledge can be incorporated intuitively. We conduct a simulation study to compare the properties and behaviour of differently defined PC priors to currently used priors in the field. The simulation study shows that the PC prior seems beneficial for the variance parameters. The use of PC priors for the correlation parameter results in more precise estimates when specified in a sensible neighbourhood around the truth. To investigate the usage of PC priors in practice, we reanalyse a meta-analysis using the telomerase marker for the diagnosis of bladder cancer and compare the results with those obtained by other commonly used modelling approaches. Copyright © 2017 John Wiley & Sons, Ltd.

Bayesian penalized spline models for the analysis of spatio-temporal count data.

In recent years, the availability of infectious disease counts in time and space has increased, and consequently, there has been renewed interest in model formulation for such data. In this paper, we describe a model that was motivated by the need to analyze hand, foot, and mouth disease surveillance data in China. The data are aggregated by geographical areas and by week, with the aims of the analysis being to gain insight into the space-time dynamics and to make short-term predictions, which will aid in the implementation of public health campaigns in those areas with a large predicted disease burden. The model we develop decomposes disease-risk into marginal spatial and temporal components and a space-time interaction piece. The latter is the crucial element, and we use a tensor product spline model with a Markov random field prior on the coefficients of the basis functions. The model can be formulated as a Gaussian Markov random field and so fast computation can be carried out using the integrated nested Laplace approximation approach. A simulation study shows that the model can pick up complex space-time structure and our analysis of hand, foot, and mouth disease data in the central north region of China provides new insights into the dynamics of the disease.

A Bayesian approach to estimate the biomass of anchovies off the coast of Perú.

The Northern Humboldt Current System (NHCS) is the world's most productive ecosystem in terms of fish. In particular, the Peruvian anchovy (Engraulis ringens) is the major prey of the main top predators, like seabirds, fish, humans, and other mammals. In this context, it is important to understand the dynamics of the anchovy distribution to preserve it as well as to exploit its economic capacities. Using the data collected by the "Instituto del Mar del Perú" (IMARPE) during a scientific survey in 2005, we present a statistical analysis that has as main goals: (i) to adapt to the characteristics of the sampled data, such as spatial dependence, high proportions of zeros and big size of samples; (ii) to provide important insights on the dynamics of the anchovy population; and (iii) to propose a model for estimation and prediction of anchovy biomass in the NHCS offshore from Perú. These data were analyzed in a Bayesian framework using the integrated nested Laplace approximation (INLA) method. Further, to select the best model and to study the predictive power of each model, we performed model comparisons and predictive checks, respectively. Finally, we carried out a Bayesian spatial influence diagnostic for the preferred model.

Bayesian analysis of RNA sequencing data by estimating multiple shrinkage priors.

Next generation sequencing is quickly replacing microarrays as a technique to probe different molecular levels of the cell, such as DNA or RNA. The technology provides higher resolution, while reducing bias. RNA sequencing results in counts of RNA strands. This type of data imposes new statistical challenges. We present a novel, generic approach to model and analyze such data. Our approach aims at large flexibility of the likelihood (count) model and the regression model alike. Hence, a variety of count models is supported, such as the popular NB model, which accounts for overdispersion. In addition, complex, non-balanced designs and random effects are accommodated. Like some other methods, our method provides shrinkage of dispersion-related parameters. However, we extend it by enabling joint shrinkage of parameters, including those for which inference is desired. We argue that this is essential for Bayesian multiplicity correction. Shrinkage is effectuated by empirically estimating priors. We discuss several parametric (mixture) and non-parametric priors and develop procedures to estimate (parameters of) those. Inference is provided by means of local and Bayesian false discovery rates. We illustrate our method on several simulations and two data sets, also to compare it with other methods. Model- and data-based simulations show substantial improvements in the sensitivity at the given specificity. The data motivate the use of the ZI-NB as a powerful alternative to the NB, which results in higher detection rates for low-count data. Finally, compared with other methods, the results on small sample subsets are more reproducible when validated on their large sample complements, illustrating the importance of the type of shrinkage.

Joint quantile disease mapping with application to malaria and G6PD deficiency.

Statistical analysis based on quantile methods is more comprehensive, flexible and less sensitive to outliers when compared to mean methods. Joint disease mapping is useful for inferring correlation between different diseases. Most studies investigate this link through multiple correlated mean regressions. We propose a joint quantile regression framework for multiple diseases where different quantile levels can be considered. We are motivated by the theorized link between the presence of malaria and the gene deficiency G6PD, where medical scientists have anecdotally discovered a possible link between high levels of G6PD and lower than expected levels of malaria initially pointing towards the occurrence of G6PD inhibiting the occurrence of malaria. Thus, the need for flexible joint quantile regression in a disease mapping framework arises. Our model can be used for linear and nonlinear effects of covariates by stochastic splines since we define it as a latent Gaussian model. We perform Bayesian inference using the R integrated nested Laplace approximation, suitable even for large datasets. Finally, we illustrate the model's applicability by considering data from 21 countries, although better data are needed to prove a significant relationship. The proposed methodology offers a framework for future studies of interrelated disease phenomena.

Non-stationary Bayesian spatial model for disease mapping based on sub-regions.

This paper aims to extend the Besag model, a widely used Bayesian spatial model in disease mapping, to a non-stationary spatial model for irregular lattice-type data. The goal is to improve the model's ability to capture complex spatial dependence patterns and increase interpretability. The proposed model uses multiple precision parameters, accounting for different intensities of spatial dependence in different sub-regions. We derive a joint penalized complexity prior to the flexible local precision parameters to prevent overfitting and ensure contraction to the stationary model at a user-defined rate. The proposed methodology can be used as a basis for the development of various other non-stationary effects over other domains such as time. An accompanying R package fbesag equips the reader with the necessary tools for immediate use and application. We illustrate the novelty of the proposal by modeling the risk of dengue in Brazil, where the stationary spatial assumption fails and interesting risk profiles are estimated when accounting for spatial non-stationary. Additionally, we model different causes of death in Brazil, where we use the new model to investigate the spatial stationarity of these causes.

Sex differences in mortality among children, adolescents, and young people aged 0-24 years: a systematic assessment of national, regional, and global trends from 1990 to 2021.

BACKGROUND: Differences in mortality exist between sexes because of biological, genetic, and social factors. Sex differentials are well documented in children younger than 5 years but have not been systematically examined for ages 5-24 years. We aimed to estimate the sex ratio of mortality from birth to age 24 years and reconstruct trends in sex-specific mortality between 1990 and 2021 for 200 countries, major regions, and the world. METHODS: We compiled comprehensive databases on the mortality sex ratio (ratio of male to female mortality rates) for individuals aged 0-4 years, 5-14 years, and 15-24 years. The databases contain mortality rates from death registration systems, full birth and sibling histories from surveys, and reports on household deaths in censuses. We modelled the sex ratio of age-specific mortality as a function of the mortality in both sexes using Bayesian hierarchical time-series models. We report the levels and trends of sex ratios and estimate the expected female mortality and excess female mortality rates (the difference between the estimated female mortality and the expected female mortality) to identify countries with outlying sex ratios. FINDINGS: Globally, the mortality sex ratio was 1·13 (ie, boys were more likely to die than girls of the same age) for ages 0-4 years (90% uncertainty interval 1·11 to 1·15) in 2021. This ratio increased with age to 1·16 (1·12 to 1·20) for 5-14 years, reaching 1·65 for 15-24 years (1·52 to 1·75). In all age groups, the global sex ratio of mortality increased between 1990 and 2021, driven by faster declines in female mortality. In 2021, the probability of a newborn male reaching age 25 years was 94·1% (93·7 to 94·4), compared with 95·1% for a newborn female (94·7 to 95·3). We found a disadvantage of females versus males (compared with countries with similar total mortality) in 2021 in five countries for ages 0-4 years (Algeria, Bangladesh, Egypt, India, and Iran), one country (Suriname) for ages 5-14 years, and 13 countries for ages 15-24 years (including Bangladesh and India). We found the reverse pattern (disadvantage of males vs females compared with countries of similar total mortality) in one country in ages 0-4 years (Vietnam) and eight countries in ages 15-24 years (including Brazil and Mexico). Globally, the number of excess female deaths from birth to age 24 years was 86 563 (-6059 to 164 000) in 2021, down from 544 636 (453 982 to 633 265) in 1990. INTERPRETATION: The global sex ratio of mortality for all age groups in the first 25 years of life increased between 1990 and 2021. Targeted interventions should focus on countries with outlying sex ratios of mortality to reduce disparities due to discrimination in health care, nutrition, and violence. FUNDING: The Bill & Melinda Gates Foundation, US Agency for International Development, and King Abdullah University of Science and Technology.

Assessing the impact of a movement network on the spatiotemporal spread of infectious diseases.

Linking information on a movement network with space-time data on disease incidence is one of the key challenges in infectious disease epidemiology. In this article, we propose and compare two statistical frameworks for this purpose, namely, parameter-driven (PD) and observation-driven (OD) models. Bayesian inference in PD models is done using integrated nested Laplace approximations, while OD models can be easily fitted with existing software using maximum likelihood. The predictive performance of both formulations is assessed using proper scoring rules. As a case study, the impact of cattle trade on the spatiotemporal spread of Coxiellosis in Swiss cows, 2004-2009, is finally investigated.

Variance partitioning in spatio-temporal disease mapping models.

Bayesian disease mapping, yet if undeniably useful to describe variation in risk over time and space, comes with the hurdle of prior elicitation on hard-to-interpret random effect precision parameters. We introduce a reparametrized version of the popular spatio-temporal interaction models, based on Kronecker product intrinsic Gaussian Markov random fields, that we name the variance partitioning model. The variance partitioning model includes a mixing parameter that balances the contribution of the main and interaction effects to the total (generalized) variance and enhances interpretability. The use of a penalized complexity prior on the mixing parameter aids in coding prior information in an intuitive way. We illustrate the advantages of the variance partitioning model using two case studies.

Bayesian inference for generalized linear mixed models.

Generalized linear mixed models (GLMMs) continue to grow in popularity due to their ability to directly acknowledge multiple levels of dependency and model different data types. For small sample sizes especially, likelihood-based inference can be unreliable with variance components being particularly difficult to estimate. A Bayesian approach is appealing but has been hampered by the lack of a fast implementation, and the difficulty in specifying prior distributions with variance components again being particularly problematic. Here, we briefly review previous approaches to computation in Bayesian implementations of GLMMs and illustrate in detail, the use of integrated nested Laplace approximations in this context. We consider a number of examples, carefully specifying prior distributions on meaningful quantities in each case. The examples cover a wide range of data types including those requiring smoothing over time and a relatively complicated spline model for which we examine our prior specification in terms of the implied degrees of freedom. We conclude that Bayesian inference is now practically feasible for GLMMs and provides an attractive alternative to likelihood-based approaches such as penalized quasi-likelihood. As with likelihood-based approaches, great care is required in the analysis of clustered binary data since approximation strategies may be less accurate for such data.

Combined effects of hydrometeorological hazards and urbanisation on dengue risk in Brazil: a spatiotemporal modelling study.

BACKGROUND: Temperature and rainfall patterns are known to influence seasonal patterns of dengue transmission. However, the effect of severe drought and extremely wet conditions on the timing and intensity of dengue epidemics is poorly understood. In this study, we aimed to quantify the non-linear and delayed effects of extreme hydrometeorological hazards on dengue risk by level of urbanisation in Brazil using a spatiotemporal model. METHODS: We combined distributed lag non-linear models with a spatiotemporal Bayesian hierarchical model framework to determine the exposure-lag-response association between the relative risk (RR) of dengue and a drought severity index. We fit the model to monthly dengue case data for the 558 microregions of Brazil between January, 2001, and January, 2019, accounting for unobserved confounding factors, spatial autocorrelation, seasonality, and interannual variability. We assessed the variation in RR by level of urbanisation through an interaction between the drought severity index and urbanisation. We also assessed the effect of hydrometeorological hazards on dengue risk in areas with a high frequency of water supply shortages. FINDINGS: The dataset included 12 895 293 dengue cases reported between 2001 and 2019 in Brazil. Overall, the risk of dengue increased between 0-3 months after extremely wet conditions (maximum RR at 1 month lag 1·56 [95% CI 1·41-1·73]) and 3-5 months after drought conditions (maximum RR at 4 months lag 1·43 [1·22-1·67]). Including a linear interaction between the drought severity index and level of urbanisation improved the model fit and showed the risk of dengue was higher in more rural areas than highly urbanised areas during extremely wet conditions (maximum RR 1·77 [1·32-2·37] at 0 months lag vs maximum RR 1·58 [1·39-1·81] at 2 months lag), but higher in highly urbanised areas than rural areas after extreme drought (maximum RR 1·60 [1·33-1·92] vs 1·15 [1·08-1·22], both at 4 months lag). We also found the dengue risk following extreme drought was higher in areas that had a higher frequency of water supply shortages. INTERPRETATION: Wet conditions and extreme drought can increase the risk of dengue with different delays. The risk associated with extremely wet conditions was higher in more rural areas and the risk associated with extreme drought was exacerbated in highly urbanised areas, which have water shortages and intermittent water supply during droughts. These findings have implications for targeting mosquito control activities in poorly serviced urban areas, not only during the wet and warm season, but also during drought periods. FUNDING: Royal Society, Medical Research Council, Wellcome Trust, National Institutes of Health, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, and Conselho Nacional de Desenvolvimento Científico e Tecnológico. TRANSLATION: For the Portuguese translation of the abstract see Supplementary Materials section.

Discrete versus continuous domain models for disease mapping.

The main goal of disease mapping is to estimate disease risk and identify high-risk areas. Such analyses are hampered by the limited geographical resolution of the available data. Typically the available data are counts per spatial unit and the common approach is the Besag-York-Mollié (BYM) model. When precise geocodes are available, it is more natural to use Log-Gaussian Cox processes (LGCPs). In a simulation study mimicking childhood leukaemia incidence using actual residential locations of all children in the canton of Zürich, Switzerland, we compare the ability of these models to recover risk surfaces and identify high-risk areas. We then apply both approaches to actual data on childhood leukaemia incidence in the canton of Zürich during 1985-2015. We found that LGCPs outperform BYM models in almost all scenarios considered. Our findings suggest that there are important gains to be made from the use of LGCPs in spatial epidemiology.

New echocardiographic imaging method based on the bandwidth of the ultrasound Doppler signal with applications in blood/tissue segmentation in the left ventricle.

A new imaging method, Bandwidth Imaging, which is related to the bandwidth of the ultrasound Doppler signal is proposed as a classification function for blood and tissue signal in transthoracial echocardiography in the left ventricle. An in vivo experiment is presented, where the apparent error rate of Bandwidth Imaging is compared with the apparent error rate of Second-Harmonic Imaging on 15 healthy men. The apparent error rates are calculated from the 16 myocardial wall segments defined in [M.D., Cerqueira, N.J. Weissman, V. Dilsizian, A.K. Jacobs, S. Kaul, W.K. Laskey, D.J. Pennell, J.A. Rumberger, T. Ryan, M.S. Verani, Standardized myocardial segmentation and nomenclature for tomographic imaging of the heart, Circulation (2002) 539-542]. A hypothesis test of Bandwidth Imaging having lower apparent error rate than Second-Harmonic Imaging is proved for a p-value of 0.94 in three segments in end diastole and in one segment in end systole. When data was averaged by a structural element of five radial, three lateral and four temporal samples the numbers of segments increased to nine in end diastole and to six in end systole. This experiment indicates that Bandwidth Imaging can supply additional information for automatic border detection routines on endocardium.

A note on intrinsic conditional autoregressive models for disconnected graphs.

In this note we discuss (Gaussian) intrinsic conditional autoregressive (CAR) models for disconnected graphs, with the aim of providing practical guidelines for how these models should be defined, scaled and implemented. We show how these suggestions can be implemented in two examples, on disease mapping.