Published papers

Factor analysis has been extensively used to reveal the dependence structures among multivariate variables, offering valuable insight in various fields. However, it cannot incorporate the spatial heterogeneity that is typically present in spatial data. To address this issue, we introduce an effective method specifically designed to discover the potential dependence structures in multivariate spatial data. Our approach assumes that spatial locations can be approximately divided into a finite number of clusters, with locations within the same cluster sharing similar dependence structures. By leveraging an iterative algorithm that combines spatial clustering with factor analysis, we simultaneously detect spatial clusters and estimate a unique factor model for each cluster. The proposed method is evaluated through comprehensive simulation studies, demonstrating its flexibility. In addition, we apply the proposed method to a dataset of railway station attributes in the Tokyo metropolitan area, highlighting its practical applicability and effectiveness in uncovering complex spatial dependencies.

ANOVA-based estimators of variance components for nested-error regression models are always constructed based on moment equations through residual variance. We consider moment equations associated with residual covariance and construct improved ANOVA-based estimators. The proposed estimators have closed-form analytic expressions, which enables easy computation. Moreover, they are shown to be consistent, asymptotically unbiased, and robust to the choice of distribution of the random effects. These estimators have comparable and often better performances than many traditional estimators of variance components like the Prasad-Rao, maximum likelihood, and the restricted maximum likelihood estimators for almost all kinds of sample allocations. Their improved performances are demonstrated analytically as well as through detailed simulation studies and applications to real data sets.

We consider modeling and prediction of Capelin distribution in the Barents Sea based on zero-inflated count observation data that vary continuously over a specified survey region. The model is a mixture of two components; a one-point distribution at the origin and a Poisson distribution with spatio-temporal intensity, where both intensity and mixing proportions are modeled by some auxiliary variables and unobserved spatio-temporal effects. The spatio-temporal effects are modeled by a dynamic linear model combined with the predictive Gaussian process. We develop an efficient posterior computational algorithm for the model using a data augmentation strategy. The performance of the proposed model is demonstrated through simulation studies, and an application to the number of Capelin caught in the Barents Sea from 2014 to 2019.

The development of molecular diagnostic tools to achieve individualized medicine requires identifying predictive biomarkers associated with subgroups of individuals who might receive beneficial or harmful effects from different available treatments. However, due to the large number of candidate biomarkers in the large-scale genetic and molecular studies, and complex relationships among clinical outcome, biomarkers, and treatments, the ordinary statistical tests for the interactions between treatments and covariates have difficulties from their limited statistical powers. In this paper, we propose an efficient method for detecting predictive biomarkers. We employ weighted loss functions of Chen et al. to directly estimate individual treatment scores and propose synthetic posterior inference for effect sizes of biomarkers. We develop an empirical Bayes approach, namely, we estimate unknown hyperparameters in the prior distribution based on data. We then provide efficient screening methods for the candidate biomarkers via optimal discovery procedure with adequate control of false discovery rate. The proposed method is demonstrated in simulation studies and an application to a breast cancer clinical study in which the proposed method was shown to detect the much larger numbers of significant biomarkers than existing standard methods.

In genetic association studies, rare variants with extremely low allele frequencies play a crucial role in complex traits. Therefore, set-based testing methods that jointly assess the effects of groups of single nucleotide polymorphisms (SNPs) were developed to increase the powers of the association tests. However, these powers are still insufficient, and precise estimations of the effect sizes of individual SNPs are largely impossible. In this article, we provide an efficient set-based statistical inference framework that addresses both of these important issues simultaneously using an empirical Bayes method with semiparametric multilevel mixture modeling. We propose to utilize the hierarchical model that incorporates variations in set-specific effects and to apply the optimal discovery procedure (ODP) that achieves the largest overall power in multiple significance testing. In addition, we provide an optimal “set-based” estimator of the empirical distribution of effect sizes. The efficiency of the proposed methods is demonstrated through application to a genome-wide association study of coronary artery disease and through simulation studies. The results demonstrated numerous rare variants with large effect sizes for coronary artery disease, and the number of significant sets detected by the ODP was much greater than those identified by existing methods.

This study proposes a novel approach to ensemble prediction, called “covariate-dependent stacking” (CDST). Unlike traditional stacking and model averaging methods, CDST allows model weights to vary flexibly as a function of covariates, thereby enhancing predictive performance in complex scenarios. We formulate the covariate-dependent weights through combinations of basis functions and estimate them by optimizing cross-validation. To analyze the theoretical properties, we establish an oracle inequality regarding the expected loss to be minimized for estimating model weights. Through comprehensive simulation studies and an application to large-scale land price prediction, we demonstrate that the CDST consistently outperforms conventional model averaging methods, particularly on datasets where base models fail to capture the underlying complexity. Our findings suggest that the CDST is especially valuable for, but not limited to, spatio-temporal prediction problems, offering a powerful tool for researchers and practitioners across a wide spectrum of data analysis fields.

Quantiles are useful characteristics of random variables that can provide substantial information on distributions compared with commonly used summary statistics such as means. In this study, we propose a Bayesian quantile trend filtering method to estimate the non-stationary trend of quantiles. We introduce general shrinkage priors to induce locally adaptive Bayesian inference on trends and mixture representation of the asymmetric Laplace likelihood. To quickly compute the posterior distribution, we develop calibrated mean-field variational approximations to guarantee that the frequentist coverage of credible intervals obtained from the approximated posterior is a specified nominal level. Simulation and empirical studies show that the proposed algorithm is computationally much more efficient than the Gibbs sampler and tends to provide stable inference results, especially for high/low quantiles.

Count data frequently arises in biomedical applications, such as the length of hospital stay. However, their discrete nature poses significant challenges for appropriately modeling conditional quantiles, which are crucial for understanding heterogeneous effects and variability in outcomes. To solve the practical difficulty, we propose a novel general Bayesian framework for quantile regression tailored to count data. We seek the regression parameter on the conditional quantile by minimizing the expected loss with respect to the distribution of the conditional quantile of the latent continuous variable associated with the observed count response variable. By modeling the unknown conditional distribution through a Bayesian nonparametric kernel mixture for the joint distribution of the count response and covariates, we obtain the posterior distribution of the regression parameter via a simple optimization. We numerically demonstrate that the proposed method improves bias and estimation accuracy of the existing crude approaches to count quantile regression. Furthermore, we analyze the length of hospital stay for acute myocardial infarction and demonstrate that the proposed method gives more interpretable and flexible results than the existing ones.

This study attempts to predict and forecast the future heterogeneous increase in the vacant house ratio among prefectures in Japan using spatial panel models with unobserved dynamic spatiotemporal effects. The study formulated models with autoregressive and random-walk spatiotemporal effects, referring to the dynamic spatiotemporal effects (DSE) models. We estimate the model parameters and latent spatiotemporal effects via Markov Chain Monte Carlo algorithm. Simulation studies demonstrated the superior performance of the DSE model in terms of future prediction when spatial and temporal correlation exists. The model is then applied to the prefecturewise ratio of vacant houses to the non-rental housing stock in Japan, and the results imply existence of spatiotemporal correlations that cannot be captured by explanatory variables. Furthermore, it is revealed that the DSE models can provide better forecasting than the existing spatial panel models.

Recently, Jewson and Rossell proposed a new approach for kernel density estimation using an exponentiated form of kernel density estimators. The density estimator contained two hyperparameters that flexibly controls the smoothness of the resulting density. We tune them in a data-driven manner by minimizing an objective function based on the Hyvärinen score to avoid the optimization involving the intractable normalizing constant caused by the exponentiation. We show the asymptotic properties of the proposed estimator and emphasize the importance of including the two hyperparameters for flexible density estimation. Our simulation studies and application to income data show that the proposed density estimator is promising when the underlying density is multi-modal or when observations contain outliers.

This paper introduces a general framework for estimating variance components in the linear mixed models via general unbiased estimating equations, which include some well-used estimators such as the restricted maximum likelihood estimator. We derive the asymptotic covariance matrices and second-order biases under general estimating equations without assuming the normality of the underlying distributions and identify a class of second-order unbiased estimators of variance components. It is also shown that the asymptotic covariance matrices and second-order biases do not depend on whether the regression coefficients are estimated by the generalized or ordinary least squares methods. We carry out numerical studies to check the performance of the proposed methods based on typical linear mixed models.

Sampling from matrix generalized inverse Gaussian (MGIG) distributions is required in Markov chain Monte Carlo (MCMC) algorithms for a variety of statistical models. However, an efficient sampling scheme for the MGIG distributions has not been fully developed. We here propose a novel blocked Gibbs sampler for the MGIG distributions based on the Cholesky decomposition. We show that the full conditionals of the entries of the diagonal and unit lower-triangular matrices are univariate generalized inverse Gaussian and multivariate normal distributions, respectively. Several variants of the Metropolis-Hastings algorithm can also be considered for this problem, but we mathematically prove that the average acceptance rates become extremely low in particular scenarios. We demonstrate the computational efficiency of the proposed Gibbs sampler through simulation studies and data analysis. Supplementary materials for this article are available online.

Generalized estimating equation (GEE) is widely adopted for regression modeling for longitudinal data, taking account of potential correlations within the same subjects. Although the standard GEE assumes common regression coefficients among all the subjects, such an assumption may not be realistic when there is potential heterogeneity in regression coefficients among subjects. In this paper, we develop a flexible and interpretable approach, called grouped GEE analysis, to modeling longitudinal data with allowing heterogeneity in regression coefficients. The proposed method assumes that the subjects are divided into a finite number of groups and subjects within the same group share the same regression coefficient. We provide a simple algorithm for grouping subjects and estimating the regression coefficients simultaneously, and show the asymptotic properties of the proposed estimator. The number of groups can be determined by the cross validation with averaging method. We demonstrate the proposed method through simulation studies and an application to a real data set.

Clustered data are ubiquitous in a variety of scientific fields. In this article, we propose a flexible and interpretable modeling approach, called grouped heterogeneous mixture modeling, for clustered data, which models cluster-wise conditional distributions by mixtures of latent conditional distributions common to all the clusters. In the model, we assume that clusters are divided into a finite number of groups and mixing proportions are the same within the same group. We provide a simple generalized EM algorithm for computing the maximum likelihood estimator, and an information criterion to select the numbers of groups and latent distributions. We also propose structured grouping strategies by introducing penalties on grouping parameters in the likelihood function. Under the settings where both the number of clusters and cluster sizes tend to infinity, we present asymptotic properties of the maximum likelihood estimator and the information criterion. We demonstrate the proposed method through simulation studies and an application to crime risk modeling in Tokyo.

Meta-analyses of diagnostic test accuracy (DTA) studies have been gathering attention in research in clinical epidemiology and health technology development, and bivariate random-effects model is becoming a standard tool. However, standard inference methods usually underestimate statistical errors and possibly provide highly overconfident results under realistic situations since they ignore the variability in the estimation of variance parameters. To overcome the difficulty, a new improved inference method, namely, an accurate confidence region for the meta-analysis of DTA, by asymptotically expanding the coverage probability of the standard confidence region. The advantage of the proposed confidence region is that it holds a relatively simple expression and does not require any repeated calculations such as Bootstrap or Monte Carlo methods to compute the region, thereby the proposed method can be easily carried out in practical applications. The effectiveness of the proposed method is demonstrated through simulation studies and an application to meta-analysis of screening test accuracy for alcohol problems.

Network meta-analysis has played an important role in evidence-based medicine for assessing the comparative effectiveness of multiple available treatments. The prediction interval has been one of the standard outputs in recent network meta-analysis as an effective measure that enables simultaneous assessment of uncertainties in treatment effects and heterogeneity among studies. To construct the prediction interval, a large-sample approximating method based on the t-distribution has generally been applied in practice; however, recent studies have shown that similar t-approximation methods for conventional pairwise meta-analyses can substantially underestimate the uncertainty under realistic situations. In this article, we performed simulation studies to assess the validity of the current standard method for network meta-analysis, and we show that its validity can also be violated under realistic situations. To address the invalidity issue, we developed two new methods to construct more accurate prediction intervals through bootstrap and Kenward–Roger-type adjustment. In simulation experiments, the two proposed methods exhibited better coverage performance and generally provided wider prediction intervals than the ordinary t-approximation. We also developed an R package, PINMA (https://cran.r-project.org/web/packages/PINMA/), to perform the proposed methods using simple commands. We illustrate the effectiveness of the proposed methods through applications to two real network meta-analyses.

This study investigates the relationship between access to domestic and international communication and economic development. It does so by constructing two indices of linguistic distance, domestic and international, capturing language acquisition costs, which are higher when acquiring linguistically more distant languages. The domestic linguistic distance index captures the constraints of communication among speakers of different mother tongues within a country, while the international linguistic distance index captures the constraints of global communication via English. This study’s results reveal a negative association between domestic linguistic distance and GDP per capita, whereas international linguistic distance has no significant association. Moreover, by investigating the mechanisms of the negative association of the domestic linguistic distance, we find that communication difficulty among different language groups hinders economic development through a channel of less employment requiring communication-based operations. Furthermore, we determine that the negative association of the domestic linguistic distance may be mainly driven by relatively poor countries such as many in Africa.

Isotonic regression or monotone function estimation is a problem of estimating function values under monotonicity constraints, which appears naturally in many scientific fields. This paper proposes a new Bayesian method with global–local shrinkage priors for estimating monotone function values. Specifically, we introduce half shrinkage priors for positive valued random variables and assign them for the first-order differences of function values. We also develop fast and simple Gibbs sampling algorithms for full posterior analysis. By incorporating advanced shrinkage priors, the proposed method is adaptive to local abrupt changes or jumps in target functions. We show this adaptive property theoretically by proving that the posterior mean estimators are robust to large differences and that asymptotic risk for unchanged points can be improved. Finally, we demonstrate the proposed methods through simulations and applications to a real data set.

Spatial trend estimation under potential heterogeneity is an important problem to extract spatial characteristics and hazards such as criminal activity. By focusing on quantiles, which provide substantial information on distributions compared with commonly used summary statistics such as means, it is often useful to estimate not only the average trend but also the high (low) risk trend additionally. In this paper, we propose a Bayesian quantile trend filtering method to estimate the non-stationary trend of quantiles on graphs and apply it to crime data in Tokyo between 2013 and 2017. By modeling multiple observation cases, we can estimate the potential heterogeneity of spatial crime trends over multiple years in the application. To induce locally adaptive Bayesian inference on trends, we introduce general shrinkage priors for graph differences. Introducing so-called shadow priors with multivariate distribution for local scale parameters and mixture representation of the asymmetric Laplace distribution, we provide a simple Gibbs sampling algorithm to generate posterior samples. The numerical performance of the proposed method is demonstrated through simulation studies.

Linear regression that employs the assumption of normality for the error distribution may lead to an undesirable posterior inference of regression coefficients due to potential outliers. A finite mixture of two components, one with thin and one with heavy tails, is considered as the error distribution in this study. For the heavily-tailed component, the novel class of distributions is introduced; their densities are log-regularly varying and have heavier tails than the Cauchy distribution. Yet, they are expressed as a scale mixture of normals which enables the efficient posterior inference when using a Gibbs sampler. The robustness of the posterior distributions is proved under the proposed models using a minimal set of assumptions, which justifies the use of shrinkage priors with unbounded densities for the coefficient vector in the presence of outliers. An extensive comparison with the existing methods via simulation study shows the improved performance of the proposed model in point and interval estimation, as well as its computational efficiency. Further, the posterior robustness of the proposed method is confirmed in an empirical study with shrinkage priors for regression coefficients.