The spatial modeling of extreme snow is important for adequate risk management in Alpine and high altitude countries. A natural approach to such modeling is through the theory of max-stable processes, an infinite-dimensional extension of multivariate extreme value theory. In this paper we describe the application of such processes in modeling the spatial dependence of extreme snow depth in Switzerland, based on data for the winters 1966–2008 at 101 stations. The models we propose rely on a climate transformation that allows us to account for the presence of climate regions and for directional effects, resulting from synoptic weather patterns. Estimation is performed through pairwise likelihood inference and the models are compared using penalized likelihood criteria. The max-stable models provide a much better fit to the joint behavior of the extremes than do independence or full dependence models.

In observational studies the assignment of units to treatments is not under control. Consequently, the estimation and comparison of treatment effects based on the empirical distribution of the responses can be biased since the units exposed to the various treatments could differ in important unknown pretreatment characteristics, which are related to the response. An important example studied in this article is the question of whether private schools offer better quality of education than public schools. In order to address this question, we use data collected in the year 2000 by OECD for the Programme for International Student Assessment (PISA). Focusing for illustration on scores in mathematics of 15-year-old pupils in Ireland, we find that the raw average score of pupils in private schools is higher than of pupils in public schools. However, application of a newly proposed method for observational studies suggests that the less able pupils tend to enroll in public schools, such that their lower scores are not necessarily an indication of bad quality of the public schools. Indeed, when comparing the average score in the two types of schools after adjusting for the enrollment effects, we find quite surprisingly that public schools perform better on average. This outcome is supported by the methods of instrumental variables and latent variables, commonly used by econometricians for analyzing and evaluating social programs.

Reproducibility is essential to reliable scientific discovery in high-throughput experiments. In this work we propose a unified approach to measure the reproducibility of findings identified from replicate experiments and identify putative discoveries using reproducibility. Unlike the usual scalar measures of reproducibility, our approach creates a curve, which quantitatively assesses when the findings are no longer consistent across replicates. Our curve is fitted by a copula mixture model, from which we derive a quantitative reproducibility score, which we call the “irreproducible discovery rate” (IDR) analogous to the FDR. This score can be computed at each set of paired replicate ranks and permits the principled setting of thresholds both for assessing reproducibility and combining replicates.

Since our approach permits an arbitrary scale for each replicate, it provides useful descriptive measures in a wide variety of situations to be explored. We study the performance of the algorithm using simulations and give a heuristic analysis of its theoretical properties. We demonstrate the effectiveness of our method in a ChIP-seq experiment.

We consider applying Bayesian Variable Selection Regression, or BVSR, to genome-wide association studies and similar large-scale regression problems. Currently, typical genome-wide association studies measure hundreds of thousands, or millions, of genetic variants (SNPs), in thousands or tens of thousands of individuals, and attempt to identify regions harboring SNPs that affect some phenotype or outcome of interest. This goal can naturally be cast as a variable selection regression problem, with the SNPs as the covariates in the regression. Characteristic features of genome-wide association studies include the following: (i) a focus primarily on identifying relevant variables, rather than on prediction; and (ii) many relevant covariates may have tiny effects, making it effectively impossible to confidently identify the complete “correct” subset of variables. Taken together, these factors put a premium on having interpretable measures of confidence for individual covariates being included in the model, which we argue is a strength of BVSR compared with alternatives such as penalized regression methods. Here we focus primarily on analysis of quantitative phenotypes, and on appropriate prior specification for BVSR in this setting, emphasizing the idea of considering what the priors imply about the total proportion of variance in outcome explained by relevant covariates. We also emphasize the potential for BVSR to estimate this proportion of variance explained, and hence shed light on the issue of “missing heritability” in genome-wide association studies. More generally, we demonstrate that, despite the apparent computational challenges, BVSR can provide useful inferences in these large-scale problems, and in our simulations produces better power and predictive performance compared with standard single-SNP analyses and the penalized regression method LASSO. Methods described here are implemented in a software package, pi-MASS, available from the Guan Lab website http://bcm.edu/cnrc/mcmcmc/pimass.

Acute respiratory diseases are transmitted over networks of social contacts. Large-scale simulation models are used to predict epidemic dynamics and evaluate the impact of various interventions, but the contact behavior in these models is based on simplistic and strong assumptions which are not informed by survey data. These assumptions are also used for estimating transmission measures such as the basic reproductive number and secondary attack rates. Development of methodology to infer contact networks from survey data could improve these models and estimation methods. We contribute to this area by developing a model of within-household social contacts and using it to analyze the Belgian POLYMOD data set, which contains detailed diaries of social contacts in a 24-hour period. We model dependency in contact behavior through a latent variable indicating which household members are at home. We estimate age-specific probabilities of being at home and age-specific probabilities of contact conditional on two members being at home. Our results differ from the standard random mixing assumption. In addition, we find that the probability that all members contact each other on a given day is fairly low: 0.49 for households with two 0–5 year olds and two 19–35 year olds, and 0.36 for households with two 12–18 year olds and two 36+ year olds. We find higher contact rates in households with 2–3 members, helping explain the higher influenza secondary attack rates found in households of this size.

We consider the problem of algorithmically recommending items to users on a Yahoo! front page module. Our approach is based on a novel multilevel hierarchical model that we refer to as a User Profile Model with Graphical Lasso (UPG). The UPG provides a personalized recommendation to users by simultaneously incorporating both user covariates and historical user interactions with items in a model based way. In fact, we build a per-item regression model based on a rich set of user covariates and estimate individual user affinity to items by introducing a latent random vector for each user. The vector random effects are assumed to be drawn from a prior with a precision matrix that measures residual partial associations among items. To ensure better estimates of a precision matrix in high-dimensions, the matrix elements are constrained through a Lasso penalty. Our model is fitted through a penalized-quasi likelihood procedure coupled with a scalable EM algorithm. We employ several computational strategies like multi-threading, conjugate gradients and heavily exploit problem structure to scale our computations in the E-step. For the M-step we take recourse to a scalable variant of the Graphical Lasso algorithm for covariance selection.

Through extensive experiments on a new data set obtained from Yahoo! front page and a benchmark data set from a movie recommender application, we show that our UPG model significantly improves performance compared to several state-of-the-art methods in the literature, esp"otherecially those based on a bilinear random effects model (BIRE). In particular, we show that the gains of UPG are significant compared to BIRE when the number of users is large and the number of items to select from is small. For large item sets and relatively small user sets the results of UPG and BIRE are comparable. The UPG leads to faster model building and produces outputs which are interpretable.

Motivated by a potential-outcomes perspective, the idea of principal stratification has been widely recognized for its relevance in settings susceptible to posttreatment selection bias such as randomized clinical trials where treatment received can differ from treatment assigned. In one such setting, we address subtleties involved in inference for causal effects when using a key covariate to predict membership in latent principal strata. We show that when treatment received can differ from treatment assigned in both study arms, incorporating a stratum-predictive covariate can make estimates of the “complier average causal effect” (CACE) derive from observations in the two treatment arms with different covariate distributions. Adopting a Bayesian perspective and using Markov chain Monte Carlo for computation, we develop posterior checks that characterize the extent to which incorporating the pretreatment covariate endangers estimation of the CACE. We apply the method to analyze a clinical trial comparing two treatments for jaw fractures in which the study protocol allowed surgeons to overrule both possible randomized treatment assignments based on their clinical judgment and the data contained a key covariate (injury severity) predictive of treatment received.

Monte Carlo approaches have recently been proposed to quantify connectivity in neuronal networks. The key problem is to sample from the conditional distribution of a single neuronal spike train, given the activity of the other neurons in the network. Dependencies between neurons are usually relatively weak; however, temporal dependencies within the spike train of a single neuron are typically strong. In this paper we develop several specialized Metropolis–Hastings samplers which take advantage of this dependency structure. These samplers are based on two ideas: (1) an adaptation of fast forward–backward algorithms from the theory of hidden Markov models to take advantage of the local dependencies inherent in spike trains, and (2) a first-order expansion of the conditional likelihood which allows for efficient exact sampling in the limit of weak coupling between neurons. We also demonstrate that these samplers can effectively incorporate side information, in particular, noisy fluorescence observations in the context of calcium-sensitive imaging experiments. We quantify the efficiency of these samplers in a variety of simulated experiments in which the network parameters are closely matched to data measured in real cortical networks, and also demonstrate the sampler applied to real calcium imaging data.

A spatial stochastic model is developed which describes the 3D nanomorphology of composite materials, being blends of two different (organic and inorganic) solid phases. Such materials are used, for example, in photoactive layers of hybrid polymer zinc oxide solar cells. The model is based on ideas from stochastic geometry and spatial statistics. Its parameters are fitted to image data gained by electron tomography (ET), where adaptive thresholding and stochastic segmentation have been used to represent morphological features of the considered ET data by unions of overlapping spheres. Their midpoints are modeled by a stack of 2D point processes with a suitably chosen correlation structure, whereas a moving-average procedure is used to add the radii of spheres. The model is validated by comparing physically relevant characteristics of real and simulated data, like the efficiency of exciton quenching, which is important for the generation of charges and their transport toward the electrodes.

Landscape classification of the well-known biodiversity hotspot, Western Ghats (mountains), on the west coast of India, is an important part of a world-wide program of monitoring biodiversity. To this end, a massive vegetation data set, consisting of 51,834 4-variate observations has been clustered into different landscapes by Nagendra and Gadgil [Current Sci. 75 (1998) 264–271]. But a study of such importance may be affected by nonuniqueness of cluster analysis and the lack of methods for quantifying uncertainty of the clusterings obtained.

Motivated by this applied problem of much scientific importance, we propose a new methodology for obtaining the global, as well as the local modes of the posterior distribution of clustering, along with the desired credible and “highest posterior density” regions in a nonparametric Bayesian framework. To meet the need of an appropriate metric for computing the distance between any two clusterings, we adopt and provide a much simpler, but accurate modification of the metric proposed in [In Felicitation Volume in Honour of Prof. B. K. Kale (2009) MacMillan]. A very fast and efficient Bayesian methodology, based on [Sankhyā Ser. B 70 (2008) 133–155], has been utilized to solve the computational problems associated with the massive data and to obtain samples from the posterior distribution of clustering on which our proposed methods of summarization are illustrated.

Clustering of the Western Ghats data using our methods yielded landscape types different from those obtained previously, and provided interesting insights concerning the differences between the results obtained by Nagendra and Gadgil [Current Sci. 75 (1998) 264–271] and us. Statistical implications of the differences are also discussed in detail, providing interesting insights into methodological concerns of the traditional clustering methods.

The vast amount of biological knowledge accumulated over the years has allowed researchers to identify various biochemical interactions and define different families of pathways. There is an increased interest in identifying pathways and pathway elements involved in particular biological processes. Drug discovery efforts, for example, are focused on identifying biomarkers as well as pathways related to a disease. We propose a Bayesian model that addresses this question by incorporating information on pathways and gene networks in the analysis of DNA microarray data. Such information is used to define pathway summaries, specify prior distributions, and structure the MCMC moves to fit the model. We illustrate the method with an application to gene expression data with censored survival outcomes. In addition to identifying markers that would have been missed otherwise and improving prediction accuracy, the integration of existing biological knowledge into the analysis provides a better understanding of underlying molecular processes.

In biomedical studies it is of substantial interest to develop risk prediction scores using high-dimensional data such as gene expression data for clinical endpoints that are subject to censoring. In the presence of well-established clinical risk factors, investigators often prefer a procedure that also adjusts for these clinical variables. While accelerated failure time (AFT) models are a useful tool for the analysis of censored outcome data, it assumes that covariate effects on the logarithm of time-to-event are linear, which is often unrealistic in practice. We propose to build risk prediction scores through regularized rank estimation in partly linear AFT models, where high-dimensional data such as gene expression data are modeled linearly and important clinical variables are modeled nonlinearly using penalized regression splines. We show through simulation studies that our model has better operating characteristics compared to several existing models. In particular, we show that there is a nonnegligible effect on prediction as well as feature selection when nonlinear clinical effects are misspecified as linear. This work is motivated by a recent prostate cancer study, where investigators collected gene expression data along with established prognostic clinical variables and the primary endpoint is time to prostate cancer recurrence. We analyzed the prostate cancer data and evaluated prediction performance of several models based on the extended c statistic for censored data, showing that (1) the relationship between the clinical variable, prostate specific antigen, and the prostate cancer recurrence is likely nonlinear, that is, the time to recurrence decreases as PSA increases and it starts to level off when PSA becomes greater than 11; (2) correct specification of this nonlinear effect improves performance in prediction and feature selection; and (3) addition of gene expression data does not seem to further improve the performance of the resultant risk prediction scores.

Research in examining the equity of service accessibility has emerged as economic and social equity advocates recognized that where people live influences their opportunities for economic development, access to quality health care and political participation. In this research paper service accessibility equity is concerned with where and when services have been and are accessed by different groups of people, identified by location or underlying socioeconomic variables. Using new statistical methods for modeling spatial-temporal data, this paper estimates demographic association patterns to financial service accessibility varying over a large geographic area (Georgia) and over a period of 13 years. The underlying model is a space–time varying coefficient model including both separable space and time varying coefficients and space–time interaction terms. The model is extended to a multilevel response where the varying coefficients account for both the within- and between-variability. We introduce an inference procedure for assessing the shape of the varying regression coefficients using confidence bands.

Interactions among multiple genes across the genome may contribute to the risks of many complex human diseases. Whole-genome single nucleotide polymorphisms (SNPs) data collected for many thousands of SNP markers from thousands of individuals under the case–control design promise to shed light on our understanding of such interactions. However, nearby SNPs are highly correlated due to linkage disequilibrium (LD) and the number of possible interactions is too large for exhaustive evaluation. We propose a novel Bayesian method for simultaneously partitioning SNPs into LD-blocks and selecting SNPs within blocks that are associated with the disease, either individually or interactively with other SNPs. When applied to homogeneous population data, the method gives posterior probabilities for LD-block boundaries, which not only result in accurate block partitions of SNPs, but also provide measures of partition uncertainty. When applied to case–control data for association mapping, the method implicitly filters out SNP associations created merely by LD with disease loci within the same blocks. Simulation study showed that this approach is more powerful in detecting multi-locus associations than other methods we tested, including one of ours. When applied to the WTCCC type 1 diabetes data, the method identified many previously known T1D associated genes, including PTPN22, CTLA4, MHC, and IL2RA. The method also revealed some interesting two-way associations that are undetected by single SNP methods. Most of the significant associations are located within the MHC region. Our analysis showed that the MHC SNPs form long-distance joint associations over several known recombination hotspots. By controlling the haplotypes of the MHC class II region, we identified additional associations in both MHC class I (HLA-A, HLA-B) and class III regions (BAT1). We also observed significant interactions between genes PRSS16, ZNF184 in the extended MHC region and the MHC class II genes. The proposed method can be broadly applied to the classification problem with correlated discrete covariates.

We propose a semiparametric model for autonomous nonlinear dynamical systems and devise an estimation procedure for model fitting. This model incorporates subject-specific effects and can be viewed as a nonlinear semiparametric mixed effects model. We also propose a computationally efficient model selection procedure. We show by simulation studies that the proposed estimation as well as model selection procedures can efficiently handle sparse and noisy measurements. Finally, we apply the proposed method to a plant growth data used to study growth displacement rates within meristems of maize roots under two different experimental conditions.

Genetic association study is an essential step to discover genetic factors that are associated with a complex trait of interest. In this paper we present a novel generalized quasi-likelihood score (GQLS) test that is suitable for a study with either a quantitative trait or a binary trait. We use a logistic regression model to link the phenotypic value of the trait to the distribution of allelic frequencies. In our model, the allele frequencies are treated as a response and the trait is treated as a covariate that allows us to leave the distribution of the trait values unspecified. Simulation studies indicate that our method is generally more powerful in comparison with the family-based association test (FBAT) and controls the type I error at the desired levels. We apply our method to analyze data on Holstein cattle for an estimated breeding value phenotype, and to analyze data from the Collaborative Study of the Genetics of Alcoholism for alcohol dependence. The results show a good portion of significant SNPs and regions consistent with previous reports in the literature, and also reveal new significant SNPs and regions that are associated with the complex trait of interest.

In order to find previously unknown subgroups in biomedical data and generate testable hypotheses, visually guided exploratory analysis can be of tremendous importance. In this paper we propose a new dissimilarity measure that can be used within the Multidimensional Scaling framework to obtain a joint low-dimensional representation of both the samples and variables of a multivariate data set, thereby providing an alternative to conventional biplots. In comparison with biplots, the representations obtained by our approach are particularly useful for exploratory analysis of data sets where there are small groups of variables sharing unusually high or low values for a small group of samples.

Joinpoint regression is used to determine the number of segments needed to adequately explain the relationship between two variables. This methodology can be widely applied to real problems, but we focus on epidemiological data, the main goal being to uncover changes in the mortality time trend of a specific disease under study. Traditionally, Joinpoint regression problems have paid little or no attention to the quantification of uncertainty in the estimation of the number of change-points. In this context, we found a satisfactory way to handle the problem in the Bayesian methodology. Nevertheless, this novel approach involves significant difficulties (both theoretical and practical) since it implicitly entails a model selection (or testing) problem. In this study we face these challenges through (i) a novel reparameterization of the model, (ii) a conscientious definition of the prior distributions used and (iii) an encompassing approach which allows the use of MCMC simulation-based techniques to derive the results. The resulting methodology is flexible enough to make it possible to consider mortality counts (for epidemiological applications) as Poisson variables. The methodology is applied to the study of annual breast cancer mortality during the period 1980–2007 in Castellón, a province in Spain.

Drug discovery is the process of identifying compounds which have potentially meaningful biological activity. A major challenge that arises is that the number of compounds to search over can be quite large, sometimes numbering in the millions, making experimental testing intractable. For this reason computational methods are employed to filter out those compounds which do not exhibit strong biological activity. This filtering step, also called virtual screening reduces the search space, allowing for the remaining compounds to be experimentally tested.

In this paper we propose several novel approaches to the problem of virtual screening based on Canonical Correlation Analysis (CCA) and on a kernel-based extension. Spectral learning ideas motivate our proposed new method called Indefinite Kernel CCA (IKCCA). We show the strong performance of this approach both for a toy problem as well as using real world data with dramatic improvements in predictive accuracy of virtual screening over an existing methodology.

Originating from a system theory and an input/output point of view, I introduce a new class of generalized distributions. A parametric nonlinear transformation converts a random variable X into a so-called Lambert W random variable Y, which allows a very flexible approach to model skewed data. Its shape depends on the shape of X and a skewness parameter γ. In particular, for symmetric X and nonzero γ the output Y is skewed. Its distribution and density function are particular variants of their input counterparts. Maximum likelihood and method of moments estimators are presented, and simulations show that in the symmetric case additional estimation of γ does not affect the quality of other parameter estimates. Applications in finance and biomedicine show the relevance of this class of distributions, which is particularly useful for slightly skewed data. A practical by-result of the Lambert W framework: data can be “unskewed.”

The R package LambertW developed by the author is publicly available (CRAN).

GRavitational lEnsing Accuracy Testing 2010 (GREAT10) is a public image analysis challenge aimed at the development of algorithms to analyze astronomical images. Specifically, the challenge is to measure varying image distortions in the presence of a variable convolution kernel, pixelization and noise. This is the second in a series of challenges set to the astronomy, computer science and statistics communities, providing a structured environment in which methods can be improved and tested in preparation for planned astronomical surveys. GREAT10 extends upon previous work by introducing variable fields into the challenge. The “Galaxy Challenge” involves the precise measurement of galaxy shape distortions, quantified locally by two parameters called shear, in the presence of a known convolution kernel. Crucially, the convolution kernel and the simulated gravitational lensing shape distortion both now vary as a function of position within the images, as is the case for real data. In addition, we introduce the “Star Challenge” that concerns the reconstruction of a variable convolution kernel, similar to that in a typical astronomical observation. This document details the GREAT10 Challenge for potential participants. Continually updated information is also available from www.greatchallenges.info.