Published papers

Unmeasured confounders pose a major challenge in accurately estimating causal effects in observational studies. To address this issue when estimating hazard ratios (HRs) using Cox proportional hazards models, several methods, including instrumental variables (IVs) approaches, have been proposed. However, these methods often face limitations, such as weak IV problems and restrictive assumptions regarding unmeasured confounder distributions. In this study, we introduce a novel nonparametric Bayesian procedure that provides accurate HR estimates while addressing these limitations. A key assumption of our approach is that unmeasured confounders exhibit a cluster structure. Under this assumption, we integrate two remarkable Bayesian techniques, the Dirichlet process mixture (DPM) and general Bayes (GB), to simultaneously (1) detect latent clusters based on the likelihood of exposure and outcome variables and (2) estimate HRs using the likelihood constructed within each cluster. Notably, leveraging DPM, our procedure eliminates the need for IVs by identifying unmeasured confounders under an alternative condition. Additionally, GB techniques remove the need for explicit modeling of the baseline hazard function, distinguishing our procedure from traditional Bayesian approaches. Simulation experiments demonstrate that the proposed Bayesian procedure outperforms existing methods in some performance metrics. Moreover, it achieves statistical efficiency comparable to the efficient estimator while accurately identifying cluster structures. These features highlight its ability to overcome challenges associated with traditional IV approaches for time-to-event data.

In this article, we introduce a new and efficient data augmentation approach to the posterior inference of the models with shape parameters when the reciprocal gamma function appears in full conditional densities. Our approach is to approximate full conditional densities of shape parameters by using Gauss’s multiplication formula and Stirling’s formula for the gamma function, where the approximation error can be made arbitrarily small. We use the techniques to construct efficient Gibbs and Metropolis–Hastings algorithms for a variety of models that involve the gamma distribution, Student’s t-distribution, the Dirichlet distribution, the negative binomial distribution, and the Wishart distribution. The proposed sampling method is numerically demonstrated through simulation studies. Supplementary materials for this article are available online.

Global-local shrinkage priors have been recognized as a useful class of priors that can strongly shrink small signals toward prior means while keeping large signals unshrunk. Although such priors have been extensively discussed under Gaussian responses, in practice, we often encounter count responses. Previous contributions on global-local shrinkage priors cannot be readily applied to count data. In this paper, we discuss global-local shrinkage priors for analyzing a sequence of counts. We provide sufficient conditions under which the posterior mean is unshrunk for very large signals, known as the tail robustness property. Then, we propose tractable priors to satisfy those conditions approximately or exactly and develop a custom posterior computation algorithm for Bayesian inference without tuning parameters. We demonstrate the proposed methods through simulation studies and an application to a real dataset.

Although robust divergence, such as density power divergence and γ-divergence, is helpful for robust statistical inference in the presence of outliers, the tuning parameter that controls the degree of robustness is chosen in a rule-of-thumb, which may lead to an inefficient inference. We here propose a selection criterion based on an asymptotic approximation of the Hyvarinen score applied to an unnormalized model defined by robust divergence. The proposed selection criterion only requires first and second-order partial derivatives of an assumed density function with respect to observations, which can be easily computed regardless of the number of parameters. We demonstrate the usefulness of the proposed method via numerical studies using normal distributions and regularized linear regression.

Various multivariate extensions to the well-known Fay-Herriot model have been proposed in the small area estimation literature. Such multivariate models are quite effective in combining information through correlations among small area survey estimates of related variables or historical survey estimates of the same variable or both. Though the literature on small area estimation is already very rich, construction of second-order efficient confidence intervals from multivariate models has received little attention. In this article, we develop a parametric bootstrap method for constructing a second-order efficient confidence interval for a general linear combination of small area means using the multivariate Fay-Herriot normal model. The proposed parametric bootstrap method replaces difficult and tedious analytical derivations by the power of efficient algorithm and high speed computer. Moreover, the proposed method is more versatile than the analytical method because the parametric bootstrap method can be easily applied to any method of model parameter estimation and any specific structure of the variance-covariance matrix of the multivariate Fay-Herriot model avoiding all the cumbersome and time-consuming calculations required in the analytical method. We apply our proposed methodology in constructing confidence intervals for the median income of four-person families for the fifty states and the District of Columbia in the United States. Our data analysis demonstrates that the proposed parametric bootstrap method, applied to both multivariate and univariate Fay-Herriot models, generally provides much shorter confidence intervals compared to the corresponding traditional direct method. Moreover, the confidence intervals obtained from the multivariate model are generally shorter than the corresponding intervals from the univariate model indicating the potential advantage of exploiting correlations of median income of four-person families with median incomes of three- and five-person families.

Robust Bayesian linear regression is a classical but essential statistical tool. Although novel robustness properties of posterior distributions have been proved recently under a certain class of error distributions, their sufficient conditions are restrictive and exclude several important situations. In this work, we revisit a classical two-component mixture model for response variables, also known as contamination model, where one component is a light-tailed regression model and the other component is heavy-tailed. The latter component is independent of the regression parameters, which is crucial in proving the posterior robustness. We obtain new sufficient conditions for posterior (non-)robustness and reveal non-trivial robustness results by using those conditions. In particular, we find that even the Student-t error distribution can achieve the posterior robustness in our framework. A numerical study is performed to check the Kullback–Leibler divergence between the posterior distribution based on full data and that based on data obtained by removing outliers.

This paper proposes a novel methodology called the mixture of Bayesian predictive syntheses (MBPS) for multiple time series count data for the chal-lenging task of predicting the numbers of COVID-19 inpatients and isolated cases in Japan and Korea at the subnational level. MBPS combines a set of predictive models and partitions the multiple time series into clusters based on their contribution to predicting the outcome. In this way MBPS lever-ages the shared information within each cluster and is suitable for predicting COVID-19 inpatients since the data exhibit similar dynamics over multiple areas. Also, MBPS avoids using a multivariate count model, which is gener-ally cumbersome to develop and implement. Our Japanese and Korean data analyses demonstrate that the proposed MBPS methodology has improved predictive accuracy and uncertainty quantification.

We introduce a new robust Bayesian change-point analysis in the presence of outliers. We employ an idea of general posterior based on density power divergence combined with horseshoe prior for differences of underlying signals. A posterior computation algorithm is proposed using Markov chain Monte Carlo. The proposed method is demonstrated through simulation and real data analysis.

Count data with zero inflation and large outliers are ubiquitous in many scientific applications. However, posterior analysis under a standard statistical model, such as Poisson or negative binomial distribution, is sensitive to such contamination. This study introduces a novel framework for Bayesian modeling of counts that is robust to both zero inflation and large outliers. In doing so, we introduce rescaled beta distribution and adopt it to absorb undesirable effects from zero and outlying counts. The proposed approach has two appealing features: the efficiency of the posterior computation via a custom Gibbs sampling algorithm and a theoretically guaranteed posterior robustness, where extreme outliers are automatically removed from the posterior distribution. We demonstrate the usefulness of the proposed method by applying it to trend filtering and spatial modeling using predictive Gaussian processes. Supplementary materials for this article are available online, including a standardized description of the materials available for reproducing the work.

Mixture modeling, which considers the potential heterogeneity in data, is widely adopted for classification and clustering problems. Mixture models can be estimated using the Expectation-Maximization algorithm, which works with the complete estimating equations conditioned by the latent membership variables of the cluster assignment based on the hierarchical expression of mixture models. However, when the mixture components have light tails such as a normal distribution, the mixture model can be sensitive to outliers. This study proposes a method of weighted complete estimating equations (WCE) for the robust fitting of mixture models. Our WCE introduces weights to complete estimating equations such that the weights can automatically downweight the outliers. The weights are constructed similarly to the density power divergence for mixture models, but in our WCE, they depend only on the component distributions and not on the whole mixture. A novel expectation-estimating-equation (EEE) algorithm is also developed to solve the WCE. For illustrative purposes, a multivariate Gaussian mixture, a mixture of experts, and a multivariate skew normal mixture are considered, and how our EEE algorithm can be implemented for these specific models is described. The numerical performance of the proposed robust estimation method was examined using simulated and real datasets.

Meta-analysis is an essential tool to comprehensively synthesize and quantitatively evaluate results of multiple clinical studies in evidence-based medicine. In many meta-analyses, the characteristics of some studies might markedly differ from those of the others, and these outlying studies can generate biases and potentially yield misleading results. In this article, we provide effective robust statistical inference methods using generalized likelihoods based on the density power divergence. The robust inference methods are designed to adjust the influences of outliers through the use of modified estimating equations based on a robust criterion, even when multiple and serious influential outliers are present. We provide the robust estimators, statistical tests, and confidence intervals via the generalized likelihoods for the fixed-effect and random-effects models of meta-analysis. We also assess the contribution rates of individual studies to the robust overall estimators that indicate how the influences of outlying studies are adjusted. Through simulations and applications to two recently published systematic reviews, we demonstrate that the overall conclusions and interpretations of meta-analyses can be markedly changed if the robust inference methods are applied and that only the conventional inference methods might produce misleading evidence. These methods would be recommended to be used at least as a sensitivity analysis method in the practice of meta-analysis. We have also developed an R package, robustmeta, that implements the robust inference methods.

Linear mixed models (LMMs) are a popular class of methods for analyzing longitudinal and clustered data. However, such models can be sensitive to outliers, and this can lead to biased inference on model parameters and inaccurate prediction of random effects if the data are contaminated. We propose a new approach to robust estimation and inference for LMMs using a hierarchical gamma-divergence, which offers an automated, data-driven approach to downweight the effects of outliers occurring in both the error and the random effects, using normalized powered density weights. For estimation and inference, we develop a computationally scalable minorization-maximization algorithm for the resulting objective function, along with a clustered bootstrap method for uncertainty quantification and a Hyvarinen score criterion for selecting a tuning parameter controlling the degree of robustness. Under suitable regularity conditions, we show the resulting robust estimates can be asymptotically controlled even under a heavy level of (covariate-dependent) contamination. Simulation studies demonstrate hierarchical gamma-divergence consistently outperforms several currently available methods for robustifying LMMs. We also illustrate the proposed method using data from a multi-center AIDS cohort study. Supplementary materials for this article are available online.

This paper proposes a flexible Bayesian approach to multiple imputation based on conditional Gaussian mixture models. We introduce shrinkage priors for covariate-dependent mixing proportions, which allow the number of mixture components used for imputation to be selected automatically. We develop an efficient Markov chain Monte Carlo algorithm for posterior computation and imputation. The proposed method extends naturally to mixed data containing both continuous and discrete variables (e.g., binary and count outcomes). We also propose an approximate Bayesian inference procedure for parameters defined through loss functions, using the posterior predictive distribution of missing observations and extending bootstrap-based Bayesian inference for complete data. Numerical studies using simulated and real data demonstrate the proposed method.

Random partition distribution is a crucial tool for model-based clustering. This study advances the field of random partition in the context of functional spatial data, focusing on the challenges posed by hourly population data across various regions and dates. We propose an extension of the generalized Dirichlet process, named the similarity-based generalized Dirichlet process (SGDP)-type distribution, to address the limitations of simple random partition distributions (e.g. those induced by the Dirichlet process), such as an overabundance of clusters. This model prevents excess cluster production and incorporates pairwise similarity information to ensure accurate and meaningful clustering. The theoretical properties of the SGDP-type distribution are studied. Then, SGDP-type random partition is applied to a real-world dataset of hourly population flow in 500 m meshes in the central part of Tokyo. In this empirical context, our method excels at detecting meaningful patterns in the data while accounting for spatial nuances. The results underscore the adaptability and utility of the method, showcasing its prowess in revealing intricate spatiotemporal dynamics. The proposed random partition will significantly contribute to urban planning, transportation, and policy-making and will be a helpful tool for understanding population dynamics and their implications.

In various applications, we deal with high-dimensional positive-valued data that often exhibits sparsity. This paper develops a new class of continuous global-local shrinkage priors tailored to analyzing gamma-distributed observations where most of the underlying means are concentrated around a certain value. Unlike existing shrinkage priors, our new prior is a shape-scale mixture of inversegamma distributions, which has a desirable interpretation of the form of posterior mean and admits flexible shrinkage. We show that the proposed prior has two desirable theoretical properties; Kullback-Leibler super-efficiency under sparsity and robust shrinkage rules for large observations.We propose an efficient sampling algorithm for posterior inference. The performance of the proposed method is illustrated through simulation and two real data examples, the average length of hospital stay for COVID-19 in South Korea and adaptive variance estimation of gene expression data.

Spatial regression and geographically weighted regression models have been widely adopted to capture the effects of auxiliary information on a response variable of interest over a region. In contrast, relationships between response and auxiliary variables are expected to exhibit complex spatial patterns in many applications. This paper proposes a new approach for spatial regression, called spatially clustered regression, to estimate possibly clustered spatial patterns of the relationships. We combine K-means-based clustering formulation and penalty function motivated from a spatial process known as Potts model for encouraging similar clustering in neighboring locations. We provide a simple iterative algorithm to fit the proposed method, scalable for large spatial datasets. Through simulation studies, the proposed method demonstrates its superior performance to existing methods even under the true structure does not admit spatial clustering. Finally, the proposed method is applied to crime event data in Tokyo and produces interpretable results for spatial patterns. The R code is available at https://github.com/sshonosuke/SCR.

The proliferation of mobile devices has led to the collection of large amounts of population data. This situation has prompted the need to utilize this rich, multidimensional data in practical applications. In response to this trend, we have integrated functional data analysis (FDA) and factor analysis to address the challenge of predicting hourly population changes across various districts in Tokyo. Specifically, by assuming a Gaussian process, we avoided the large covariance matrix parameters of the multivariate normal distribution. In addition, the data were both time and spatially dependent between districts. To capture various characteristics, a Bayesian factor model was introduced, which modeled the time series of a small number of common factors and expressed the spatial structure through factor loading matrices. Furthermore, the factor loading matrices were made identifiable and sparse to ensure the interpretability of the model. We also proposed a Bayesian shrinkage method as a systematic approach for factor selection. Through numerical experiments and data analysis, we investigated the predictive accuracy and interpretability of our proposed method. We concluded that the flexibility of the method allows for the incorporation of additional time series features, thereby improving its accuracy.

This study proposes a method for aggregating/synthesizing global and local sub-models for fast and flexible spatial regression modeling. Eigenvector spatial filtering (ESF) was used to model spatially varying coefficients and spatial dependence in the residuals by sub-model, while the generalized product-of-experts method was used to aggregate these sub-models. The major advantages of the proposed method are as follows: (i) it is highly scalable for large samples in terms of accuracy and computational efficiency; (ii) it is easily implemented by estimating sub-models independently first and aggregating/averaging them thereafter; and (iii) likelihood-based inference is available because the marginal likelihood is available in closed-form. The accuracy and computational efficiency of the proposed method are confirmed using Monte Carlo simulation experiments. This method was then applied to residential land price analysis in Japan. The results demonstrate the usefulness of this method for improving the interpretability of spatially varying coefficients. The proposed method is implemented in an R package spmoran.

Shrinkage priors are a popular Bayesian paradigm to handle sparsity in high-dimensional regression. Still limited, however, is a flexible class of shrinkage priors to handle grouped sparsity, where covariates exhibit some natural grouping structure. This paper proposes a novel extension of the R2-induced Dirichlet Decomposition (R2D2) prior to accommodate grouped variable selection in linear regression models. The proposed method, called the gR2D2 prior, employs a Dirichlet prior distribution on the coefficient of determination for each group, allowing for a flexible and adaptive shrinkage that operates at both group and individual variable levels. We present several theoretical properties of this proposed prior distribution while also developing a Markov Chain Monte Carlo algorithm. Through simulation studies and real-data analysis, we demonstrate that our method outperforms traditional shrinkage priors in terms of both estimation accuracy, inference, and prediction.

Despite increasing accessibility to function data, effective methods for flexibly estimating underlying functional trend are still scarce. We thereby develop a functional version of trend filtering for estimating trend of functional data indexed by time or on general graph by extending the conventional trend filtering, a powerful nonparametric trend estimation technique for scalar data. We formulate the new trend filtering by introducing penalty terms based on (Figure presented.) -norm of the differences of adjacent trend functions. We develop an efficient iteration algorithm for optimizing the objective function obtained by orthonormal basis expansion. Furthermore, we introduce additional penalty terms to eliminate redundant basis functions, which leads to automatic adaptation of the number of basis functions. The tuning parameter in the proposed method is selected via cross validation. We demonstrate the proposed method is locally adaptive and can identify change points through simulation studies and applications to real-world datasets.