Modelling and smoothing parameter estimation with multiple quadratic penalties

Penalized likelihood methods provide a range of practical modelling tools, including spline smoothing, generalized additive models and variants of ridge regression. Selecting the correct weights for penalties is a critical part of using these methods and in the single-penalty case the analyst has several well-founded techniques to choose from. However, many modelling problems suggest a formulation employing multiple penalties, and here general methodology is lacking. A wide family of models with multiple penalties can be fitted to data by iterative solution of the generalized ridge regression problem minimize || W ^1/2 ( Xp − y ) ||²ρ+Σ _{i =1} ^m  θ _i p ' S _i p ( p is a parameter vector, X a design matrix, S _i a non-negative definite coefficient matrix defining the i th penalty with associated smoothing parameter θ _i , W a diagonal weight matrix, y a vector of data or pseudodata and ρ an 'overall' smoothing parameter included for computational efficiency). This paper shows how smoothing parameter selection can be performed efficiently by applying generalized cross-validation to this problem and how this allows non-linear, generalized linear and linear models to be fitted using multiple penalties, substantially increasing the scope of penalized modelling methods. Examples of non-linear modelling, generalized additive modelling and anisotropic smoothing are given.     

A bounded influence regression estimator based on the statistics of the hat matrix

Summary. Many geophysical regression problems require the analysis of large (more than 10⁴ values) data sets, and, because the data may represent mixtures of concurrent natural processes with widely varying statistical properties, contamination of both response and predictor variables is common. Existing bounded influence or high breakdown point estimators frequently lack the ability to eliminate extremely influential data and/or the computational efficiency to handle large data sets. A new bounded influence estimator is proposed that combines high asymptotic efficiency for normal data, high breakdown point behaviour with contaminated data and computational simplicity for large data sets. The algorithm combines a standard M -estimator to downweight data corresponding to extreme regression residuals and removal of overly influential predictor values (leverage points) on the basis of the statistics of the hat matrix diagonal elements. For this, the exact distribution of the hat matrix diagonal elements p _ii for complex multivariate Gaussian predictor data is shown to be β ( p _ii ,  m ,  N − m ), where N is the number of data and m is the number of parameters. Real geophysical data from an auroral zone magnetotelluric study which exhibit severe outlier and leverage point contamination are used to illustrate the estimator's performance. The examples also demonstrate the utility of looking at both the residual and the hat matrix distributions through quantile–quantile plots to diagnose robust regression problems.     

Estimating smooth monotone functions

Many situations call for a smooth strictly monotone function f of arbitrary flexibility. The family of functions defined by the differential equation D  ² f  = w Df , where w is an unconstrained coefficient function comprises the strictly monotone twice differentiable functions. The solution to this equation is f = C ₀ + C ₁  D ⁻¹{exp( D ⁻¹ w )}, where C ₀ and C ₁ are arbitrary constants and D ⁻¹ is the partial integration operator. A basis for expanding w is suggested that permits explicit integration in the expression of f . In fitting data, it is also useful to regularize f by penalizing the integral of w ² since this is a measure of the relative curvature in f . Applications are discussed to monotone nonparametric regression, to the transformation of the dependent variable in non-linear regression and to density estimation.     

Non-parametric tests for distributional treatment effect for randomly censored responses

Summary.  For a binary treatment ν =0, 1 and the corresponding 'potential response' Y ⁰ for the control group ( ν =0) and Y ¹ for the treatment group ( ν =1), one definition of no treatment effect is that Y ⁰ and Y ¹ follow the same distribution given a covariate vector X . Koul and Schick have provided a non-parametric test for no distributional effect when the realized response (1− ν ) Y ⁰+ ν Y ¹ is fully observed and the distribution of X is the same across the two groups. This test is thus not applicable to censored responses, nor to non-experimental (i.e. observational) studies that entail different distributions of X across the two groups. We propose ' X -matched' non-parametric tests generalizing the test of Koul and Schick following an idea of Gehan. Our tests are applicable to non-experimental data with randomly censored responses. In addition to these motivations, the tests have several advantages. First, they have the intuitive appeal of comparing all available pairs across the treatment and control groups, instead of selecting a number of matched controls (or treated) in the usual pair or multiple matching. Second, whereas most matching estimators or tests have a non-overlapping support (of X ) problem across the two groups, our tests have a built-in protection against the problem. Third, Gehan's idea allows the tests to make good use of censored observations. A simulation study is conducted, and an empirical illustration for a job training effect on the duration of unemployment is provided.     

MINIMAX ESTIMATION OF A MULTIVARIATE NORMAL MEAN UNDER A CONVEX LOSS FUNCTION

Let X = (X₁, - X_p)prime; ˜ N_p (μ, Σ) where μ= (μ₁, -, μ_p)' and Σ= diag (Σ²₁, -, Σ²_p) are both unknown and p3. Let (n_i - 2) w_i/Σ²_i! X²_ni, independent. of w_i (I ≠ j = 1, -, p). Assume that (w₁, -, w_p) and X are independent. Define W = diag (w₁, -, w_p) and ¶ X ¶²_w= X'W^-1Q^-1W^-1X where Q = diag (q₁, -,n q_p), q_i > 0, i = 1, -, p. In this paper, the minimax estimator of Berger & Bock (1976), given by δ (X, W) = [I_p - r(X, W) ¶ X ¶^-2_w Q^-1W^-1] X, is shown to be minimax relative to the convex loss (δ - μ)'[αQ + (1 - α) Σ^-1] δ - μ)/C, where C =α tr (Σ) + (1 - α)p and 0 α 1, under certain conditions on r(X, W). This generalizes the above mentioned result of Berger & Bock.     

On Statistical Models for d-Dimensional Stable Processes, and Some Generalizations

In statistical models where jumps of a d -dimensional stable process ( S _t ) _{t ≥0} are observed in windows with certain asymptotic properties, and where parameters appearing in the Levy measure of S are to be estimated, we have asymptotically efficient estimators. If Poisson random measure μ on (0, ∞) × ( R ^d \{0}) with intensity dt Λ( dx ) replaces the jump measure of S , where Λ is a ε-finite measure on R ^d \{0} admitting tail parameters in a suitable sense, we specify a notion of neighbourhood which allows to treat efficiency in statistical experiments of the second type by switching to accompanying sequences of the stable process type considered first.     

Computation of the Generalized F Distribution

Exact expressions for the cumulative distribution function of a random variable of the form ( α ₁ X ₁+ α ₂ X ₂)/ Y are given where X ₁, X ₂ and Y are independent chi-squared random variables. The expressions are applied to the detection of joint outliers and Hotelling's mis-specified T ² distribution.     

Approximate Representation of Estimators in Constrained Regression Problems

The estimators of inequality-constrained regression problems can be computed by iterative algorithms of mathematical programming, but they do not have analytical expressions in terms of the given data. This situation brings obstacles to further studies on the constrained regression. In this paper we derive approximate representations of the estimators with a remainder of magnitude ( N ⁻¹ log log N )^1/2. From these representations one can clearly see the concrete structure of the estimators of these problems. It will be very helpful for further regression analysis.     

The NPMLE for Doubly Censored Current Status Data

In biostatistical applications interest often focuses on the estimation of the distribution of time T between two consecutive events. If the initial event time is observed and the subsequent event time is only known to be larger or smaller than an observed point in time, then the data is described by the well understood singly censored current status model, also known as interval censored data, case I. Jewell et al. (1994) extended this current status model by allowing the initial time to be unobserved, but with its distribution over an observed interval ' A, B ' known to be uniformly distributed; the data is referred to as doubly censored current status data. These authors used this model to handle application in AIDS partner studies focusing on the NPMLE of the distribution G of T . The model is a submodel of the current status model, but the distribution G is essentially the derivative of the distribution of interest F in the current status model. In this paper we establish that the NPMLE of G is uniformly consistent and that the resulting estimators for the n ^1/2-estimable parameters are efficient. We propose an iterative weighted pool-adjacent-violator-algorithm to compute the estimator. It is also shown that, without smoothness assumptions, the NPMLE of F converges at rate n ^−2/5 in L ²-norm while the NPMLE of F in the non-parametric current status data model converges at rate n ^−1/3 in L ²-norm, which shows that there is a substantial gain in using the submodel information.     

Likelihood Ratio Tests Under Local and Fixed Alternatives in Monotone Function Problems

Abstract.  We focus on a class of non-standard problems involving non-parametric estimation of a monotone function that is characterized by n ^1/3 rate of convergence of the maximum likelihood estimator, non-Gaussian limit distributions and the non-existence of     -regular estimators. We have shown elsewhere that under a null hypothesis of the type ψ ( z ₀) =  θ ₀ ( ψ being the monotone function of interest) in non-standard problems of the above kind, the likelihood ratio statistic has a 'universal' limit distribution that is free of the underlying parameters in the model. In this paper, we illustrate its limiting behaviour under local alternatives of the form ψ _n ( z ), where ψ _n (·) and ψ (·) vary in O ( n ^−1/3) neighbourhoods around z ₀ and ψ _n converges to ψ at rate n ^1/3 in an appropriate metric. Apart from local alternatives, we also consider the behaviour of the likelihood ratio statistic under fixed alternatives and establish the convergence in probability of an appropriately scaled version of the same to a constant involving a Kullback–Leibler distance.     

Parametric models for response-biased sampling

Suppose that subjects in a population follow the model f   ( y ^* x ^*; ) where y ^* denotes a response, x ^* denotes a vector of covariates and is the parameter to be estimated. We consider response-biased sampling, in which a subject is observed with a probability which is a function of its response. Such response-biased sampling frequently occurs in econometrics, epidemiology and survey sampling. The semiparametric maximum likelihood estimate of is derived, along with its asymptotic normality, efficiency and variance estimates. The estimate proposed can be used as a maximum partial likelihood estimate in stratified response-selective sampling. Some computation algorithms are also provided.     

On a Posterior Predictive Density Sample Size Criterion

Abstract.  Let Ω be a space of densities with respect to some σ -finite measure μ and let Π be a prior distribution having support Ω with respect to some suitable topology. Conditional on f , let X ⁿ  = ( X ₁ ,…, X _n ) be an independent and identically distributed sample of size n from f . This paper introduces a Bayesian non-parametric criterion for sample size determination which is based on the integrated squared distance between posterior predictive densities. An expression for the sample size is obtained when the prior is a Dirichlet mixture of normal densities.     

Optimal Bootstrap Sample Size in Construction of Percentile Confidence Bounds

In traditional bootstrap applications the size of a bootstrap sample equals the parent sample size, n say. Recent studies have shown that using a bootstrap sample size different from n may sometimes provide a more satisfactory solution. In this paper we apply the latter approach to correct for coverage error in construction of bootstrap confidence bounds. We show that the coverage error of a bootstrap percentile method confidence bound, which is of order O ( n ^−2/2) typically, can be reduced to O ( n ⁻¹) by use of an optimal bootstrap sample size. A simulation study is conducted to illustrate our findings, which also suggest that the new method yields intervals of shorter length and greater stability compared to competitors of similar coverage accuracy.     

Estimation of the Quantile Function of an IFRA Distribution

Let F and G be lifetime distributions and consider the problem of estimating F ⁻¹ when it is known that G ⁻¹ F is star-shaped. Estimators of F ⁻¹ are considered here which are shown to be uniformly strongly consistent. The case of censored data is also presented. Asymptotic confidence intervals and bands for F ⁻¹ are provided. The result are applicable, for example, to the estimation of quantile functions of k -out-of- n systems in reliability. The special case of an IFRA distribution follows immediately from the more general case presented here     

Testing a Regression Model When We Have Smooth Alternatives in Mind

Goodness-of-fit tests based on residual sums of squares are standard procedures used when fitting regression models. Often we have a smooth alternative in mind, a qualitative feature that the χ²-test does not take into account. We show that the power of detecting a smooth alternative increases when we smooth the current model as well. The proposed test is shown to be able to detect any continuous local alternative tending to zero slower than n ^−1/2. Theoretical results also address minimax non-parametric hypothesis testing in Sobolev spaces. A simulation study is presented, and the procedure is applied to expenditure curve estimation.     

A NOTE ON THE MULTIVARIATE LINEAR MODEL WITH CONSTRAINTS ON THE DEPENDENT VECTOR

It is shown that the least squares estimators of B and Σ in the multivariate linear model {E Y _i= X ₁ B , D ( Y _i) =Σ, 1 ≤ i ≤ n , Y ₁ Y _n uncorrelated} subject to the constraints Y _i M = X _i N are just the usual least squares estimators = ( X'X )-¹ X'Y and ΣC = 1/n( Y-X )( Y-X ) in the unconstrained model where Σ has full rank. Tests of hypotheses concerning B are discussed for situations in which each Y _i has a multivariate normal distribution, and examples of the applicability of the model reviewed.     

Bayesian model selection for join point regression with application to age-adjusted cancer rates

Summary.  The method of Bayesian model selection for join point regression models is developed. Given a set of K +1 join point models M ₀,  M ₁, …,  M _K with 0, 1, …,  K join points respec-tively, the posterior distributions of the parameters and competing models M _k are computed by Markov chain Monte Carlo simulations. The Bayes information criterion BIC is used to select the model M _k with the smallest value of BIC as the best model. Another approach based on the Bayes factor selects the model M _k with the largest posterior probability as the best model when the prior distribution of M _k is discrete uniform. Both methods are applied to analyse the observed US cancer incidence rates for some selected cancer sites. The graphs of the join point models fitted to the data are produced by using the methods proposed and compared with the method of Kim and co-workers that is based on a series of permutation tests. The analyses show that the Bayes factor is sensitive to the prior specification of the variance σ ², and that the model which is selected by BIC fits the data as well as the model that is selected by the permutation test and has the advantage of producing the posterior distribution for the join points. The Bayesian join point model and model selection method that are presented here will be integrated in the National Cancer Institute's join point software ( http://www.srab.cancer.gov/joinpoint/ ) and will be available to the public.     

Effectiveness of potent antiretroviral therapy on progression of human immunodeficiency virus: Bayesian modelling and model checking via counterfactual replicates

Summary.  We analyse data from a seroincident cohort of 457 homosexual men who were infected with the human immunodeficiency virus, followed within the multicentre Italian Seroconversion Study. These data include onset times to acquired immune deficiency syndrome (AIDS), longitudinal measurements of CD4⁺ T-cell counts taken on each subject during the AIDS-free period of observation and the period of administration of a highly active antiretro- viral therapy (HAART), for the subset of individuals who received it. The aim of the study is to assess the effect of HAART on the course of the disease. We analyse the data by a Bayesian model in which the sequence of longitudinal CD4⁺ cell count observations and the associated time to AIDS are jointly modelled at an individual subject's level as depending on the treatment. We discuss the inferences obtained about the efficacy of HAART, as well as modelling and computation difficulties that were encountered in the analysis. These latter motivate a model criticism stage of the analysis, in which the model specification of CD4⁺ cell count progression and of the effect of treatment are checked. Our approach to model criticism is based on the notion of a counterfactual replicate data set Z ^c . This is a data set with the same shape and size as the observed data, which we might have observed by rerunning the study in exactly the same conditions as the actual study if the treated patients had not been treated at all. We draw samples of Z ^c from a null model M ₀, which assumes absence of treatment effect, conditioning on data collected in each subject before initiation of treatment. Model checking is performed by comparing the observed data with a set of samples of Z ^c drawn from M ₀.     

A PROBABILISTIC DERIVATION OF THE NON-CENTRAL χ2 DISTRIBUTION

A simple derivation of the non-central χ² distribution is presented. This requires no advanced mathematical knowledge, and is suitable for use in elementary courses.
The non-central χ² distribution is of great importance in statistical theory, both in its own right, and as a step in the derivation of other distributions; however, it tends to be neglected in statistical courses largely, it seems, because the standard derivations are too difficult. A recent review of derivations by Guenther (1964) shows that most require some knowledge of n -dimensional geometry or the equivalent matrix theory. The alternative is the use of generating functions, which is straightforward, apart from the inversion back to the density function. An objection to such methods is that they give no insight into the probabilistic nature of the proble.     

Estimation of an Ergodic Diffusion from Discrete Observations

We consider a one-dimensional diffusion process X , with ergodic property, with drift b ( x , θ) and diffusion coefficient a ( x , σ) depending on unknown parameters θ and σ. We are interested in the joint estimation of (θ, σ). For that purpose, we dispose of a discretized trajectory, observed at n equidistant times tⁿ_i = ih_n , 1 ≤ i ≤ n . We assume that h_n ← 0 and nh_n ←∞. Under the condition nh_n^p ← 0 for an arbitrary integer p , we exhibit a contrast dependent on p which provides us with an asymptotically normal and efficient estimator of (θ, σ).