Regression Discontinuity Design (RDD) is a popular framework for estimating a causal effect in settings where treatment is assigned if an observed covariate exceeds a fixed threshold. We consider estimation and inference in the common setting where the sample consists of multiple known sub-populations with potentially heterogeneous treatment effects. In the applied literature, it is common to account for heterogeneity by either fitting a parametric model or considering each sub-population separately. In contrast, we develop a Bayesian hierarchical model using Gaussian process regression which allows for non-parametric regression while borrowing information across sub-populations. We derive the posterior distribution, prove posterior consistency, and develop a Metropolis-Hastings within Gibbs sampling algorithm. In extensive simulations, we show that the proposed procedure outperforms existing methods in both estimation and inferential tasks. Finally, we apply our procedure to U.S. Senate election data and discover an incumbent party advantage which is heterogeneous over different time periods.

Cosmic demographics—the statistical study of populations of astrophysical objects—has long relied on tools from multivariate statistics for analyzing data comprising fixed-length vectors of properties of objects, as might be compiled in a tabular astronomical catalog (say, with sky coordinates, and brightness measurements in a fixed number of spectral passbands). But beginning with the emergence of automated digital sky surveys, ca. 2000, astronomers began producing large collections of data with more complex structures: light curves (brightness time series) and spectra (brightness vs. wavelength). These comprise what statisticians call functional data—measurements of populations of functions. Upcoming automated sky surveys will soon provide astronomers with a flood of functional data. New methods are needed to accurately and optimally analyze large ensembles of light curves and spectra, accumulating information both along individual measured functions and across a population of such functions. Functional data analysis (FDA) provides tools for statistical modeling of functional data. Astronomical data presents several challenges for FDA methodology, e.g., sparse, irregular, and asynchronous sampling, and heteroscedastic measurement error. Bayesian FDA uses hierarchical Bayesian models for function populations, and is well suited to addressing these challenges. We provide an overview of astronomical functional data and some key Bayesian FDA modeling approaches, including functional mixed effects models, and stochastic process models. We briefly describe a Bayesian FDA framework combining FDA and machine learning methods to build low-dimensional parametric models for galaxy spectra.

Functional principal components analysis is a popular tool for inference on functional data. Standard approaches rely on an eigendecomposition of a smoothed covariance surface in order to extract the orthonormal functions representing the major modes of variation. This approach can be a computationally intensive procedure, especially in the presence of large datasets with irregular observations. In this article, we develop a Bayesian approach, which aims to determine the Karhunen-Lo\`eve decomposition directly without the need to smooth and estimate a covariance surface. More specifically, we develop a variational Bayesian algorithm via message passing over a factor graph, which is more commonly referred to as variational message passing. Message passing algorithms are a powerful tool for compartmentalizing the algebra and coding required for inference in hierarchical statistical models. Recently, there has been much focus on formulating variational inference algorithms in the message passing framework because it removes the need for rederiving approximate posterior density functions if there is a change to the model. Instead, model changes are handled by changing specific computational units, known as fragments, within the factor graph. We extend the notion of variational message passing to functional principal components analysis. Indeed, this is the first article to address a functional data model via variational message passing. Our approach introduces two new fragments that are necessary for Bayesian functional principal components analysis. We present the computational details, a set of simulations for assessing accuracy and speed and an application to United States temperature data.

This paper addresses the deconvolution problem of estimating a square-integrable probability density from observations contaminated with additive measurement errors having a known density. The estimator begins with a density estimate of the contaminated observations and minimizes a reconstruction error penalized by an integrated squared $m$-th derivative. Theory for deconvolution has mainly focused on kernel- or wavelet-based techniques, but other methods including spline-based techniques and this smoothness-penalized estimator have been found to outperform kernel methods in simulation studies. This paper fills in some of these gaps by establishing asymptotic guarantees for the smoothness-penalized approach. Consistency is established in mean integrated squared error, and rates of convergence are derived for Gaussian, Cauchy, and Laplace error densities, attaining some lower bounds already in the literature. The assumptions are weak for most results; the estimator can be used with a broader class of error densities than the deconvoluting kernel. Our application example estimates the density of the mean cytotoxicity of certain bacterial isolates under random sampling; this mean cytotoxicity can only be measured experimentally with additive error, leading to the deconvolution problem. We also describe a method for approximating the solution by a cubic spline, which reduces to a quadratic program.

ABSTRACT More than 5000 extrasolar planets have already been detected. JWST and near-term ground-based telescopes like the Extremely Large Telescope (ELT), Giant Magellan Telescope (GMT), Thirty Meter Telescope (TMT), and upcoming telescopes such as the Nancy Grace Roman Space Telescope, Xuntian, and Ariel are designed to characterize the atmosphere of directly imaged Jovian planets. Here, we used five diverse machine learning algorithms to investigate how well broad-band filter photometric fluxes could initially characterize giant exoplanets. We use an established grid of 8813 reflected light model spectra of different metallicities, planet–star distances, and cloud properties to assess the performance of several machine learning algorithms on both noiseless and noisy data to provide classification and regression results as a function of signal to noise of the data. In all cases, the algorithms were tested on noisy validation data. The results show that the use of machine learning to characterize giant planets from reflected broad-band filter photometry provides a promising tool for initial characterization, with over 65 per cent accuracy in characterizing metallicity for signal-to-noise ratios (S/N) ≳ 30, over 80 per cent for cloud coverage for S/N ≳ 30. This approach will allow initial characterization for large surveys of giant exoplanets and prioritization for spectroscopy observations of a subset of these worlds.

We propose a projection-based test to check logistic regression models when the dimension of the covariate vector may be divergent. The proposed test achieves a reduction in dimension, and the proposed method behaves as if only a single covariate is present. The test is shown to be consistent and can detect root-n local alternatives. We derive the asymptotic distribution of the proposed test under the null hypothesis and establish the test's asymptotic behavior under the local and global alternatives. The numerical performance is remarkably attractive comparing to the existing methods. Real examples are presented for illustration. Supplementary materials for this article are available online.

In high-dimensional prediction problems, we propose subsampling the predictors prior to the analysis. Specifically, we draw features using random sampling, and then fit a model and make predictions based on the sampled feature subset. This greatly reduces the dimension, storage, and computational bottlenecks. We explore this "subset regression" strategy under a linear regression framework. We propose an ensemble method that combines multiple subset regressions, called the ensemble subset regression (ENSURE) that reduces the uncertainty due to feature sampling. We provide a theoretical upper bound on the excess risk of the predictions computed in the subset regression, and provide theoretical support that the ensemble can improve the performance of the subset regression. Detailed empirical studies demonstrate that ENSURE performs well, better than methods that use all features.

We present a new functional Bayes classifier that uses principal component (PC) or partial least squares (PLS) scores from the common covariance function, that is, the covariance function marginalized over groups. When the groups have different covariance functions, the PC or PLS scores need not be independent or even uncorrelated. We use copulas to model the dependence. Our method is semiparametric; the marginal densities are estimated nonparametrically by kernel smoothing and the copula is modeled parametrically. We focus on Gaussian and t-copulas, but other copulas could be used. The strong performance of our methodology is demonstrated through simulation, real data examples, and asymptotic properties.

Background Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) is a complex, heterogenous disease characterized by unexplained persistent fatigue and other features including cognitive impairment, myalgias, post-exertional malaise, and immune system dysfunction. Cytokines are present in plasma and encapsulated in extracellular vesicles (EVs), but there have been only a few reports of EV characteristics and cargo in ME/CFS. Several small studies have previously described plasma proteins or protein pathways that are associated with ME/CFS. Methods We prepared extracellular vesicles (EVs) from frozen plasma samples from a cohort of Myalgic Encephalomyelitis/Chronic Fatigue Syndrome (ME/CFS) cases and controls with prior published plasma cytokine and plasma proteomics data. The cytokine content of the plasma-derived extracellular vesicles was determined by a multiplex assay and differences between patients and controls were assessed. We then performed multi-omic statistical analyses that considered not only this new data, but extensive clinical data describing the health of the subjects. Results ME/CFS cases exhibited greater size and concentration of EVs in plasma. Assays of cytokine content in EVs revealed IL2 was significantly higher in cases. We observed numerous correlations among EV cytokines, among plasma cytokines, and among plasma proteins from mass spectrometry proteomics. Significant correlations between clinical data and protein levels suggest roles of particular proteins and pathways in the disease. For example, higher levels of the pro-inflammatory cytokines Granulocyte-Monocyte Colony-Stimulating Factor (CSF2) and Tumor Necrosis Factor (TNFα) were correlated with greater physical and fatigue symptoms in ME/CFS cases. Higher serine protease SERPINA5, which is involved in hemostasis, was correlated with higher SF-36 general health scores in ME/CFS. Machine learning classifiers were able to identify a list of 20 proteins that could discriminate between cases and controls, with XGBoost providing the best classification with 86.1% accuracy and a cross-validated AUROC value of 0.947. Random Forest distinguished cases from controls with 79.1% accuracy and an AUROC value of 0.891 using only 7 proteins. Conclusions These findings add to the substantial number of objective differences in biomolecules that have been identified in individuals with ME/CFS. The observed correlations of proteins important in immune responses and hemostasis with clinical data further implicates a disturbance of these functions in ME/CFS.

Survey-based measurements of the spectral energy distributions (SEDs) of galaxies have flux density estimates on badly misaligned grids in rest-frame wavelength. The shift to rest frame wavelength also causes estimated SEDs to have differing support. For many galaxies, there are sizeable wavelength regions with missing data. Finally, dim galaxies dominate typical samples and have noisy SED measurements, many near the limiting signal-to-noise level of the survey. These limitations of SED measurements shifted to the rest frame complicate downstream analysis tasks, particularly tasks requiring computation of functionals (e.g., weighted integrals) of the SEDs, such as synthetic photometry, quantifying SED similarity, and using SED measurements for photometric redshift estimation. We describe a hierarchical Bayesian framework, drawing on tools from functional data analysis, that models SEDs as a random superposition of smooth continuum basis functions (B-splines) and line features, comprising a finite-rank, nonstationary Gaussian process, measured with additive Gaussian noise. We apply this *Splines 'n Lines* (SnL) model to a collection of 678,239 galaxy SED measurements comprising the Main Galaxy Sample from the Sloan Digital Sky Survey, Data Release 17, demonstrating capability to provide continuous estimated SEDs that reliably denoise, interpolate, and extrapolate, with quantified uncertainty, including the ability to predict line features where there is missing data by leveraging correlations between line features and the entire continuum.

We construct the maximally predictable portfolio (MPP) of stocks using machine learning. Solving for the optimal constrained weights in the multi-asset MPP gives portfolios with a high monthly coefficient of determination, given the sample covariance matrix of predicted return errors from a machine learning model. Various models for the covariance matrix are tested. The MPPs of S&P 500 index constituents with estimated returns from Elastic Net, Random Forest, and Support Vector Regression models can outperform or underperform the index depending on the time period. Portfolios that take advantage of the high predictability of the MPP's returns and employ a Kelly criterion style strategy consistently outperform the benchmark.

In this article, we develop uniform inference methods for the conditional mode based on quantile regression. Specifically, we propose to estimate the conditional mode by minimizing the derivative of the estimated conditional quantile function defined by smoothing the linear quantile regression estimator, and develop two bootstrap methods, a novel pivotal bootstrap and the nonparametric bootstrap, for our conditional mode estimator. Building on high-dimensional Gaussian approximation techniques, we establish the validity of simultaneous confidence rectangles constructed from the two bootstrap methods for the conditional mode. We also extend the preceding analysis to the case where the dimension of the covariate vector is increasing with the sample size. Finally, we conduct simulation experiments and a real data analysis using the U.S. wage data to demonstrate the finite sample performance of our inference method. The supplemental materials include the wage dataset, R codes and an appendix containing proofs of the main results, additional simulation results, discussion of model misspecification and quantile crossing, and additional details of the numerical implementation.

Regression models that ignore measurement error in predictors may produce highly biased estimates leading to erroneous inferences. It is well known that it is extremely difficult to take measurement error into account in Gaussian non-parametric regression. This problem becomes even more difficult when considering other families such as binary, Poisson and negative binomial regression. We present a novel method aiming to correct for measurement error when estimating regression functions. Our approach is sufficiently flexible to cover virtually all distributions and link functions regularly considered in generalised linear models. This approach depends on approximating the first and the second moment of the response after integrating out the true unobserved predictors in any semi-parametric generalised regression model. By the latter is meant a model with both linear and non-parametric effects that are connected to the mean response by a link function and with a response distribution in an exponential family or quasi-likelihood model. Unlike previous methods, the method we now propose is not restricted to truncated splines and can utilise various basis functions. Moreover, it can operate without making any distributional assumption about the unobserved predictor. Through extensive simulation studies, we study the performance of our method under many scenarios.

We develop a generalized partially additive model to build a single semiparametric risk scoring system for physical activity across multiple populations. A score comprised of distinct and objective physical activity measures is a new concept that offers challenges due to the nonlinear relationship between physical behaviors and various health outcomes. We overcome these challenges by modeling each score component as a smooth term, an extension of generalized partially linear single-index models. We use penalized splines and propose two inferential methods, one using profile likelihood and a nonparametric bootstrap, the other using a full Bayesian model, to solve additional computational problems. Both methods exhibit similar and accurate performance in simulations. These models are applied to the National Health and Nutrition Examination Survey and quantify nonlinear and interpretable shapes of score components for all-cause mortality.

We find economically and statistically significant gains when using machine learning for portfolio allocation between the market index and risk-free asset. Optimal portfolio rules for time-varying expected returns and volatility are implemented with two Random Forest models. One model is employed in forecasting the sign probabilities of the excess return with payout yields. The second is used to construct an optimized volatility estimate. Reward-risk timing with machine learning provides substantial improvements over the buy-and-hold in utility, risk-adjusted returns, and maximum drawdowns. This paper presents a new theoretical basis and unifying framework for machine learning applied to both return- and volatility-timing.

This paper studies a \textit{partial functional partially linear single-index model} that consists of a functional linear component as well as a linear single-index component. This model generalizes many well-known existing models and is suitable for more complicated data structures. However, its estimation inherits the difficulties and complexities from both components and makes it a challenging problem, which calls for new methodology. We propose a novel profile B-spline method to estimate the parameters by approximating the unknown nonparametric link function in the single-index component part with B-spline, while the linear slope function in the functional component part is estimated by the functional principal component basis. The consistency and asymptotic normality of the parametric estimators are derived, and the global convergence of the proposed estimator of the linear slope function is also established. More excitingly, the latter convergence is optimal in the minimax sense. A two-stage procedure is implemented to estimate the nonparametric link function, and the resulting estimator possesses the optimal global rate of convergence. Furthermore, the convergence rate of the mean squared prediction error for a predictor is also obtained. Empirical properties of the proposed procedures are studied through Monte Carlo simulations. A real data example is also analyzed to illustrate the power and flexibility of the proposed methodology.

Regression models that ignore measurement error in predictors may produce highly biased estimates leading to erroneous inferences. It is well known that it is extremely difficult to take measurement error into account in Gaussian nonparametric regression. This problem becomes tremendously more difficult when considering other families such as logistic regression, Poisson and negative-binomial. For the first time, we present a method aiming to correct for measurement error when estimating regression functions flexibly covering virtually all distributions and link functions regularly considered in generalized linear models. This approach depends on approximating the first and the second moment of the response after integrating out the true unobserved predictors in a semiparametric generalized linear model. Unlike previous methods, this method is not restricted to truncated splines and can utilize various basis functions. Through extensive simulation studies, we study the performance of our method under many scenarios.

Under "measurement constraints," responses are expensive to measure and initially unavailable on most of records in the dataset, but the covariates are available for the entire dataset. Our goal is to sample a relatively small portion of the dataset where the expensive responses will be measured and the resultant sampling estimator is statistically efficient. Measurement constraints require the sampling probabilities can only depend on a very small set of the responses. A sampling procedure that uses responses at most only on a small pilot sample will be called "response-free." We propose a response-free sampling procedure optimal sampling under measurement constraints (OSUMC) for generalized linear models. Using the A-optimality criterion, that is, the trace of the asymptotic variance, the resultant estimator is statistically efficient within a class of sampling estimators. We establish the unconditional asymptotic distribution of a general class of response-free sampling estimators. This result is novel compared with the existing conditional results obtained by conditioning on both covariates and responses. Under our unconditional framework, the subsamples are no longer independent and new martingale techniques are developed for our asymptotic theory. We further derive the A-optimal response-free sampling distribution. Since this distribution depends on population level quantities, we propose the OSUMC algorithm to approximate the theoretical optimal sampling. Finally, we conduct an intensive empirical study to demonstrate the advantages of OSUMC algorithm over existing methods in both statistical and computational perspectives. We find that OSUMC's performance is comparable to that of sampling algorithms that use complete responses. This shows that, provided an efficient algorithm such as OSUMC is used, there is little or no loss in accuracy due to the unavailability of responses because of measurement constraints. Supplementary materials for this article are available online.