Search

The advent of wearable and sensor technologies now leads to functional predictors which are intrinsically infinite dimensional. While the existing approaches for functional data and survival outcomes lean on the well-established Cox model, the proportional hazard (PH) assumption might not always be suitable in real-world applications. Motivated by physiological signals encountered in digital medicine, we develop a more general and flexible functional time-transformation model for estimating the conditional survival function with both functional and scalar covariates. A partially functional regression model is used to directly model the survival time on the covariates through an unknown monotone transformation and a known error distribution. We use Bernstein polynomials to model the monotone transformation function and the smooth functional coefficients. A sieve method of maximum likelihood is employed for estimation. Numerical simulations illustrate a satisfactory performance of the proposed method in estimation and inference. We demonstrate the application of the proposed model through two case studies involving wearable data i) Understanding the association between diurnal physical activity pattern and all-cause mortality based on accelerometer data from the National Health and Nutrition Examination Survey (NHANES) 2011-2014 and ii) Modelling Time-to-Hypoglycemia events in a cohort of diabetic patients based on distributional representation of continuous glucose monitoring (CGM) data. The results provide important epidemiological insights into the direct association between survival times and the physiological signals and also exhibit superior predictive performance compared to traditional summary based biomarkers in the CGM study.

Fundamental to our understanding of planetary bulk compositions is the relationship between their masses and radii, two properties that are often not simultaneously known for most exoplanets. However, while many previous studies have modeled the two-dimensional relationship between planetary mass and radii, this approach largely ignores the dependencies on other properties that may have influenced the formation and evolution of the planets. In this work, we extend the existing nonparametric and probabilistic framework of \texttt{MRExo} to jointly model distributions beyond two dimensions. Our updated framework can now simultaneously model up to four observables, while also incorporating asymmetric measurement uncertainties and upper limits in the data. We showcase the potential of this multi-dimensional approach to three science cases: (i) a 4-dimensional joint fit to planetary mass, radius, insolation, and stellar mass, hinting of changes in planetary bulk density across insolation and stellar mass; (ii) a 3-dimensional fit to the California Kepler Survey sample showing how the planet radius valley evolves across different stellar masses; and (iii) a 2-dimensional fit to a sample of Class-II protoplanetary disks in Lupus while incorporating the upper-limits in dust mass measurements. In addition, we employ bootstrap and Monte-Carlo sampling to quantify the impact of the finite sample size as well as measurement uncertainties on the predicted quantities. We update our existing open-source user-friendly \texttt{MRExo} \texttt{Python} package with these changes, which allows users to apply this highly flexible framework to a variety of datasets beyond what we have shown here., Comment: Accepted in ApJ. Updated MRExo package and sample scripts available here: https://github.com/shbhuk/mrexo/tree/v1.0dev. Package will be released on PyPI (pip) along with full documentation upon publication in ApJ

Motivated by the need to model the dependence between regions of interest in functional neuroconnectivity for efficient inference, we propose a new sampling-based Bayesian clustering approach for covariance structures of high-dimensional Gaussian outcomes. The key technique is based on a Dirichlet process that clusters covariance sub-matrices into independent groups of outcomes, thereby naturally inducing sparsity in the whole brain connectivity matrix. A new split-merge algorithm is employed to achieve convergence of the Markov chain that is shown empirically to recover both uniform and Dirichlet partitions with high accuracy. We investigate the empirical performance of the proposed method through extensive simulations. Finally, the proposed approach is used to group regions of interest into functionally independent groups in the Autism Brain Imaging Data Exchange participants with autism spectrum disorder and and co-occurring attention-deficit/hyperactivity disorder., Comment: 23 pages, 8 figures, 2 tables

Modern clinical and epidemiological studies widely employ wearables to record parallel streams of real-time data on human physiology and behavior. With recent advances in distributional data analysis, these high-frequency data are now often treated as distributional observations resulting in novel regression settings. Motivated by these modelling setups, we develop a distributional outcome regression via quantile functions (DORQF) that expands existing literature with three key contributions: i) handling both scalar and distributional predictors, ii) ensuring jointly monotone regression structure without enforcing monotonicity on individual functional regression coefficients, iii) providing statistical inference via asymptotic projection-based joint confidence bands and a statistical test of global significance to quantify uncertainty of the estimated functional regression coefficients. The method is motivated by and applied to Actiheart component of Baltimore Longitudinal Study of Aging that collected one week of minute-level heart rate (HR) and physical activity (PA) data on 781 older adults to gain deeper understanding of age-related changes in daily life heart rate reserve, defined as a distribution of daily HR, while accounting for daily distribution of physical activity, age, gender, and body composition. Intriguingly, the results provide novel insights in epidemiology of daily life heart rate reserve.

Bayesian statistical inference for Generalized Linear Models (GLMs) with parameters lying on a constrained space is of general interest (e.g., in monotonic or convex regression), but often constructing valid prior distributions supported on a subspace spanned by a set of linear inequality constraints can be challenging, especially when some of the constraints might be binding leading to a lower dimensional subspace. For the general case with canonical link, it is shown that a generalized truncated multivariate normal supported on a desired subspace can be used. Moreover, it is shown that such prior distribution facilitates the construction of a general purpose product slice sampling method to obtain (approximate) samples from corresponding posterior distribution, making the inferential method computationally efficient for a wide class of GLMs with an arbitrary set of linear inequality constraints. The proposed product slice sampler is shown to be uniformly ergodic, having a geometric convergence rate under a set of mild regularity conditions satisfied by many popular GLMs (e.g., logistic and Poisson regressions with constrained coefficients). One of the primary advantages of the proposed Bayesian estimation method over classical methods is that uncertainty of parameter estimates is easily quantified by using the samples simulated from the path of the Markov Chain of the slice sampler. Numerical illustrations using simulated data sets are presented to illustrate the superiority of the proposed methods compared to some existing methods in terms of sampling bias and variances. In addition, real case studies are presented using data sets for fertilizer-crop production and estimating the SCRAM rate in nuclear power plants.

The presence of missing values within high-dimensional data is an ubiquitous problem for many applied sciences. A serious limitation of many available data mining and machine learning methods is their inability to handle partially missing values and so an integrated approach that combines imputation and model estimation is vital for down-stream analysis. A computationally fast algorithm, called EMFlow, is introduced that performs imputation in a latent space via an online version of Expectation-Maximization (EM) algorithm by using a normalizing flow (NF) model which maps the data space to a latent space. The proposed EMFlow algorithm is iterative, involving updating the parameters of online EM and NF alternatively. Extensive experimental results for high-dimensional multivariate and image datasets are presented to illustrate the superior performance of the EMFlow compared to a couple of recently available methods in terms of both predictive accuracy and speed of algorithmic convergence. We provide code for all our experiments.

The relationship between mass and radius (M-R relation) is the key for inferring the planetary compositions and thus valuable for the studies of formation and migration models. However, the M-R relation alone is not enough for planetary characterization due to the dependence of it on other confounding variables. This paper provides a non-trivial extension of the M-R relation by including the incident flux as an additional variable. By using Bayesian hierarchical modeling (BHM) that leverages the flexibility of finite mixture models, a probabilistic mass-radius-flux relationship (M-R-F relation) is obtained based on a sample of 319 exoplanets. We find that the flux has nonnegligible impact on the M-R relation, while such impact is strongest for hot-Jupiters. On the population level, the planets with higher level of flux tend to be denser, and high flux could trigger significant mass loss for plants with radii larger than $13R_{\oplus}$. As a result, failing to account for the flux in mass prediction would cause systematic over or under-estimation. With the recent advent of computing power, although a lot of complex statistical models can be fitted using Monte Carlo methods, it has largely remain illusive how to validate these complex models when the data are observed with large measurement errors. We present two novel methods to examine model assumptions, which can be used not only for the models we present in this paper but can also be adapted for other statistical models., Comment: Accepted for publication in MNRAS

Estimating parameters of mixture model has wide applications ranging from classification problems to estimating of complex distributions. Most of the current literature on estimating the parameters of the mixture densities are based on iterative Expectation Maximization (EM) type algorithms which require the use of either taking expectations over the latent label variables or generating samples from the conditional distribution of such latent labels using the Bayes rule. Moreover, when the number of components is unknown, the problem becomes computationally more demanding due to well-known label switching issues \cite{richardson1997bayesian}. In this paper, we propose a robust and quick approach based on change-point methods to determine the number of mixture components that works for almost any location-scale families even when the components are heavy tailed (e.g., Cauchy). We present several numerical illustrations by comparing our method with some of popular methods available in the literature using simulated data and real case studies. The proposed method is shown be as much as 500 times faster than some of the competing methods and are also shown to be more accurate in estimating the mixture distributions by goodness-of-fit tests.

Data augmentation (DA) turns seemingly intractable computational problems into simple ones by augmenting latent missing data. In addition to computational simplicity, it is now well-established that DA equipped with a deterministic transformation can improve the convergence speed of iterative algorithms such as an EM algorithm or Gibbs sampler. In this article, we outline a framework for the transformation-based DA, which we call data transforming augmentation (DTA), allowing augmented data to be a deterministic function of latent and observed data, and unknown parameters. Under this framework, we investigate a novel DTA scheme that turns heteroscedastic models into homoscedastic ones to take advantage of simpler computations typically available in homoscedastic cases. Applying this DTA scheme to fitting linear mixed models, we demonstrate simpler computations and faster convergence rates of resulting iterative algorithms, compared with those under a non-transformation-based DA scheme. We also fit a Beta-Binomial model using the proposed DTA scheme, which enables sampling approximate marginal posterior distributions that are available only under homoscedasticity. An R package, Rdta, is publicly available at CRAN.

Inferring unknown conic sections on the basis of noisy data is a challenging problem with applications in computer vision. A major limitation of the currently available methods for conic sections is that estimation methods rely on the underlying shape of the conics (being known to be ellipse, parabola or hyperbola). A general purpose Bayesian hierarchical model is proposed for conic sections and corresponding estimation method based on noisy data is shown to work even when the specific nature of the conic section is unknown. The model, thus, provides probabilistic detection of the underlying conic section and inference about the associated parameters of the conic section. Through extensive simulation studies where the true conics may not be known, the methodology is demonstrated to have practical and methodological advantages relative to many existing techniques. In addition, the proposed method provides probabilistic measures of uncertainty of the estimated parameters. Furthermore, we observe high fidelity to the true conics even in challenging situations, such as data arising from partial conics in arbitrarily rotated and non-standard form, and where a visual inspection is unable to correctly identify the type of conic section underlying the data.

The proportional hazards (PH), proportional odds (PO) and accelerated failure time (AFT) models have been widely used in different applications of survival analysis. Despite their popularity, these models are not suitable to handle lifetime data with crossing survival curves. In 2005, Yang and Prentice proposed a semiparametric two-sample strategy (YP model), including the PH and PO frameworks as particular cases, to deal with this type of data. Assuming a general regression setting, the present paper proposes an unified approach to fit the YP model by employing Bernstein polynomials to manage the baseline hazard and odds under both the frequentist and Bayesian frameworks. The use of the Bernstein polynomials has some advantages: it allows for uniform approximation of the baseline distribution, it leads to closed-form expressions for all baseline functions, it simplifies the inference procedure, and the presence of a continuous survival function allows a more accurate estimation of the crossing survival time. Extensive simulation studies are carried out to evaluate the behavior of the models. The analysis of a clinical trial data set, related to non-small-cell lung cancer, is also developed as an illustration. Our findings indicate that assuming the usual PH model, ignoring the existing crossing survival feature in the real data, is a serious mistake with implications for those patients in the initial stage of treatment.

In recent years there have been a lot of interest to test for similarity between biological drug products, commonly known as biologics. Biologics are large and complex molecule drugs that are produced by living cells and hence these are sensitive to the environmental changes. In addition, biologics usually induce antibodies which raises the safety and efficacy issues. The manufacturing process is also much more complicated and often costlier than the small-molecule generic drugs. Because of these complexities and inherent variability of the biologics, the testing paradigm of the traditional generic drugs cannot be directly used to test for biosimilarity. Taking into account some of these concerns we propose a functional distance based methodology that takes into consideration the entire time course of the study and is based on a class of flexible semi-parametric models. The empirical results show that the proposed approach is more sensitive than the classical equivalence tests approach which are usually based on arbitrarily chosen time point. Bootstrap based methodologies are also presented for statistical inference.

The well-known Bayes theorem assumes that a posterior distribution is a probability distribution. However, the posterior distribution may no longer be a probability distribution if an improper prior distribution (non-probability measure) such as an unbounded uniform prior is used. Improper priors are often used in the astronomical literature to reflect a lack of prior knowledge, but checking whether the resulting posterior is a probability distribution is sometimes neglected. It turns out that 23 articles out of 75 articles (30.7%) published online in two renowned astronomy journals (ApJ and MNRAS) between Jan 1, 2017 and Oct 15, 2017 make use of Bayesian analyses without rigorously establishing posterior propriety. A disturbing aspect is that a Gibbs-type Markov chain Monte Carlo (MCMC) method can produce a seemingly reasonable posterior sample even when the posterior is not a probability distribution (Hobert and Casella, 1996). In such cases, researchers may erroneously make probabilistic inferences without noticing that the MCMC sample is from a non-existing probability distribution. We review why checking posterior propriety is fundamental in Bayesian analyses, and discuss how to set up scientifically motivated proper priors.

A Gaussian measurement error assumption, i.e., an assumption that the data are observed up to Gaussian noise, can bias any parameter estimation in the presence of outliers. A heavy tailed error assumption based on Student's t distribution helps reduce the bias. However, it may be less efficient in estimating parameters if the heavy tailed assumption is uniformly applied to all of the data when most of them are normally observed. We propose a mixture error assumption that selectively converts Gaussian errors into Student's t errors according to latent outlier indicators, leveraging the best of the Gaussian and Student's t errors; a parameter estimation can be not only robust but also accurate. Using simulated hospital profiling data and astronomical time series of brightness data, we demonstrate the potential for the proposed mixture error assumption to estimate parameters accurately in the presence of outliers.

Measurements of radial velocity variations from the spectroscopic monitoring of stars and their companions are essential for a broad swath of astrophysics, providing access to the fundamental physical properties that dictate all phases of stellar evolution and facilitating the quantitative study of planetary systems. The conversion of those measurements into both constraints on the orbital architecture and individual component spectra can be a serious challenge, however, especially for extreme flux ratio systems and observations with relatively low sensitivity. Gaussian processes define sampling distributions of flexible, continuous functions that are well-motivated for modeling stellar spectra, enabling proficient search for companion lines in time-series spectra. We introduce a new technique for spectral disentangling, where the posterior distributions of the orbital parameters and intrinsic, rest-frame stellar spectra are explored simultaneously without needing to invoke cross-correlation templates. To demonstrate its potential, this technique is deployed on red-optical time-series spectra of the mid-M dwarf binary LP661-13. We report orbital parameters with improved precision compared to traditional radial velocity analysis and successfully reconstruct the primary and secondary spectra. We discuss potential applications for other stellar and exoplanet radial velocity techniques and extensions to time-variable spectra. The code used in this analysis is freely available as an open source Python package., Comment: submitted to ApJ

In recent years, the Hamiltonian Monte Carlo (HMC) algorithm has been found to work more efficiently compared to other popular Markov Chain Monte Carlo (MCMC) methods (such as random walk Metropolis-Hastings) in generating samples from a posterior distribution. A general framework for HMC based on the use of graphical processing units (GPUs) is shown to greatly reduce the computing time needed for Bayesian inference. The most expensive computational tasks in HMC are the evaluation of the posterior kernel and computing its gradient with respect to the parameters of interest. One of primary goals of this article to show that by expressing each of these tasks in terms of simple matrix or element-wise operations and maintaining persistent objects in GPU memory, the computational time can be drastically reduced. By using GPU objects to perform the entire HMC simulation, most of the latency penalties associated with transferring data from main to GPU memory can be avoided. Thus, the proposed computational framework is conceptually very simple, but also is general enough to be applied to most problems that use HMC sampling. For clarity of exposition, the effectiveness of the proposed approach is demonstrated in the high-dimensional setting on a standard statistical model - multinomial regression. Using GPUs, analyses of data sets that were previously intractable for fully Bayesian approaches due to the prohibitively high computational cost are now feasible.