Monte Carlo EM Research Articles

Missing data imputation in multivariate t distribution with unknown degrees of freedom using expectation maximization algorithm and its stochastic variants

Many researchers encounter the missing data problem. The phenomenon may be occasioned by data omission, non-response, death of respondents, recording errors, among others. It is important to find an appropriate data imputation technique to fill in the missing positions. In this study, the Expectation Maximization (EM) algorithm and two of its stochastic variants, stochastic EM (SEM) and Monte Carlo EM (MCEM), are employed in missing data imputation and parameter estimation in multivariate t distribution with unknown degrees of freedom. The imputation efficiencies of the three methods are then compared using mean square error (MSE) criterion. SEM yields the lowest MSE, making it the most efficient method in data imputation when the data assumes the multivariate t distribution. The algorithm’s stochastic nature enables it to avoid local saddle points and achieve global maxima; ultimately increasing its efficiency. The EM and MCEM techniques yield almost similar results. Large sample draws in the MCEM’s E-step yield more or less the same results as the deterministic EM. In parameter estimation, it is observed that the parameter estimates for EM and MCEM are relatively close to the simulated data’s maximum likelihood (ML) estimates. This is not the case in SEM, owing to the random nature of the algorithm.

Linear censored regression models with skew scale mixtures of normal distributions

A special source of difficulty in the statistical analysis is the possibility that some subjects may not have a complete observation of the response variable. Such incomplete observation of the response variable is called censoring. Censorship can occur for a variety of reasons, including limitations of measurement equipment, design of the experiment, and non-occurrence of the event of interest until the end of the study. In the presence of censoring, the dependence of the response variable on the explanatory variables can be explored through regression analysis. In this paper, we propose to examine the censorship problem in context of the class of asymmetric, i.e., we have proposed a linear regression model with censored responses based on skew scale mixtures of normal distributions. We develop a Monte Carlo EM (MCEM) algorithm to perform maximum likelihood inference of the parameters in the proposed linear censored regression models with skew scale mixtures of normal distributions. The MCEM algorithm has been discussed with an emphasis on the skew-normal, skew Student-t-normal, skew-slash and skew-contaminated normal distributions. To examine the performance of the proposed method, we present some simulation studies and analyze a real dataset.

A multistate joint model for interval-censored event-history data subject to within-unit clustering and informative missingness, with application to neurocysticercosis research.

We propose a multistate joint model to analyze interval-censored event-history data subject to within-unit clustering and nonignorable missing data. The model is motivated by a study of the neurocysticercosis (NC) cyst evolution at the cyst-level, taking into account the multiple cysts phases with intermittent missing data and loss to follow-up, as well as the intra-brain clustering of observations made on a predefined data collection schedule. Of particular interest in this study is the description of the process leading to cyst resolution, and whether this process varies by antiparasitic treatment. The model uses shared random effects to account for within-brain correlation and to explain the hidden heterogeneity governing the missing data mechanism. We developed a likelihood-based method using a Monte Carlo EM algorithm for the inference. The practical utility of the methods is illustrated using data from a randomized controlled trial on the effect of antiparasitic treatment with albendazole on NC cysts among patients from six hospitals in Ecuador. Simulation results demonstrate that the proposed methods perform well in the finite sample and misspecified models that ignore the data complexities could lead to substantial biases.

Robust bivariate random-effects model for accommodating outlying and influential studies in meta-analysis of diagnostic test accuracy studies.

Due to the inevitable inter-study correlation between test sensitivity (Se) and test specificity (Sp), mostly because of threshold variability, hierarchical or bivariate random-effects models are widely used to perform a meta-analysis of diagnostic test accuracy studies. Conventionally, these models assume that the random-effects follow the bivariate normal distribution. However, the inference made using the well-established bivariate random-effects models, when outlying and influential studies are present, may lead to misleading conclusions, since outlying or influential studies can extremely influence parameter estimates due to their disproportional weight. Therefore, we developed a new robust bivariate random-effects model that accommodates outlying and influential observations and gives robust statistical inference by down-weighting the effect of outlying and influential studies. The marginal model and the Monte Carlo expectation-maximization algorithm for our proposed model have been derived. A simulation study has been carried out to validate the proposed method and compare it against the standard methods. Regardless of the parameters varied in our simulations, the proposed model produced robust point estimates of Se and Sp compared to the standard models. Moreover, our proposed model resulted in precise estimates as it yielded the narrowest confidence intervals. The proposed model also generated a similar point and interval estimates of Se and Sp as the standard models when there are no outlying and influential studies. Two published meta-analyses have also been used to illustrate the methods.

Combining Social Media and FDA Adverse Event Reporting System to Detect Adverse Drug Reactions.

IntroductionAdverse drug reactions (ADRs) are unintended reactions caused by a drug or combination of drugs taken by a patient. The current safety surveillance system relies on spontaneous reporting systems (SRSs) and more recently on observational health data; however, ADR detection may be delayed and lack geographic diversity. The broad scope of social media conversations, such as those on Twitter, can include health-related topics. Consequently, these data could be used to detect potentially novel ADRs with less latency. Although research regarding ADR detection using social media has made progress, findings are based on single information sources, and no study has yet integrated drug safety evidence from both an SRS and Twitter.ObjectiveThe aim of this study was to combine signals from an SRS and Twitter to facilitate the detection of safety signals and compare the performance of the combined system with signals generated by individual data sources.MethodsWe extracted potential drug–ADR posts from Twitter, used Monte Carlo expectation maximization to generate drug safety signals from both the US FDA Adverse Event Reporting System and posts from Twitter, and then integrated these signals using a Bayesian hierarchical model. The results from the integrated system and two individual sources were evaluated using a reference standard derived from drug labels. Area under the receiver operating characteristics curve (AUC) was computed to measure performance.ResultsWe observed a significant improvement in the AUC of the combined system when comparing it with Twitter alone, and no improvement when comparing with the SRS alone. The AUCs ranged from 0.587 to 0.637 for the combined SRS and Twitter, from 0.525 to 0.534 for Twitter alone, and from 0.612 to 0.642 for the SRS alone. The results varied because different preprocessing procedures were applied to Twitter.ConclusionThe accuracy of signal detection using social media can be improved by combining signals with those from SRSs. However, the combined system cannot achieve better AUC performance than data from FAERS alone, which may indicate that Twitter data are not ready to be integrated into a purely data-driven combination system.

NIMBLE for Bayesian Disease Mapping

This tutorial describes the basic implementation of Bayesian hierarchical models for spatial health data using the R package nimble. To quote the nimble R description:A system for writing hierarchical statistical models largely compatible with 'BUGS' and 'JAGS', writing nimbleFunctions to operate models and do basic R-style math, and compiling both models and nimbleFunctions via custom-generated C++. 'NIMBLE' includes default methods for MCMC, particle filtering, Monte Carlo Expectation Maximization, and some other tools. The nimbleFunction system makes it easy to do things like implement new MCMC samplers from R, customize the assignment of samplers to different parts of a model from R, and compile the new samplers automatically via C++ alongside the samplers 'NIMBLE' provides.Examples of the use of the package for a small range of Bayesian Disease Mapping (BDM) models is explored and focus on different approaches to model fitting and analysis are discussed. Examples of publicly available small area health data is used throughout.

Audio-Visual Speech Enhancement Using Conditional Variational Auto-Encoders

Variational auto-encoders (VAEs) are deep generative latent variable models that can be used for learning the distribution of complex data. VAEs have been successfully used to learn a probabilistic prior over speech signals, which is then used to perform speech enhancement. One advantage of this generative approach is that it does not require pairs of clean and noisy speech signals at training. In this paper, we propose audio-visual variants of VAEs for single-channel and speaker-independent speech enhancement. We develop a conditional VAE (CVAE) where the audio speech generative process is conditioned on visual information of the lip region. At test time, the audio-visual speech generative model is combined with a noise model based on nonnegative matrix factorization, and speech enhancement relies on a Monte Carlo expectation-maximization algorithm. Experiments are conducted with the recently published NTCD-TIMIT dataset as well as the GRID corpus. The results confirm that the proposed audio-visual CVAE effectively fuses audio and visual information, and it improves the speech enhancement performance compared with the audio-only VAE model, especially when the speech signal is highly corrupted by noise. We also show that the proposed unsupervised audio-visual speech enhancement approach outperforms a state-of-the-art supervised deep learning method.

Efficient inference in state-space models through adaptive learning in online Monte Carlo expectation maximization

Expectation maximization (EM) is a technique for estimating maximum-likelihood parameters of a latent variable model given observed data by alternating between taking expectations of sufficient statistics, and maximizing the expected log likelihood. For situations where sufficient statistics are intractable, stochastic approximation EM (SAEM) is often used, which uses Monte Carlo techniques to approximate the expected log likelihood. Two common implementations of SAEM, Batch EM (BEM) and online EM (OEM), are parameterized by a “learning rate”, and their efficiency depend strongly on this parameter. We propose an extension to the OEM algorithm, termed Introspective Online Expectation Maximization (IOEM), which removes the need for specifying this parameter by adapting the learning rate to trends in the parameter updates. We show that our algorithm matches the efficiency of the optimal BEM and OEM algorithms in multiple models, and that the efficiency of IOEM can exceed that of BEM/OEM methods with optimal learning rates when the model has many parameters. Finally we use IOEM to fit two models to a financial time series. A Python implementation is available at https://github.com/luntergroup/IOEM.git.

Joint model of accelerated failure time and mechanistic nonlinear model for censored covariates, with application in HIV/AIDS

For a time-to-event outcome with censored time-varying covariates, a joint Cox model with a linear mixed effects model is the standard modeling approach. In some applications such as AIDS studies, mechanistic nonlinear models are available for some covariate process such as viral load during anti-HIV treatments, derived from the underlying data-generation mechanisms and disease progression. Such a mechanistic nonlinear covariate model may provide better-predicted values when the covariates are left censored or mismeasured. When the focus is on the impact of the time-varying covariate process on the survival outcome, an accelerated failure time (AFT) model provides an excellent alternative to the Cox proportional hazard model since an AFT model is formulated to allow the influence of the outcome by the entire covariate process. In this article, we consider a nonlinear mixed effects model for the censored covariates in an AFT model, implemented using a Monte Carlo EM algorithm, under the framework of a joint model for simultaneous inference. We apply the joint model to an HIV/AIDS data to gain insights for assessing the association between viral load and immunological restoration during antiretroviral therapy. Simulation is conducted to compare model performance when the covariate model and the survival model are misspecified.

An approach to model clustered survival data with dependent censoring.

In this study we introduce a likelihood-based method, via the Weibull and piecewise exponential distributions, capable of accommodating the dependence between failure and censoring times. The methodology is developed for the analysis of clustered survival data and it assumes that failure and censoring times are mutually independent conditional on a latent frailty. The dependent censoring mechanism is accounted through the frailty effect and this is accomplished by means of a key parameter accommodating the correlation between failure and censored observations. The full specification of the likelihood in our work simplifies the inference procedures with respect to Huang and Wolfe since it reduces the computation burden of working with the profile likelihood. In addition, the assumptions made for the baseline distributions lead to models with continuous survival functions. In order to carry out inferences, we devise a Monte Carlo EM algorithm. The performance of the proposed models is investigated through a simulation study. Finally, we explore a real application involving patients from the Dialysis Outcomes and Practice Patterns Study observed between 1996 and2015.

Joint mean-covariance random effect model for longitudinal data.

In this paper, we consider the inherent association between mean and covariance in the joint mean-covariance modeling and propose a joint mean-covariance random effect model based on the modified Cholesky decomposition for longitudinal data. Meanwhile, we apply M-H algorithm to simulate the posterior distributions of model parameters. Besides, a computationally efficient Monte Carlo expectation maximization (MCEM) algorithm is developed for carrying out maximum likelihood estimation. Simulation studies show that the model taking into account the inherent association between mean and covariance has smaller standard deviations of the estimators of parameters, which makes the statistical inferences much more reliable. In the real data analysis, the estimation of parameters in the mean and covariance structure is highlyefficient.

Bayesian Estimation for the Extended t-process Regression Models with Independent Errors

The extended t-process regression model is developed to robustly model functional data with outlier functional curves. This paper applies Bayesian estimation to propose an estimation procedure for the model with independent errors. A Monte Carlo EM method is built to estimate parameters involved in the model. Simulation studies and real examples show the proposed method performs well against outliers.

Modelling the aggregate loss for insurance claims with dependence

In this paper, we propose a new model to relax the impractical independence assumption between the counts and the amounts of insurance claims, which is commonly made in the existing literature for mathematical convenience. When considering the dependence between the claim counts and the claim amounts, we treat the number of claims as an explanatory variable in the model for claim sizes. Besides, generalized linear models (GLMs) are employed to fit the claim counts in a given time period. To describe the claim amounts which are repeatedly measured on a group of subjects over time, we adopt generalized linear mixed models (GLMMs) to incorporate the dependence among the related observations on the same subject. In addition, a Monte Carlo Expectation-Maximization (MCEM) algorithm is proposed by using a Metropolis-Hastings algorithm sampling scheme to obtain the maximum likelihood estimates of the parameters for the linear predictor and variance component. Finally, we conduct a simulation to illustrate the feasibility of our proposed model.

A latent class approach for joint modeling of a time-to-event outcome and multiple longitudinal biomarkers subject to limits of detection.

Multiple biomarkers on different biological pathways are often measured over time to investigate the complex mechanism of disease development and progression. Identification of informative subpopulation patterns of longitudinal biomarkers and clinical endpoint may assist in risk stratification and provide insights into new therapeutic targets. Motivated by a multicenter study to assess the inflammatory markers of sepsis in patients with community-acquired pneumonia, we propose a joint latent class analysis of multiple biomarkers and a time-to-event outcome while accounting for censored biomarker measurements due to detection limits. The interrelationship between biomarker trajectories and clinical endpoint is fully captured by a latent class structure, which reveals the subpopulation profiles of biomarkers and clinical outcome. The estimation of joint latent class models becomes more complicated when biomarkers are subject to detection limits. Based on a Metropolis-Hastings method, we develop a Monte Carlo Expectation-Maximization (MCEM) algorithm to estimate model parameters. We demonstrate the satisfactory performance of our MCEM algorithm using simulation studies, and apply our method to the motivating study to examine the heterogeneous patterns of cytokine responses to pneumonia and associated mortality risks.

Fast Iterative Combinatorial Auctions via Bayesian Learning

Iterative combinatorial auctions (CAs) are often used in multibillion dollar domains like spectrum auctions, and speed of convergence is one of the crucial factors behind the choice of a specific design for practical applications. To achieve fast convergence, current CAs require careful tuning of the price update rule to balance convergence speed and allocative efficiency. Brero and Lahaie (2018) recently introduced a Bayesian iterative auction design for settings with singleminded bidders. The Bayesian approach allowed them to incorporate prior knowledge into the price update algorithm, reducing the number of rounds to convergence with minimal parameter tuning. In this paper, we generalize their work to settings with no restrictions on bidder valuations. We introduce a new Bayesian CA design for this general setting which uses Monte Carlo Expectation Maximization to update prices at each round of the auction. We evaluate our approach via simulations on CATS instances. Our results show that our Bayesian CA outperforms even a highly optimized benchmark in terms of clearing percentage and convergence speed.

A multivariate Poisson-log normal mixture model for clustering transcriptome sequencing data

BackgroundHigh-dimensional data of discrete and skewed nature is commonly encountered in high-throughput sequencing studies. Analyzing the network itself or the interplay between genes in this type of data continues to present many challenges. As data visualization techniques become cumbersome for higher dimensions and unconvincing when there is no clear separation between homogeneous subgroups within the data, cluster analysis provides an intuitive alternative. The aim of applying mixture model-based clustering in this context is to discover groups of co-expressed genes, which can shed light on biological functions and pathways of gene products.ResultsA mixture of multivariate Poisson-log normal (MPLN) model is developed for clustering of high-throughput transcriptome sequencing data. Parameter estimation is carried out using a Markov chain Monte Carlo expectation-maximization (MCMC-EM) algorithm, and information criteria are used for model selection.ConclusionsThe mixture of MPLN model is able to fit a wide range of correlation and overdispersion situations, and is suited for modeling multivariate count data from RNA sequencing studies. All scripts used for implementing the method can be found at https://github.com/anjalisilva/MPLNClust.

Gaussian process regression for survival time prediction with genome-wide gene expression.

Predicting the survival time of a cancer patient based on his/her genome-wide gene expression remains a challenging problem. For certain types of cancer, the effects of gene expression on survival are both weak and abundant, so identifying non-zero effects with reasonable accuracy is difficult. As an alternative to methods that use variable selection, we propose a Gaussian process accelerated failure time model to predict survival time using genome-wide or pathway-wide gene expression data. Using a Monte Carlo expectation-maximization algorithm, we jointly impute censored log-survival time and estimate model parameters. We demonstrate the performance of our method and its advantage over existing methods in both simulations and real data analysis. The real data that we analyze were collected from 513 patients with kidney renal clear cell carcinoma and include survival time, demographic/clinical variables, and expression of more than 20 000 genes. In addition to the right-censored survival time, our method can also accommodate left-censored or interval-censored outcomes; and it provides a natural way to combine multiple types of high-dimensional -omics data. An R package implementing our method is available in the Supplementary material available at Biostatistics online.

Genetic parameters for female fertility in Nordic Holstein and Red Cattle dairy breeds

Genetic evaluation of female fertility in Danish, Finnish, and Swedish dairy cows was updated in 2015 to multiple-trait animal model evaluation, where heifer and cow fertility up to third parity are considered as separate traits. A model for conception rate was also developed, which required variance component estimation for Nordic Holstein and Nordic Red Dairy Cattle. We used a multiple-trait multiple-lactation sire model to determine variance components for interval from calving to first insemination, length of service period, and conception rate. Monte Carlo Expectation Maximization REML allowed estimation of all 11 traits simultaneously. Study data were sampled from Swedish Holstein (n = 140,040) and Red Dairy Cattle (n = 101,315) heifers and cows. Conception rate observations are binomial observations with various numbers of failures preceding an observation of success. Using a simulation study, we confirmed that including a service number effect into the conception rate model allowed us to model the change in expectation of successful AI with increasing number of services. Heifers outperformed cows in all fertility traits according to the phenotypic means in the records. Heritabilities for the traits varied from 3 to 7% for interval from calving to first insemination, from 1 to 5% for length of service period, and from 1 to 3% for conception rate. Genetic correlations within traits (i.e., between parities) were favorable, ranging from moderate to high; genetic correlations between heifer and cow traits were lower than between cow traits in different parities. Lowest genetic correlations between traits were for interval from calving to first insemination and conception rate, intermediate for interval from calving to first insemination and length of service period, and highest for length of service period and conception rate. The variance components estimated in this study have been used in Nordic fertility breeding value evaluations since 2016.

Efficient Monte Carlo algorithm for restricted maximum likelihood estimation of genetic parameters.

Monte Carlo (MC) methods have been found useful in estimation of variance parameters for large data and complex models with many variance components (VC), with respect to both computer memory and computing time. A disadvantage has been a fluctuation in round-to-round values of estimates that makes the estimation of convergence challenging. Furthermore, with Newton-type algorithms, the approximate Hessian matrix might have sufficient accuracy, but the inaccuracy in the gradient vector exaggerates the round-to-round fluctuation to intolerable. In this study, the reuse of the same random numbers within each MC sample was used to remove the MC fluctuation. Simulated data with six VC parameters were analysed by four different MC REML methods: expectation-maximization (EM), Newton-Raphson (NR), average information (AI) and Broyden's method (BM). In addition, field data with 96 VC parameters were analysed by MC EM REML. In all the analyses with reused samples, the MC fluctuations disappeared, but the final estimates by the MC REML methods differed from the analytically calculated values more than expected especially when the number of MC samples was small. The difference depended on the random numbers generated, and based on repeated MC AI REML analyses, the VC estimates were on average non-biased. The advantage of reusing MC samples is more apparent in the NR-type algorithms. Smooth convergence opens the possibility to use the fast converging Newton-type algorithms. However, a disadvantage from reusing MC samples is a possible "bias" in the estimates. To attain acceptable accuracy, sufficient number of MC samples need to be generated.

Simplex regression models with measurement error

This paper considers the simplex regression model when there is measurement error in the covariate. We consider a structural approach where the measurement error follows a normal or gamma distribution. We apply a Monte Carlo EM algorithm to estimate the parameters using a pseudo-likelihood function. A simulation study is used to investigate the impact of ignoring the measurement error. Finally, the results are illustrated with a data set.