Showing papers in "Journal of The Royal Statistical Society Series B-statistical Methodology in 1999"

Probabilistic Principal Component Analysis

Abstract: Principal component analysis (PCA) is a ubiquitous technique for data analysis and processing, but one which is not based upon a probability model. In this paper we demonstrate how the principal axes of a set of observed data vectors may be determined through maximum-likelihood estimation of parameters in a latent variable model closely related to factor analysis. We consider the properties of the associated likelihood function, giving an EM algorithm for estimating the principal subspace iteratively, and discuss the advantages conveyed by the definition of a probability density function for PCA.

Statistical applications of the multivariate skew normal distribution

Abstract: Azzalini and Dalla Valle have recently discussed the multivariate skew normal distribution which extends the class of normal distributions by the addition of a shape parameter. The first part of the present paper examines further probabilistic properties of the distribution, with special emphasis on aspects of statistical relevance. Inferential and other statistical issues are discussed in the following part, with applications to some multivariate statistics problems, illustrated by numerical examples. Finally, a further extension is described which introduces a skewing factor of an elliptical density.

Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm

Abstract: Summary. Two new implementations of the EM algorithm are proposed for maximum likelihood fitting of generalized linear mixed models. Both methods use random (independent and identically distributed) sampling to construct Monte Carlo approximations at the E-step. One approach involves generating random samples from the exact conditional distribution of the random effects (given the data) by rejection sampling, using the marginal distribution as a candidate. The second method uses a multivariate t importance sampling approximation. In many applications the two methods are complementary. Rejection sampling is more efficient when sample sizes are small, whereas importance sampling is better with larger sample sizes. Monte Carlo approximation using random samples allows the Monte Carlo error at each iteration to be assessed by using standard central limit theory combined with Taylor series methods. Specifically, we construct a sandwich variance estimate for the maximizer at each approximate E-step. This suggests a rule for automatically increasing the Monte Carlo sample size after iterations in which the true EM step is swamped by Monte Carlo error. In contrast, techniques for assessing Monte Carlo error have not been developed for use with alternative implementations of Monte Carlo EM algorithms utilizing Markov chain Monte Carlo E-step approximations. Three different data sets, including the infamous salamander data of McCullagh and Nelder, are used to illustrate the techniques and to compare them with the alternatives. The results show that the methods proposed can be considerably more efficient than those based on Markov chain Monte Carlo algorithms. However, the methods proposed may break down when the intractable integrals in the likelihood function are of high dimension.

Inference in generalized additive mixed modelsby using smoothing splines

Abstract: Generalized additive mixed models are proposed for overdispersed and correlated data, which arise frequently in studies involving clustered, hierarchical and spatial designs. This class of models allows flexible functional dependence of an outcome variable on covariates by using nonparametric regression, while accounting for correlation between observations by using random effects. We estimate nonparametric functions by using smoothing splines and jointly estimate smoothing parameters and variance components by using marginal quasi-likelihood. Because numerical integration is often required by maximizing the objective functions, double penalized quasi-likelihood is proposed to make approximate inference. Frequentist and Bayesian inferences are compared. A key feature of the method proposed is that it allows us to make systematic inference on all model components within a unified parametric mixed model framework and can be easily implemented by fitting a working generalized linear mixed model by using existing statistical software. A bias correction procedure is also proposed to improve the performance of double penalized quasi-likelihood for sparse data. We illustrate the method with an application to infectious disease data and we evaluate its performance through simulation.

Gibbs sampling for Bayesian non‐conjugate and hierarchical models by using auxiliary variables

Abstract: Summary. We demonstrate the use of auxiliary (or latent) variables for sampling non-standard densities which arise in the context of the Bayesian analysis of non-conjugate and hierarchical models by using a Gibbs sampler. Their strategic use can result in a Gibbs sampler having easily sampled full conditionals. We propose such a procedure to simplify or speed up the Markov chain Monte Carlo algorithm. The strength of this approach lies in its generality and its ease of implementation. The aim of the paper, therefore, is to provide an alternative sampling algorithm to rejection-based methods and other sampling approaches such as the Metropolis-Hastings algorithm.

Missing covariates in generalized linear models when the missing data mechanism is non-ignorable

Abstract: We propose a method for estimating parameters in generalized linear models with missing covariates and a non-ignorable missing data mechanism. We use a multinomial model for the missing data indicators and propose a joint distribution for them which can be written as a sequence of one-dimensional conditional distributions, with each one-dimensional conditional distribution consisting of a logistic regression. We allow the covariates to be either categorical or continuous. The joint covariate distribution is also modelled via a sequence of one-dimensional conditional distributions, and the response variable is assumed to be completely observed. We derive the E- and M-steps of the EM algorithm with non-ignorable missing covariate data. For categorical covariates, we derive a closed form expression for the E- and M-steps of the EM algorithm for obtaining the maximum likelihood estimates (MLEs). For continuous covariates, we use a Monte Carlo version of the EM algorithm to obtain the MLEs via the Gibbs sampler. Computational techniques for Gibbs sampling are proposed and implemented. The parametric form of the assumed missing data mechanism itself is not `testable' from the data, and thus the non-ignorable modelling considered here can be viewed as a sensitivity analysis concerning a more complicated model. Therefore, although a model may have `passed' the tests for a certain missing data mechanism, this does not mean that we have captured, even approximately, the correct missing data mechanism. Hence, model checking for the missing data mechanism and sensitivity analyses play an important role in this problem and are discussed in detail. Several simulations are given to demonstrate the methodology. In addition, a real data set from a melanoma cancer clinical trial is presented to illustrate the methods proposed.

Bayesian Nonparametric Inference for Random Distributions and Related Functions

Abstract: In recent years, Bayesian nonparametric inference, both theoretical and computational, has witnessed considerable advances. However, these advances have not received a full critical and comparative analysis of their scope, impact and limitations in statistical modelling; many aspects of the theory and methods remain a mystery to practitioners and many open guestions remain. In this paper, we discuss and illustrate the rich modelling and analytic possibilities that are available to the statistician within the Bayesian nonparametric and/or semiparametric framework.

Direct calculation of the information matrix via the EM

Abstract: A simple explicit formula is provided for the matrix of second derivatives of the observed data log-likelihood in terms of derivatives of the criterion function (conditional expectation of the complete data log-likelihood given the observed data) invoked by the EM algorithm.

Semiparametric methods for response-selective and missing data problems in regression

Abstract: Summary. Suppose that data are generated according to the model f(ylx; 0) g(x), where y is a response and x are covariates. We derive and compare semiparametric likelihood and pseudolikelihood methods for estimating 0 for situations in which units generated are not fully observed and in which it is impossible or undesirable to model the covariate distribution. The probability that a unit is fully observed may depend on y, and there may be a subset of covariates which is observed only for a subsample of individuals. Our key assumptions are that the probability that a unit has missing data depends only on which of a finite number of strata that (y, x) belongs to and that the stratum membership is observed for every unit. Applications include case-control studies in epidemiology, field reliability studies and broad classes of missing data and measurement error problems. Our results make fully efficient estimation of 0 feasible, and they generalize and provide insight into a variety of methods that have been proposed for specific problems.

Bayesian analysis of agricultural field experiments

Abstract: The paper describes Bayesian analysis for agricultural field experiments, a topic that has received very little previous attention, despite a vast frequentist literature. Adoption of the Bayesian paradigm simplifies the interpretation of the results, especially in ranking and selection. Also, complex formulations can be analysed with comparative ease, by using Markov chain Monte Carlo methods. A key ingredient in the approach is the need for spatial representations of the unobserved fertility patterns. This is discussed in detail. Problems caused by outliers and by jumps in fertility are tackled via hierarchical-f formulations that may find use in other contexts. The paper includes three analyses of variety trials for yield and one example involving binary data; none is entirely straightforward. Some numerical comparisons with frequentist analyses are made.

On the time to extinction in recurrent epidemics

Abstract: Summary. An approximation is derived for the expected time to extinction in a stochastic model for recurrent epidemics. Numerical illustrations indicate that the approximation is crude but that it has the correct order of magnitude. The quasi-stationary distribution plays an important role in the derivation. Approximations for the critical community size and for the persistence threshold are derived. Comments are made on the classical study by Bartlett (1956-1960).

Conformal normal curvature and assessment of local influence

Abstract: In 1986, R. D. Cook proposed differential geometry to assess local influence of minor perturbations of statistical models. We construct a conformally invariant curvature, the conformal normal curvature, for the same purpose. This curvature provides a measure of local influence ranging from 0 to 1, with objective bench-marks to judge largeness. We study various approaches to using the conformal normal curvature and the relationships between these approaches.

Fractal analysis of surface roughness by using spatial data

Abstract: Summary. We develop fractal methodology for data taking the form of surfaces. An advantage of fractal analysis is that it partitions roughness characteristics of a surface into a scale-free component (fractal dimension) and properties that depend purely on scale. Particular emphasis is given to anisotropy where we show that, for many surfaces, the fractal dimension of line transects across a surface must either be constant in every direction or be constant in each direction except one. This virtual direction invariance of fractal dimension provides another canonical feature of fractal analysis, complementing its scale invariance properties and enhancing its attractiveness as a method for summarizing properties of roughness. The dependence of roughness on direction may be explained in terms of scale rather than dimension and can vary with orientation. Scale may be described by a smooth periodic function and may be estimated nonparametrically. Our results and techniques are applied to analyse data on the surfaces of soil and plastic food wrapping. For the soil data, interest centres on the effect of surface roughness on retention of rain-water, and data are recorded as a series of digital images over time. Our analysis captures the way in which both the fractal dimension and the scale change with rainfall, or equivalently with time. The food wrapping data are on a much finer scale than the soil data and are particularly anisotropic. The analysis allows us to determine the manufacturing process which produces the smoothest wrapping, with least tendency for micro-organisms to adhere.

Convergence of slice sampler Markov chains

Abstract: We analyse theoretical properties of the slice sampler. We find that the algorithm has extremely robust geometric ergodicity properties. For the case of just one auxiliary variable, we demonstrate that the algorithm is stochastically monotone, and we deduce analytic bounds on the total variation distance from stationarity of the method by using Foster–Lyapunov drift condition methodology.

Intentionally biased bootstrap methods

Abstract: Summary. A class of weighted bootstrap techniques, called biased bootstrap or b-bootstrap methods, is introduced. It is motivated by the need to adjust empirical methods, such as the ‘uniform’ bootstrap, in a surgical way to alter some of their features while leaving others unchanged. Depending on the nature of the adjustment, the b-bootstrap can be used to reduce bias, or to reduce variance or to render some characteristic equal to a predetermined quantity. Examples of the last application include a b-bootstrap approach to hypothesis testing in nonparametric contexts, where the b-bootstrap enables simulation ‘under the null hypothesis’, even when the hypothesis is false, and a b-bootstrap competitor to Tibshirani’s variance stabilization method. An example of the bias reduction application is adjustment of Nadaraya‐Watson kernel estimators to make them competitive with local linear smoothing. Other applications include density estimation under constraints, outlier trimming, sensitivity analysis, skewness or kurtosis reduction and shrinkage.

Minimum aberration and model robustness for two-level fractional factorial designs

Abstract: Summary. The performance of minimum aberration two-level fractional factorial designs is studied under two criteria of model robustness. Simple sufficient conditions for a design to dominate another design with respect to each of these two criteria are derived. It is also shown that a minimum aberration design of resolution IlIl or higher maximizes the number of two-factor interactions which are not aliases of main effects and, subject to that condition, minimizes the sum of squares of the sizes of alias sets of two-factor interactions. This roughly says that minimum aberration designs tend to make the sizes of the alias sets very uniform. It follows that minimum aberration is a good surrogate for the two criteria of model robustness that are studied here. Examples are given to show that minimum aberration designs are indeed highly efficient.

Statistical studies of infectious disease incidence

Abstract: Summary. Methods for the analysis of data on the incidence of an infectious disease are reviewed, with an emphasis on important objectives that such analyses should address and identifying areas where further work is required. Recent statistical work has adapted methods for constructing estimating functions from martingale theory, methods of data augmentation and methods developed for studying the human immunodeficiency virus-acquired immune deficiency syndrome epidemic. Infectious disease data seem particularly suited to analysis by Markov chain Monte Carlo methods. Epidemic modellers have recently made substantial progress in allowing for community structure and heterogeneity among individuals when studying the requirements for preventing major epidemics. This has stimulated interest in making statistical inferences about crucial parameters from infectious disease data for such community settings.

Effect of frailty on marginal regression estimates in survival analysis

Abstract: Unexplained heterogeneity in univariate survival data and association in multivariate survival can both be modelled by the inclusion of frailty effects. This paper investigates the consequences of ignoring frailty in analysis, fitting misspecified Cox proportional hazards models to the marginal distributions. Regression coefficients are biased towards 0 by an amount which depends in magnitude on the variability of the frailty terms and the form of frailty distribution. The bias is reduced when censoring is present. Fitted marginal survival curves can also differ substantially from the true marginals.

The Covariance Inflation Criterion for Adaptive Model Selection

Abstract: We propose a new criterion for model selection in prediction problems. The covariance inflation criterion adjusts the training error by the average covariance of the predictions and responses, when the prediction rule is applied to permuted versions of the data set. This criterion can be applied to general prediction problems (e.g. regression or classification) and to general prediction rules (e.g. stepwise regression, tree-based models and neural nets). As a by-product we obtain a measure of the effective number of parameters used by an adaptive procedure. We relate the covariance inflation criterion to other model selection procedures and illustrate its use in some regression and classification problems. We also revisit the conditional bootstrap approach to model selection.

Prior elicitation, variable selection and Bayesian computation for logistic regression models

Abstract: Summary. Bayesian selection of variables is often difficult to carry out because of the challenge in specifying prior distributions for the regression parameters for all possible models, specifying a prior distribution on the model space and computations. We address these three issues for the logistic regression model. For the first, we propose an informative prior distribution for variable selection. Several theoretical and computational properties of the prior are derived and illustrated with several examples. For the second, we propose a method for specifying an informative prior on the model space, and for the third we propose novel methods for computing the marginal distribution of the data. The new computational algorithms only require Gibbs samples from the full model to facilitate the computation of the prior and posterior model probabilities for all possible models. Several properties of the algorithms are also derived. The prior specification for the first challenge focuses on the observables in that the elicitation is based on a prior prediction yo for the response vector and a quantity ao quantifying the uncertainty in yo. Then, yo and ao are used to specify a prior for the regression coefficients semi-automatically. Examples using real data are given to demonstrate the methodology.

Covariance structure of wavelet coefficients: theory and models in a Bayesian perspective

Abstract: We present theoretical results on the random wavelet coefficients covariance structure. We use simple properties of the coefficients to derive a recursive way to compute the within- and across-scale covariances. We point out a usetul link between the algorithm proposed and the two-dimensional discrete wavelel transform. We then focus on Bayesian wavelet shrinkage for estimating a function from noisy data A prior distribution is imposed on the coefficients of the unknown function. We show how our findings on the covariance structure make it possible to specify priors that take into account the full correlation between coefficients through a parsimonious number of hyperparameters. We use Markov chain Monte Carlo rnethods to estirnate the parameters and illustrate our method on bench-mark simulated signals.

Multivariate bandwidth selection for local linear regression

Abstract: The existence and properties of optimal bandwidths for multivariate local linear regression are established, using either a scalar bandwidth for all regressors or a diagonal bandwidth vector that has a different bandwidth for each regressor. Both involve functionals of the derivatives of the unknown multivariate regression function. Estimating these functionals is difficult primarily because they contain multivariate derivatives. In this paper, an estimator of the multivariate second derivative is obtained via local cubic regression with most cross-terms left out. This estimator has the optimal rate of convergence but is simpler and uses much less computing time than the full local estimator. Using this as a pilot estimator, we obtain plug-in formulae for the optimal bandwidth, both scalar and diagonal, for multivariate local linear regression. As a simpler alternative, we also provide rule-of-thumb bandwidth selectors. All these bandwidths have satisfactory performance in our simulation study.

Perfect simulation of conditionally specified models

Abstract: We discuss how the ideas of producing perfect simulations based on coupling from the past for finite state space models naturally extend to multivariate distributions with infinite or uncountable state spaces such as auto-gamma, auto-Poisson and autonegative binomial models, using Gibbs sampling in combination with sandwiching methods originally introduced for perfect simulation of point processes.

One‐step local quasi‐likelihood estimation

Abstract: Summary. Local quasi-likelihood estimation is a useful extension of local least squares methods, but its computational cost and algorithmic convergence problems make the procedure less appealing, particularly when it is iteratively used in methods such as the back-fitting algorithm, cross-validation and bootstrapping. A one-step local quasi-likelihood estimator is introduced to overcome the computational drawbacks of the local quasi-likelihood method. We demonstrate that as long as the initial estimators are reasonably good the one-step estimator has the same asymptotic behaviour as the local quasi-likelihood method. Our simulation shows that the one-step estimator performs at least as well as the local quasi-likelihood method for a wide range of choices of bandwidths. A data-driven bandwidth selector is proposed for the one-step estimator based on the pre-asymptotic substitution method of Fan and Gijbels. It is then demonstrated that the datadriven one-step local quasi-likelihood estimator performs as well as the maximum local quasilikelihood estimator by using the ideal optimal bandwidth.

Priors and component structures in autoregressive time series models

Abstract: New approaches to prior specification and structuring in autoregressive time series models are introduced and developed. We focus on defining classes of prior distributions for parameters and latent variables related to latent components of an autoregressive model for an observed time series. These new priors naturally permit the incorporation of both qualitative and quantitative prior information about the number and relative importance of physically meaningful components that represent low frequency trends, quasi-periodic subprocesses and high frequency residual noise components of observed series. The class of priors also naturally incorporates uncertainty about model order and hence leads in posterior analysis to model order assessment and resulting posterior and predictive inferences that incorporate full uncertainties about model order as well as model parameters. Analysis also formally incorporates uncertainty and leads to inferences about unknown initial values of the time series, as it does for predictions of future values. Posterior analysis involves easily implemented iterative simulation methods, developed and described here. One motivating field of application is climatology, where the evaluation of latent structure, especially quasi-periodic structure, is of critical importance in connection with issues of global climatic variability. We explore the analysis of data from the southern oscillation index, one of several series that has been central in recent high profile debates in the atmospheric sciences about recent apparent trends in climatic indicators.

Density and hazard rate estimation for right‐censored data by using wavelet methods

Abstract: This paper describes a wavelet method for the estimation of density and hazard rate functions from randomly right-censored data. We adopt a nonparametric approach in assuming that the density and hazard rate have no specific parametric form. The method is based on dividing the time axis into a dyadic number of intervals and then counting the number of events within each interval. The number of events and the survival function of the observations are then separately smoothed over time via linear wavelet smoothers, and then the hazard rate function estimators are obtained by taking the ratio. We prove that the estimators have pointwise and global mean-square consistency, obtain the best possible asymptotic mean integrated squared error convergence rate and are also asymptotically normally distributed. We also describe simulation experiments that show that these estimators are reasonably reliable in practice. The method is illustrated with two real examples. The first uses survival time data for patients with liver metastases from a colorectal primary tumour without other distant metastases. The second is concerned with times of unemployment for women and the wavelet estimate, through its flexibility, provides a new and interesting interpretation.

Weighted empirical adaptive variance estimators for correlated data regression

Abstract: Estimating equations based on marginal generalized linear models are useful for regression modelling of correlated data, but inference and testing require reliable estimates of standard errors. We introduce a class of variance estimators based on the weighted empirical variance of the estimating functions and show that an adaptive choice of weights allows reliable estimation both asymptotically and by simulation in finite samples. Connections with previous bootstrap and jackknife methods are explored. The effect of reliable variance estimation is illustrated in data on health effects of air pollution in King County, Washington.

Understanding exponential smoothing via kernel regression

Abstract: Exponential smoothing is the most common model-free means of forecasting a future realization of a time series. It requires the specification of a smoothing factor which is usually chosen from the data to minimize the average squared residual of previous one-step-ahead forecasts. In this paper we show that exponential smoothing can be put into a nonparametric regression framework and gain some interesting insights into its performance through this interpretation. We also use theoretical developments from the kernel regression field to derive, for the first time, asymptotic properties of exponential smoothing forecasters.

Analysing competing risks data with transformation models

Abstract: We present a flexible class of marginal models for the cumulative incidence function. The semiparametric transformation model is utilized in a decomposition for the marginal failure probabilities which extends previous work on Farewell's cure model. Novel estimation, inference and prediction procedures are developed, with large sample properties derived from the theory of martingales and U-statistics. A small simulation study demonstrates that the methods are appropriate for practical use. The methods are illustrated with a thorough analysis of a prostate cancer clinical trial. Simple graphical displays are used to check for the goodness of fit.

The achievable region approach to the optimal control of stochastic systems

Abstract: Summary. The achievable region approach seeks solutions to stochastic optimization problems by characterizing the space of all possible performances (the achievable region) of the system of interest and optimizing the overall system-wide performance objective over this space. This is radically different from conventional formulations based on dynamic programming. The approach is explained with reference to a simple two-class queuing system. Powerful new methodologies due to the authors and co-workers are deployed to analyse a general multiclass queuing system with parallel servers and then to develop an approach to optimal load distribution across a network of interconnected stations. Finally, the approach is used for the first time to analyse a class of intensity control problems.