George Karabatsos

Many statistical problems can be addressed by applying a multiple testing procedure (MTP) that controls either the Family-Wise Error Rate (FWER) or False Discovery Rate (FDR) under unknown arbitrarily interdependent p-values, without explicitly modeling these inter-correlations. They include the FWER-controlling Bonferroni MTP and Holm MTP; the FDR-controlling Benjamini-Yekutieli MTP; and the DP-MTP, based on a Dirichlet process (DP) prior distribution supporting the entire space of MTPs that control either the FWER or the FDR. For such an MTP, this study introduces a new and congenial method for Bayesian predictive power analysis, for power calculation and sample size determination for any given planned future (e.g., replication or interim) study. This novel MTP predictive power analysis method is based on a joint prior distribution defining a scale matrix mixture of asymmetric multivariate normal mean-variance mixture distributions, factorized as a general prior distribution for effect sizes (e.g., obtained from expert judgment or results of prior studies), and a uniform prior distribution for correlation matrices representing arbitrary dependencies between p-values of test statistics of given multiple hypothesis tests under their alternative hypotheses. The new MTP power analysis method also results in p-value weights which can be used to minimize the relative impacts of and assess for significance-chasing biases (e.g., publication bias, p-hacking, etc.) in multiple testing, without needing to assume that p-values (effect sizes) are independent. Previous MTP power analysis methods are conditional in that they assume fixed effect sizes and a fixed correlation matrix for test statistics (p-values) which may be difficult to specify; while, incongenially, not fully accounting for their uncertainty, especially for MTPs that control either the FWER or FDR for arbitrarily correlated p-values. The new simulation-based MTP predictive power analysis method is illustrated through the analysis of p-values obtained by a famous study of lead exposure and re-analyzed by the previous MTP literature, using R package bnpMTP (v1.0).

This paper introduces a sensitivity analysis method for multiple testing procedures (MTPs) using marginal ‐values. The method is based on the Dirichlet process (DP) prior distribution, specified to support the entire space of MTPs, where each MTP controls either the family‐wise error rate (FWER) or the false discovery rate (FDR) under arbitrary dependence between ‐values. This DP‐MTP sensitivity analysis method provides uncertainty quantification for MTPs, by accounting for uncertainty in the selection of such MTPs and their respective threshold‐based decisions regarding which number of smallest ‐values are significant discoveries, from a given set of null hypothesis tested, while measuring each ‐value's probability of significance over the DP prior predictive distribution of this space of all MTPs, and reducing the possible conservativeness of using only one such MTP for multiple testing. The DP‐MTP sensitivity analysis method is illustrated through the analysis of over 28,000 ‐values arising from hypothesis tests performed on a 2022 dataset of a representative sample of three million U.S. high school students observed on 239 variables. They include tests which, respectively, relate variables about the disruption caused by school closures during the COVID‐19 pandemic, with various mathematical cognition, academic achievement, and student background variables. R software code for the DP‐MTP sensitivity analysis method is provided in the Code and Data Supplement (CDS) of this paper.

The article makes important new contributions to the growing and large literature on the analysis (and correction for) response style bias(es) on Likert scale items, as it clearly describes. For this analytical purpose, the article describes a nonlin...

Ongoing modern computational advancements continue to make it easier to collect increasingly large and complex datasets, which can often only be realistically analyzed using models defined by intractable likelihood functions. This Stats invited feature article introduces and provides an extensive simulation study of a new approximate Bayesian computation (ABC) framework for estimating the posterior distribution and the maximum likelihood estimate (MLE) of the parameters of models defined by intractable likelihoods, that unifies and extends previous ABC methods proposed separately. This framework, copulaABCdrf, aims to accurately estimate and describe the possibly skewed and high-dimensional posterior distribution by a novel multivariate copula-based meta-t distribution based on univariate marginal posterior distributions that can be accurately estimated by distribution random forests (drf), while performing automatic summary statistics (covariates) selection, based on robustly estimated copula dependence parameters. The copulaABCdrf framework also provides a novel multivariate mode estimator to perform MLE and posterior mode estimation and an optional step to perform model selection from a given set of models using posterior probabilities estimated by drf. The posterior distribution estimation accuracy of the ABC framework is illustrated and compared with previous standard ABC methods through several simulation studies involving low- and high-dimensional models with computable posterior distributions, which are either unimodal, skewed, or multimodal; and exponential random graph and mechanistic network models, each defined by an intractable likelihood from which it is costly to simulate large network datasets. This paper also proposes and studies a new solution to the simulation cost problem in ABC involving the posterior estimation of parameters from datasets simulated from the given model that are smaller compared to the potentially large size of the dataset being analyzed. This proposal is motivated by the fact that, for many models defined by intractable likelihoods, such as the network models when they are applied to analyze massive networks, the repeated simulation of large datasets (networks) for posterior-based parameter estimation can be too computationally costly and vastly slow down or prohibit the use of standard ABC methods. The copulaABCdrf framework and standard ABC methods are further illustrated through analyses of large real-life networks of sizes ranging between 28,000 and 65.6 million nodes (between 3 million and 1.8 billion edges), including a large multilayer network with weighted directed edges. The results of the simulation studies show that, in settings where the true posterior distribution is not highly multimodal, copulaABCdrf usually produced similar point estimates from the posterior distribution for low-dimensional parametric models as previous ABC methods, but the copula-based method can produce more accurate estimates from the posterior distribution for high-dimensional models, and, in both dimensionality cases, usually produced more accurate estimates of univariate marginal posterior distributions of parameters. Also, posterior estimation accuracy was usually improved when pre-selecting the important summary statistics using drf compared to ABC employing no pre-selection of the subset of important summaries. For all ABC methods studied, accurate estimation of a highly multimodal posterior distribution was challenging. In light of the results of all the simulation studies, this article concludes by discussing how the copulaABCdrf framework can be improved for future research.

Tukey's halfspace depth is a theoretically-ideal depth function which defines multivariate order statistics by measuring the degree of centrality of each multivariate datapoint in its underlying distribution, and provides a broadly-applicable tool for multivariate nonparametric and robust data analysis. The exact computation of the halfspace depth is NP-hard, but it can be estimated by the random Tukey halfspace depth, based on drawing uniform random samples of direction vectors and resulting univariate projections of datapoints. However there is no method to determine the number of direction samples needed for accurate estimation based on any specific dataset given at hand. This article proposes to estimate the halfspace depth through an algorithm which computes the Monte Carlo confidence interval of the halfspace depth from overlapping-batch subsamples of random direction samples,while generating an increasingly larger number of direction samples, until the halfspace depth is estimated to the desired confidence level. The algorithm is easily extended to the simultaneous estimation of the halfspace depths of multiple multivariate data points. The viability of this subsampling halfspace depth algorithm is illustrated through the analysis of many real and simulated multivariate datasets. Software code for the algorithm and related methods are provided as online supplementary materials.

Approximate Bayesian Computation (ABC) can provide inferences from the (approximate) posterior distribution based on intractable likelihoods. The quality of ABC inferences relies on the choice of tolerance for the distance between the observed data summary statistics, and the pseudo-data summary statistics simulated from the likelihood, used within the context of an algorithm which samples from the approximate posterior. However, the ABC literature does not provide an automatic method to select the best tolerance level for the given dataset at hand, and in ABC practice finding the best tolerance level can be time consuming. This note introduces a fast automatic estimator of the tolerance, based on the parametric bootstrap. After the tolerance estimate is calculated, it can then be input into any suitable importance sampling or MCMC algorithm to approximate from the target approximate posterior distribution. This tolerance estimator is illustrated through ABC analyses of simulated and real datasets involving several intractable likelihood models. This includes the analysis of a real 23,000-node network dataset involving stochastic search model selection.

We propose Dirichlet Process Mixture (DPM) models for prediction and cluster-wise variable selection, based on two choices of shrinkage baseline prior distributions for the linear regression coefficients, namely the Horseshoe prior and Normal-Gamma prior. We show in a simulation study that each of the two proposed DPM models tend to outperform the standard DPM model based on the non-shrinkage normal prior, in terms of predictive, variable selection, and clustering accuracy. This is especially true for the Horseshoe model, and when the number of covariates exceeds the within-cluster sample size. A real data set is analyzed to illustrate the proposed modeling methodology, where both proposed DPM models again attained better predictive accuracy.

We propose Dirichlet process mixture (DPM) models for prediction and cluster‐wise variable selection, based on two choices of shrinkage baseline prior distributions for the linear regression coefficients, namely, the Horseshoe prior and the Normal‐Gamma prior. We show in a simulation study that each of the two proposed DPM models tends to outperform the standard DPM model based on the non‐shrinkage normal prior, in terms of predictive, variable selection, and clustering accuracy. This is especially true for the Horseshoe model and when the number of covariates exceeds the within‐cluster sample size. A real data set is analysed to illustrate the proposed modelling methodology, where both proposed DPM models again attained better predictive accuracy.

Bayesian nonparametric (BNP) infinite-mixture models provide flexible and accurate density estimation, cluster analysis, and regression. However, for the posterior inference of such a model, MCMC algorithms are complex, often need to be tailor-made for different BNP priors, and are intractable for large data sets. We introduce a BNP classification annealing EM algorithm which employs importance sampling estimation. This new fast-search algorithm, for virtually any given BNP mixture model, can quickly and accurately calculate the posterior predictive density estimate (by posterior averaging) and the maximum a-posteriori clustering estimate (by simulated annealing), even for data sets containing millions of observations. The algorithm can handle a wide range of BNP priors because it primarily relies on the ability to generate prior samples. The algorithm can be fast because in each iteration, it performs a sampling step for the (missing) clustering of the data points, instead of a costly E-step; and then performs direct posterior calculations in the M-step, given the sampled (imputed) clustering. The new algorithm is illustrated and evaluated through BNP Gaussian mixture model analyses of benchmark simulated data and real data sets. MATLAB code for the new algorithm is provided in Supplementary Materials.