Variable Selection Diagnostics Measures for High-Dimensional Regression

Markov chain Monte Carlo (MCMC) is nowadays a standard approach to numerical computation of integrals of the posterior density π of the parameter vector η. Unfortunately, Bayesian inference using MCMC is computationally intractable when the posterior density π is expensive to evaluate. In many such problems, it is possible to identify a minimal subvector β of η responsible for the expensive computation in the evaluation of π. We propose two approaches, DOSKA and INDA, that approximate π by interpolation in ways that exploit this computational structure to mitigate the curse of dimensionality. DOSKA interpolates π directly while INDA interpolates π indirectly by interpolating functions, for example, a regression function, upon which π depends. Our primary contribution is derivation of a Gaussian processes interpolant that provably improves over some of the existing approaches by reducing the effective dimension of the interpolation problem from dim(η) to dim(β). This allows a dramatic reduction of the number of expensive evaluations necessary to construct an accurate approximation of π when dim(η) is high but dim(β) is low.

We illustrate the proposed approaches in a case study for a spatio-temporal linear model for air pollution data in the greater Boston area.

Supplemental materials include proofs, details, and software implementation of the proposed procedures.     

Local Derivative-Free Approximation of Computationally Expensive Posterior Densities

Bayesian inference using Markov chain Monte Carlo (MCMC) is computationally prohibitive when the posterior density of interest, π, is computationally expensive to evaluate. We develop a derivative-free algorithm GRIMA to accurately approximate π by interpolation over its high-probability density (HPD) region, which is initially unknown. Our local approach reduces the waste of computational budget on approximation of π in the low-probability region, which is inherent in global experimental designs. However, estimation of the HPD region is nontrivial when derivatives of π are not available or are not informative about the shape of the HPD region. Without relying on derivatives, GRIMA iterates (a) sequential knot selection over the estimated HPD region of π to refine the surrogate posterior and (b) re-estimation of the HPD region using an MCMC sample from the updated surrogate density, which is inexpensive to obtain. GRIMA is applicable to approximation of general unnormalized posterior densities. To determine the range of tractable problem dimensions, we conduct simulation experiments on test densities with linear and nonlinear component-wise dependence, skewness, kurtosis and multimodality. Subsequently, we use GRIMA in a case study to calibrate a computationally intensive nonlinear regression model to real data from the Town Brook watershed. Supplemental materials for this article are available online.     

Markov chain Monte Carlo Using an Approximation

This article presents a method for generating samples from an unnormalized posterior distribution f(·) using Markov chain Monte Carlo (MCMC) in which the evaluation of f(·) is very difficult or computationally demanding. Commonly, a less computationally demanding, perhaps local, approximation to f(·) is available, say f^**_x(·). An algorithm is proposed to generate an MCMC that uses such an approximation to calculate acceptance probabilities at each step of a modified Metropolis–Hastings algorithm. Once a proposal is accepted using the approximation, f(·) is calculated with full precision ensuring convergence to the desired distribution. We give sufficient conditions for the algorithm to converge to f(·) and give both theoretical and practical justifications for its usage. Typical applications are in inverse problems using physical data models where computing time is dominated by complex model simulation. We outline Bayesian inference and computing for inverse problems. A stylized example is given of recovering resistor values in a network from electrical measurements made at the boundary. Although this inverse problem has appeared in studies of underground reservoirs, it has primarily been chosen for pedagogical value because model simulation has precisely the same computational structure as a finite element method solution of the complete electrode model used in conductivity imaging, or “electrical impedance tomography.” This example shows a dramatic decrease in CPU time, compared to a standard Metropolis–Hastings algorithm.     

An Efficient Algorithm for Bayesian Nearest Neighbours

K-Nearest Neighbours (k-NN) is a popular classification and regression algorithm, yet one of its main limitations is the difficulty in choosing the number of neighbours. We present a Bayesian algorithm to compute the posterior probability distribution for k given a target point within a data-set, efficiently and without the use of Markov Chain Monte Carlo (MCMC) methods or simulation—alongside an exact solution for distributions within the exponential family. The central idea is that data points around our target are generated by the same probability distribution, extending outwards over the appropriate, though unknown, number of neighbours. Once the data is projected onto a distance metric of choice, we can transform the choice of k into a change-point detection problem, for which there is an efficient solution: we recursively compute the probability of the last change-point as we move towards our target, and thus de facto compute the posterior probability distribution over k. Applying this approach to both a classification and a regression UCI data-sets, we compare favourably and, most importantly, by removing the need for simulation, we are able to compute the posterior probability of k exactly and rapidly. As an example, the computational time for the Ripley data-set is a few milliseconds compared to a few hours when using a MCMC approach.

     

Improving the Data Augmentation Algorithm in the Two-Block Setup

The data augmentation (DA) approach to approximate sampling from an intractable probability density f_X is based on the construction of a joint density, f_{X, Y}, whose conditional densities, f_X|Y and f_Y|X, can be straightforwardly sampled. However, many applications of the DA algorithm do not fall in this “single-block” setup. In these applications, X is partitioned into two components, X = (U, V), in such a way that it is easy to sample from f_Y|X, f_{U|V, Y}, and f_{V|U, Y}. We refer to this alternative version of DA, which is effectively a three-variable Gibbs sampler, as “two-block” DA. We develop two methods to improve the performance of the DA algorithm in the two-block setup. These methods are motivated by the Haar PX-DA algorithm, which has been developed in previous literature to improve the performance of the single-block DA algorithm. The Haar PX-DA algorithm, which adds a computationally inexpensive extra step in each iteration of the DA algorithm while preserving the stationary density, has been shown to be optimal among similar techniques. However, as we illustrate, the Haar PX-DA algorithm does not lead to the required stationary density f_X in the two-block setup. Our methods incorporate suitable generalizations and modifications to this approach, and work in the two-block setup. A theoretical comparison of our methods to the two-block DA algorithm, a much harder task than the single-block setup due to nonreversibility and structural complexities, is provided. We successfully apply our methods to applications of the two-block DA algorithm in Bayesian robit regression and Bayesian quantile regression. Supplementary materials for this article are available online.     

Recentered Importance Sampling With Applications to Bayesian Model Validation

Since its introduction in the early 1990s, the idea of using importance sampling (IS) with Markov chain Monte Carlo (MCMC) has found many applications. This article examines problems associated with its application to repeated evaluation of related posterior distributions with a particular focus on Bayesian model validation. We demonstrate that, in certain applications, the curse of dimensionality can be reduced by a simple modification of IS. In addition to providing new theoretical insight into the behavior of the IS approximation in a wide class of models, our result facilitates the implementation of computationally intensive Bayesian model checks. We illustrate the simplicity, computational savings, and potential inferential advantages of the proposed approach through two substantive case studies, notably computation of Bayesian p-values for linear regression models and simulation-based model checking. Supplementary materials including the Appendix and the R code for Section are available online.     

Maximum Likelihood for Generalized Linear Models with Nested Random Effects via High-Order,Multivariate Laplace Approximation

Abstract

Nested random effects models are often used to represent similar processes occurring in each of many clusters. Suppose that, given cluster-specific random effects b, the data y are distributed according to f(y|b, Θ), while b follows a density p(b|Θ). Likelihood inference requires maximization of ∫ f(y|b, Θ)p(b|Θdb with respect to Θ. Evaluation of this integral often proves difficult, making likelihood inference difficult to obtain. We propose a multivariate Taylor series approximation of the log of the integrand that can be made as accurate as desired if the integrand and all its partial derivatives with respect to b are continuous in the neighborhood of the posterior mode of b|Θ,y. We then apply a Laplace approximation to the integral and maximize the approximate integrated likelihood via Fisher scoring. We develop computational formulas that implement this approach for two-level generalized linear models with canonical link and multivariate normal random effects. A comparison with approximations based on penalized quasi-likelihood, Gauss—Hermite quadrature, and adaptive Gauss-Hermite quadrature reveals that, for the hierarchical logistic regression model under the simulated conditions, the sixth-order Laplace approach is remarkably accurate and computationally fast.     

Thank God That Regressing Y on X is Not the Same as Regressing X on Y: Direct and Indirect Residual Augmentations

What does regressing Y on X versus regressing X on Y have to do with Markov chain Monte Carlo (MCMC)? It turns out that many strategies for speeding up data augmentation (DA) type algorithms can be understood as fostering independence or “de-correlation” between a regression function and the corresponding residual, thereby reducing or even eliminating dependence among MCMC iterates. There are two general classes of algorithms, those corresponding to regressing parameters on augmented data/auxiliary variables and those that operate the other way around. The interweaving strategy of Yu and Meng provides a general recipe to automatically take advantage of both, and it is the existence of two different types of residuals that makes the interweaving strategy seemingly magical in some cases and promising in general. The concept of residuals—which depends on actual data—also highlights the potential for substantial improvements when DA schemes are allowed to depend on the observed data. At the same time, there is an intriguing phase transition type of phenomenon regarding choosing (partially) residual augmentation schemes, reminding us once more of the prevailing issue of trade-off between robustness and efficiency. This article reports on these latest theoretical investigations (using a class of normal/independence models) and empirical findings (using a posterior sampling for a probit regression) in the search for effective residual augmentations—and ultimately more MCMC algorithms—that meet the 3-S criterion: simple, stable, and speedy. Supplementary materials for the article are available online.     

FUNCTORIAL PROPERTIES OF THE HOPF INVARIANT I

Abstract

Homotopy operations Θ: [ΣY, U] → [ΣY, V] which are natural in Y are considered. In particular a technique used in the definition of the Hopf invariant (as treated by Berstein-Hilton) shows that any fibration p: E → B with fiber V, when provided with a homotopy section of Ωp, determines such a homotopy operation [ΣY, E] → [ΣY, V]. More generally, starting from a track class of homotopies α º f ? β º g we adapt this fibration technique to construct a homotopy operation [ΣY, M(f,g)] → [ΣY, F _α * F _β] called a Hopf invariant. The intervening fibration in the definition of this Hopf invariant arises via the fiberwise join construction.     

Legal coloring of graphs

The following computational problem was initiated by Manber and Tompa (22nd FOCS Conference, 1981): Given a graphG=(V, E) and a real functionf:V→R which is a proposed vertex coloring. Decide whetherf is a proper vertex coloring ofG. The elementary steps are taken to be linear comparisons. Lower bounds on the complexity of this problem are derived using the chromatic polynomial ofG. It is shown how geometric parameters of a space partition associated withG influence the complexity of this problem. Existing methods for analyzing such space partitions are suggested as a powerful tool for establishing lower bounds for a variety of computational problems.     

Confidence Intervals for Discrete Approximations to Ill-Posed Problems

We consider the linear model Y = Xβ + ε that is obtained by discretizing a system of first-kind integral equations describing a set of physical measurements. The n vector β represents the desired quantities, the m x n matrix X represents the instrument response functions, and the m vector Y contains the measurements actually obtained. These measurements are corrupted by random measuring errors ε drawn from a distribution with zero mean vector and known variance matrix. Solution of first-kind integral equations is an ill-posed problem, so the least squares solution for the above model is a highly unstable function of the measurements, and the classical confidence intervals for the solution are too wide to be useful. The solution can often be stabilized by imposing physically motivated nonnegativity constraints. In a previous article (O'Leary and Rust 1986) we developed a method for computing sets of nonnegatively constrained simultaneous confidence intervals. In this article we briefly review the simultaneous intervals and then show how to compute nonnegativity constrained one-at-a-time confidence intervals. The technique gives valid confidence intervals even for problems with m < n. We demonstrate the methods using both an overdetermined and an underdetermined problem obtained by discretizing an equation of Phillips (Phillips 1962).     

Generalized time-dependent conditional linear models under left truncation and right censoring

Consider the model f(S(z|X)){\phi(S(z|X))} = \pmbb(z) [(X)\vec]{\pmb{\beta}(z) {\vec{X}}}, where f{\phi} is a known link function, S(·|X) is the survival function of a response Y given a covariate X, [(X)\vec]{\vec{X}} = (1, X, X ² , . . . , X ^p ) and \pmbb(z){\pmb{\beta}(z)} is an unknown vector of time-dependent regression coefficients. The response is subject to left truncation and right censoring. Under this model, which reduces for special choices of f{\phi} to e.g. Cox proportional hazards model or the additive hazards model with time dependent coefficients, we study the estimation of the vector \pmbb(z){\pmb{\beta}(z)} . A least squares approach is proposed and the asymptotic properties of the proposed estimator are established. The estimator is also compared with a competing maximum likelihood based estimator by means of simulations. Finally, the method is applied to a larynx cancer data set.     

Bivariate Factorizations via Galois Theory, with Application to Exceptional Polynomials

We present a method for factoring polynomials of the shapef(X) − f(Y), wherefis a univariate polynomial over a fieldk. We then apply this method in the case whenfis a member of the infinite family of exceptional polynomials we discovered jointly with H. Lenstra in 1995; factoringf(X) − f(Y) in this case was posed as a problem by S. Cohen shortly after the discovery of these polynomials.     

NONLINEAR WAVELET SMOOTHING OF ERROR DENSITY IN A SEMIPARAMETRIC REGRESSION MODEL

In this paper,a semlparametrie resresaion model in which errors are i. i. d random variables from an unknown density f( ) is considered. Based on Hall et al. (1995),a nonlinear wavelet estimation of f( ) without restrictions of continuity everywhere on f( ) is given,and the convergence rate of the estimators in L2 is obtained.     

Boundary Bias Correction for Nonparametric Deconvolution

In this paper we consider the deconvolution problem in nonparametric density estimation. That is, one wishes to estimate the unknown density of a random variable X, say f _X , based on the observed variables Y's, where Y = X + with being the error. Previous results on this problem have considered the estimation of f _X at interior points. Here we study the deconvolution problem for boundary points. A kernel-type estimator is proposed, and its mean squared error properties, including the rates of convergence, are investigated for supersmooth and ordinary smooth error distributions. Results of a simulation study are also presented.     

Asymptotically Optimal Tests and Optimal Designs for Testing the Mean in Regression Models with Applications to Change-Point Problems

Let a linear regression model be given with an experimental region [a, b] R and regression functions f ₁, ..., f _d+1 : [a, b] R. In practice it is an important question whether a certain regression function f _d+1, say, does or does not belong to the model. Therefore, we investigate the test problem H ₀ : "f _d+1 does not belong to the model" against K : "f _d+1 belong to the model" based on the least-squares residuals of the observations made at design points of the experimental region [a, b]. By a new functional central limit theorem given in Bischoff (1998, Ann. Statist. 26, 1398–1410), we are able to determine optimal tests in an asymptotic way. Moreover, we introduce the problem of experimental design for the optimal test statistics. Further, we compare the asymptotically optimal test with the likelihood ratio test (F-test) under the assumption that the error is normally distributed. Finally, we consider real change-point problems as examples and investigate by simulations the behavior of the asymptotic test for finite sample sizes. We determine optimal designs for these examples.     

Bayesian calibration in the estimation of the age of rhinoceros

In this paper the Bayesian approach for nonlinear multivariate calibration will be illustrated. This goal will be achieved by applying the Gibbs sampler to the rhinoceros data given by Clarke (1992, Biometrics, 48(4), 1081–1094). It will be shown that the point estimates obtained from the profile likelihoods and those calculated from the marginal posterior densities using improper priors will in most cases be similar.     

Some Maximal Function Fields and Additive Polynomials

We derive explicit equations for the maximal function fields F over 𝔽 _{q ²ⁿ} given by F = 𝔽 _{q ²ⁿ} (X, Y) with the relation A(Y) = f(X), where A(Y) and f(X) are polynomials with coefficients in the finite field 𝔽 _{q ²ⁿ} , and where A(Y) is q-additive and deg(f) = q ⁿ  + 1. We prove in particular that such maximal function fields F are Galois subfields of the Hermitian function field H over 𝔽 _{q ²ⁿ} (i.e., the extension H/F is Galois).     

A GENERALIZATION OF AUSLANDER–BUCHSBAUM AND BASS FORMULAS

Abstract

Given a contravariant functor F : 𝒞 → 𝒮ets for some category 𝒞, we say that F (𝒞) (or F) is generated by a pair (X, x) where X is an object of 𝒞 and x ∈ F(X) if for any object Y of 𝒞 and any y ∈ F(Y), there is a morphism f : Y → X such that F(f)(x) = y. Furthermore, when Y = X and y = x, any f : X → X such that F(f)(x) = x is an automorphism of X, we say that F is minimally generated by (X, x). This paper shows that if the ring R is left noetherian, then there exists a minimal generator for the functor ?xt (?, M) : ? → 𝒮ets, where M is a left R-module and ? is the class (considered as full subcategory of left R-modules) of injective left R-modules.     

Generalized essential variation and quasi-continuity

Summary For functions f : E →Y, where E ˆ R is a semi-open set and Y is a metric space, the essential φ-variation V_φ,&weierp; (f, E) depending on some σ-ideal in R is defined. In the case where f is quasi-continuous, V_{φ, &weierp;} (f, E) is equal to the φ-variation V_φ (f, E). If in addition Y is complete, then each function f of bounded φ -variation is graph quasi-continuous.