Combination chemotherapy treatment regimens created for patients diagnosed with childhood acute lymphoblastic leukemia have had great success in improving cure rates. Unfortunately, patients prescribed these types of treatment regimens have displayed susceptibility to the onset of osteonecrosis. Some have suggested that this is due to pharmacokinetic interaction between two agents in the treatment regimen (asparaginase and dexamethasone) and other physiological variables. Determining which physiological variables to consider when searching for interactions in scenarios like these, minus a priori guidance, has proved to be a challenging problem, particularly if interactions influence the response distribution in ways beyond shifts in expectation or dispersion only. In this paper we propose an exploratory technique that is able to discover associations between covariates and responses in a general way. The procedure connects covariates to responses flexibly through dependent random partition distributions and then employs machine learning techniques to highlight potential associations found in each cluster. We provide a simulation study to show utility and apply the method to data produced from a study dedicated to learning which physiological predictors influence severity of osteonecrosis multiplicatively.

Accurate estimates of historical changes in sea surface temperatures (SSTs) and their uncertainties are important for documenting and understanding historical changes in climate. A source of uncertainty that has not previously been quantified in historical SST estimates stems from position errors. A Bayesian inference framework is proposed for quantifying errors in reported positions and their implications on SST estimates. The analysis framework is applied to data from the International Comprehensive Ocean-Atmosphere Data Set (ICOADS3.0) in 1885, a time when astronomical and chronometer estimation of position was common but predated the use of radio signals. Focus is upon a subset of 943 ship tracks from ICOADS3.0 that report their position every two hours to a precision of ${0.01}^{\circ }$ longitude and latitude. These data are interpreted as positions determined by dead reckoning that are periodically updated by celestial correction techniques. The posterior medians of uncertainties in celestial correction are 33.1 km (${0.30}^{\circ }$ on the equator) in longitude and 24.4 km (${0.22}^{\circ }$) in latitude, respectively. Celestial navigation uncertainties being smaller in latitude than longitude is qualitatively consistent with the relative difficulty of obtaining astronomical estimates. The posterior medians for two-hourly dead reckoning uncertainties are $19.2\mathrm{\%}$ for ship speed and ${13.2}^{\circ }$ for ship heading, leading to random position uncertainties with median ${0.18}^{\circ }$ (20 km on the equator) in longitude and ${0.15}^{\circ }$ (17 km) in latitude. Reported ship tracks also contain systematic position uncertainties relating to precursor dead-reckoning positions not being updated after obtaining celestial position estimates, indicating that more accurate positions can be provided for SST observations. Finally, we translate position errors into SST uncertainties by sampling an ensemble of SSTs from the Multiscale Ultrahigh Resolution Sea Surface Temperature (MURSST) data set. Evolving technology for determining ship position, heterogeneous reporting and archiving of position information, and seasonal and spatial changes in navigational uncertainty and SST gradients together imply that accounting for positional error in SST estimates over the span of the instrumental record will require substantial additional effort.

Conventional analysis of neuroscience data involves computing average neural activity over a group of trials and/or a period of time. This approach may be particularly problematic when assessing the response patterns of neurons to more than one simultaneously presented stimulus. In such cases the brain must represent each individual component of the stimuli bundle, but trial-and-time-pooled averaging methods are fundamentally unequipped to address the means by which multiitem representation occurs. We introduce and investigate a novel statistical analysis framework that relates the firing pattern of a single cell, exposed to a stimuli bundle, to the ensemble of its firing patterns under each constituent stimulus. Existing statistical tools focus on what may be called “first order stochasticity” in trial-to-trial variation in the form of unstructured noise around a fixed firing rate curve associated with a given stimulus. Our analysis is based upon the theoretical premise that exposure to a stimuli bundle induces additional stochasticity in the cell’s response pattern in the form of a stochastically varying recombination of its single stimulus firing rate curves. We discuss challenges to statistical estimation of such “second order stochasticity” and address them with a novel dynamic admixture point process (DAPP) model. DAPP is a hierarchical point process model that decomposes second order stochasticity into a Gaussian stochastic process and a random vector of interpretable features and facilitates borrowing of information on the latter across repeated trials through latent clustering. We illustrate the utility and accuracy of the DAPP analysis with synthetic data simulation studies. We present real-world evidence of second order stochastic variation with an analysis of monkey inferior colliculus recordings under auditory stimuli.

The biomarker networks measured by different modalities of data (e.g., structural magnetic resonance imaging (sMRI), diffusion tensor imaging (DTI)) may share the same true underlying biological model. In this work we propose a nodewise biomarker graphical model to leverage the shared mechanism between multimodality data to provide a more reliable estimation of the target modality network and account for the heterogeneity in networks due to differences between subjects and networks of external modality. Latent variables are introduced to represent the shared unobserved biological network, and the information from the external modality is incorporated to model the distribution of the underlying biological network. We propose an efficient approximation to the posterior expectation of the latent variables that reduces computational cost by at least 50%. The performance of the proposed method is demonstrated by extensive simulation studies and an application to construct gray matter brain atrophy network of Huntington’s disease by using sMRI data and DTI data. The identified network connections are more consistent with clinical literature and better improve prediction in follow-up clinical outcomes and separate subjects into clinically meaningful subgroups with different prognosis than alternative methods.

Opioid overdose rates have increased in the United States over the past decade and reflect a major public health crisis. Modeling and prediction of drug and opioid hotspots, where a high percentage of events fall in a small percentage of space–time, could help better focus limited social and health services. In this work we present a spatial-temporal point process model for drug overdose clustering. The data input into the model comes from two heterogeneous sources: (1) high volume emergency medical calls for service (EMS) records containing location and time but no information on the type of nonfatal overdose, and (2) fatal overdose toxicology reports from the coroner containing location and high-dimensional information from the toxicology screen on the drugs present at the time of death. We first use nonnegative matrix factorization to cluster toxicology reports into drug overdose categories, and we then develop an EM algorithm for integrating the two heterogeneous data sets, where the mark corresponding to overdose category is inferred for the EMS data and the high volume EMS data is used to more accurately predict drug overdose death hotspots. We apply the algorithm to drug overdose data from Indianapolis, showing that the point process defined on the integrated data out-performs point processes that use only coroner data (AUC improvement 0.81 to 0.85). We also investigate the extent to which overdoses are contagious, as a function of the type of overdose, while controlling for exogenous fluctuations in the background rate that might also contribute to clustering. We find that drug and opioid overdose deaths exhibit significant excitation with branching ratio ranging from 0.72 to 0.98.

The National Resources Inventory, a longitudinal survey of characteristics related to natural resources and agriculture on nonfederal U.S. land, has increasingly received requests for substate estimates in recent years. We consider estimation of erosion in subdomains of the Boone-Raccoon River Watershed. This region is of interest for its proximity to intensively cropped areas as well as important waterbodies. The NRI application requires a small area prediction approach that can handle nonlinear relationships and appropriately incorporate survey weights that may have nontrivial relationships to the response variable. Because of the informative design, the conditional distribution required to define a standard empirical Bayes predictor is unknown. We develop a prediction approach that utilizes the approximate distribution of survey weighted score equations arising from a specified two-level superpopulation model. We apply the method to construct estimates of mean erosion in small watersheds. We investigate the robustness of the procedure to an assumption of a constant dispersion parameter and validate the properties of the procedure through simulation.

While electronic health records present a rich and promising data source for observational research, they are highly susceptible to missing data. For settings like these, Seaman et al. (Biometrics 68 (2012) 129–137) proposed a strategy wherein one handles missingness in some variables using inverse-probability weighting and others using multiple imputation. Seaman et al. (Biometrics 68 (2012) 129–137) show that Rubin’s variance estimator for averaging results across datasets is asymptotically valid when the analysis and imputation models are correctly specified and the weights are either known or correctly specified. Modeled after the approach of Robins and Wang (Biometrika 87 (2000) 113–124), we propose a method for asymptotically valid inference that is robust to violation of these conditions. Following a simulation study in which we demonstrate that a proposed variance estimator can reduce bias due to model misspecification, we illustrate this approach in an electronic health records-based study investigating whether differences in long-term weight loss between bariatric surgery techniques are associated with chronic kidney disease at baseline. We observe that the weight loss advantage after five years of Roux-en-Y gastric bypass surgery, compared to vertical sleeve gastrectomy, is less pronounced among patients with chronic kidney disease at baseline compared to those without.

Biclustering methods simultaneously group samples and their associated features. In this way, biclustering methods differ from traditional clustering methods, which utilize the entire set of features to distinguish groups of samples. Motivating applications for biclustering include genomics data, where the goal is to cluster patients or samples by their gene expression profiles; and recommender systems, which seek to group customers based on their product preferences. Biclusters of interest often manifest as rank-1 submatrices of the data matrix. This submatrix detection problem can be viewed as a factor analysis problem in which both the factors and loadings are sparse. In this paper, we propose a new biclustering method called Spike-and-Slab Lasso Biclustering (SSLB) which utilizes the Spike-and-Slab Lasso of Ročková and George (J. Amer. Statist. Assoc. 113 (2018) 431–444) to find such a sparse factorization of the data matrix. SSLB also incorporates an Indian Buffet Process prior to automatically choose the number of biclusters. Many biclustering methods make assumptions about the size of the latent biclusters; either assuming that the biclusters are all of the same size, or that the biclusters are very large or very small. In contrast, SSLB can adapt to find biclusters which have a continuum of sizes. SSLB is implemented via a fast EM algorithm with a variational step. In a variety of simulation settings, SSLB outperforms other biclustering methods. We apply SSLB to both a microarray dataset and a single-cell RNA-sequencing dataset and highlight that SSLB can recover biologically meaningful structures in the data. The SSLB software is available as an R/C++ package at https://github.com/gemoran/SSLB.

Phase I dose-finding studies aim at identifying the maximum tolerated dose (MTD). Often, several dose-finding studies are conducted with some variation in the administration mode or dose panel. For instance, sorafenib (BAY 43-900) was used as monotherapy in 36 phase I trials, according to a recent clinicaltrials.gov search. Since the toxicity may not be directly related to the specific indication, synthesizing the information from several studies might be worthwhile. However, this is rarely done in practice and only a fixed-effect meta-analysis framework was proposed to date. We developed a Bayesian random-effects meta-analysis methodology to pool several phase I trials and suggest the MTD. A curve free hierarchical model on the logistic scale with random effects, accounting for between-trial heterogeneity, is used to model the probability of toxicity across the investigated doses. An Ornstein–Uhlenbeck Gaussian process is adopted for the random effects structure. Prior distributions for the curve-free model are based on a latent Gamma process. An extensive simulation study showed good performance of the proposed method also under model deviations. Sharing information between phase I studies can improve the precision of MTD selection, at least when the number of trials is reasonably large.

Graphical model is a powerful and popular approach to study high-dimensional omic data, such as genome-wide gene expression data. Nonlinear relations between genes are widely documented. However, partly due to sparsity of data points in high-dimensional space (i.e., curse of dimensionality) and computational challenges, most available methods construct graphical models by testing linear relations. We propose to address this challenge by a two-step approach: first, use a model-free approach to prioritize the neighborhood of each gene; then, apply a nonparametric conditional independence testing method to refine such neighborhood estimation. Our method, named as “mofreds” (MOdel FRee Estimation of DAG Skeletons), seeks to estimate the skeleton of a directed acyclic graph (DAG) by this two-step approach. We studied the theoretical properties of mofreds and evaluated its performance in extensive simulation settings. We found mofreds has substantially better performance than the state-of-the art method which is designed to detect linear relations of Gaussian graphical models. We applied mofreds to analyze gene expression data of breast cancer patients from The Cancer Genome Atlas (TCGA). We found that it discovers nonlinear relationships among gene pairs that are missed by the Gaussian graphical model methods.

Several programs of research have sought to assess the replicability of scientific findings in different fields, including economics and psychology. These programs attempt to replicate several findings and use the results to say something about large-scale patterns of replicability in a field. However, little work has been done to understand the analytic methods used to do this, including what they are assessing and what their statistical properties are. This article examines several methods that have been used to study patterns of replicability in the social sciences. We describe in concrete terms how each method operationalizes the idea of “replication” and examine various statistical properties, including bias, precision and statistical power. We find that some analytic methods rely on an operational definition of replication that can be misleading. Other methods involve more sound definitions of replication, but most of these have limitations, such as large bias and uncertainty or low power. The findings suggest that we should use caution interpreting the results of such analyses and that work on more accurate methods may be useful to future replication research efforts.

Inferring concerted changes among biological traits along an evolutionary history remains an important yet challenging problem. Besides adjusting for spurious correlation induced from the shared history, the task also requires sufficient flexibility and computational efficiency to incorporate multiple continuous and discrete traits as data size increases. To accomplish this, we jointly model mixed-type traits by assuming latent parameters for binary outcome dimensions at the tips of an unknown tree informed by molecular sequences. This gives rise to a phylogenetic multivariate probit model. With large sample sizes, posterior computation under this model is problematic, as it requires repeated sampling from a high-dimensional truncated normal distribution. Current best practices employ multiple-try rejection sampling that suffers from slow-mixing and a computational cost that scales quadratically in sample size. We develop a new inference approach that exploits: (1) the bouncy particle sampler (BPS) based on piecewise deterministic Markov processes to simultaneously sample all truncated normal dimensions, and (2) novel dynamic programming that reduces the cost of likelihood and gradient evaluations for BPS to linear in sample size. In an application with 535 HIV viruses and 24 traits that necessitates sampling from a 12,840-dimensional truncated normal, our method makes it possible to estimate the across-trait correlation and detect factors that affect the pathogen’s capacity to cause disease. This inference framework is also applicable to a broader class of covariance structures beyond comparative biology.

The Vaccine Adverse Event Reporting System (VAERS) plays a vital role in vaccine safety surveillance. One of the main missions of VAERS is to monitor increases in reporting rate of adverse events, as such signals can indicate safety issues caused by update of vaccines or change in vaccine practices. Existing methods can rarely be used to monitor the temporal variation of reporting adverse events. In this paper we propose a composite likelihood based variance component model to study the temporal variation of reporting adverse events using VAERS data. The proposed method is devised to identify safety signals by testing the heterogeneity of reporting rates of adverse events across years. The proposed method accounts for the well-known underreporting of adverse events and the zero-inflation observations in passive surveillance reporting systems. We applied the proposed method to VAERS reports of trivalent influenza virus vaccine and identified 14 adverse events with significantly heterogeneous reporting rates over years and two of them have increasing trend of reporting rates from 1990 to 2013. Our findings provide early warning signals that can be more rigorously investigated in future studies of the vaccine.

Complex tissues are composed of a large number of different types of cells, each involved in a multitude of biological processes. Consequently, an important component to understanding such processes is understanding the cell-type composition of the tissues. Estimating cell-type composition using high-throughput gene expression data is known as cell-type deconvolution. In this paper we first summarize the extensive deconvolution literature by identifying a common regression-like approach to deconvolution. We call this approach the unified deconvolution-as-regression (UDAR) framework. While methods that fall under this framework all use a similar model, they fit using data on different scales. Two popular scales for gene expression data are logarithmic and linear. Unfortunately, each of these scales has problems in the UDAR framework. Using log-scale gene expressions proposes a biologically implausible model and using linear-scale gene expressions will lead to statistically inefficient estimators. To explore ways to address these issues, in this paper we consider how deconvolution methods may use an adjusted model that is a hybrid of the two scales. In analysis on simulations as well as a collection of eleven real benchmark datasets, we find a prototypical hybrid-scale adjustment to the UDAR framework improves statistical efficiency and robustness. More broadly, we believe these hybrid-scale modeling principles may be incorporated into many existing deconvolution methods.

The estimation of a surface throughout a continuum of points in a study area is addressed by means of two-phase sampling strategies. To this aim, a family of two-phase inverse distance weighting interpolators is introduced, and their design-based asymptotic properties are derived when the surface remains fixed and the number of sample points approaches infinity. In particular, conditions ensuring asymptotic unbiasedness and consistency are derived and are proven to hold for some of the most widely applied environmental sampling schemes. Furthermore, a computationally simple mean squared error estimator is proposed. Finally, a simulation study is performed to assess the theoretical results. The proposed strategy is adopted to provide the map of basal area in a forested region of Casentino Valley (Central Italy).

Regression discontinuity (RD) is a widely used quasi-experimental design for causal inference. In the standard RD the assignment to treatment is determined by a continuous pretreatment variable (i.e., running variable) falling above or below a prefixed threshold. Recent applications increasingly feature ordered categorical or ordinal running variables which pose challenges to RD estimation due to the lack of a meaningful measure of distance. This paper proposes an RD approach for ordinal running variables under the local randomization framework. The proposal first estimates an ordered probit model for the ordinal running variable. The estimated probability of being assigned to treatment is then adopted as a latent continuous running variable and used to identify a covariate-balanced subsample around the threshold. Assuming local unconfoundedness of the treatment in the subsample, an estimate of the effect of the program is obtained by employing a weighted estimator of the average treatment effect. Two weighting estimators—overlap weights and ATT weights—as well as their augmented versions are considered. We apply the method to evaluate the causal effects of the corporate sector purchase programme (CSPP) of the European Central Bank which involves large-scale purchases of securities issued by corporations in the euro area. We find a statistically significant and negative effect of the CSPP on corporate bond spreads at issuance.

We introduce spatial (DLfuse) and spatiotemporal (DLfuseST) distributed lag data fusion methods for predicting point-level ambient air pollution concentrations, using, as input, gridded average pollution estimates from a deterministic numerical air quality model. The methods incorporate predictive information from grid cells surrounding the prediction location of interest and are shown to collapse to existing downscaling approaches when this information adds no benefit. The spatial lagged parameters are allowed to vary spatially/spatiotemporally to accommodate the setting where surrounding geographic information is useful in one area/time but not in another. We apply the new methods to predict ambient concentrations of eight-hour maximum ozone and 24-hour average ${\text{PM}}_{2.5}$ at unobserved spatial locations and times, and compare the predictions with those from several state-of-the-art data fusion approaches. Results show that DLfuse and DLfuseST often provide improved model fit and predictive accuracy when the lagged information is shown to be beneficial. Code to apply the methods is available in the R package $\mathtt{DLfuse}$.

CRISPR technology has enabled cell lineage tracing for complex multicellular organisms through insertion-deletion mutations of synthetic genomic barcodes during organismal development. To reconstruct the cell lineage tree from the mutated barcodes, current approaches apply general-purpose computational tools that are agnostic to the mutation process and are unable to take full advantage of the data’s structure. We propose a statistical model for the CRISPR mutation process and develop a procedure to estimate the resulting tree topology, branch lengths and mutation parameters by iteratively applying penalized maximum likelihood estimation. By assuming the barcode evolves according to a molecular clock, our method infers relative ordering across parallel lineages, whereas existing techniques only infer ordering for nodes along the same lineage. When analyzing transgenic zebrafish data from (Science 353 (2016) aaf7907), we find that our method recapitulates known aspects of zebrafish development and the results are consistent across samples.

The coronavirus disease 2019 (COVID-19) has quickly grown from a regional outbreak in Wuhan, China, to a global pandemic. Early estimates of the epidemic growth and incubation period of COVID-19 may have been biased due to sample selection. Using detailed case reports from 14 locations in and outside mainland China, we obtained 378 Wuhan-exported cases who left Wuhan before an abrupt travel quarantine. We developed a generative model we call BETS for four key epidemiological events—Beginning of exposure, End of exposure, time of Transmission, and time of Symptom onset (BETS)—and derived explicit formulas to correct for the sample selection. We gave a detailed illustration of why some early and highly influential analyses of the COVID-19 pandemic were severely biased. All our analyses, regardless of which subsample and model were being used, point to an epidemic doubling time of $\mathit{two}$ to 2.5 days during the early outbreak in Wuhan. A Bayesian nonparametric analysis further suggests that about 5% of the symptomatic cases may not develop symptoms within 14 days of infection and that men may be much more likely than women to develop symptoms within two days of infection.

Gene expression deconvolution is a powerful tool for exploring the microenvironment of complex tissues comprised of multiple cell groups using transcriptomic data. Characterizing cell activities for a particular condition has been regarded as a primary mission against diseases. For example, cancer immunology aims to clarify the role of the immune system in the progression and development of cancer through analyzing the immune cell components of tumors. To that end, many deconvolution methods have been proposed for inferring cell subpopulations within tissues. Nevertheless, two problems limit the practicality of current approaches. First, most approaches use external purified data to preselect cell type-specific genes that contribute to deconvolution. However, some types of cells cannot be found in purified profiles, and the genes specifically over- or under-expressed in them cannot be identified. This is particularly a problem in cancer studies. Hence, a preselection strategy that is independent from deconvolution is inappropriate. The second problem is that existing approaches do not recover the expression profiles of unknown cells present in bulk tissues when the reference set of purified cell-specific profiles is incomplete which results in biased estimation of unknown cell proportions. Furthermore, it causes the shift-invariant property of deconvolution to fail which then affects the estimation performance. To address these two problems, we propose a novel semireference-based deconvolution approach, BayICE which employs hierarchical Bayesian modeling with stochastic search variable selection. We develop a comprehensive Markov chain Monte Carlo procedure through Gibbs sampling to estimate proportions, expression profiles and signature genes for a set of known reference cell types as well as an unknown cell type. Simulation and validation studies illustrate that BayICE outperforms existing semireference-based deconvolution approaches in estimating cell proportions. We further show that BayICE is applicable to single-cell RNA-seq data. Subsequently, we demonstrate an application of BayICE in the RNA sequencing of patients with nonsmall cell lung cancer. The model is implemented in the R package “BayICE,” and the algorithm is available for download.

Causal analysis of observational studies requires data that comprise a set of covariates, a treatment assignment indicator and the observed outcomes. However, data confidentiality restrictions or the nature of data collection may distribute these variables across two or more datasets. In the absence of unique identifiers to link records across files, probabilistic record linkage algorithms can be leveraged to merge the datasets. Current applications of record linkage are concerned with estimation of associations between variables that are exclusive to one file and not causal relationships. We propose a Bayesian framework for record linkage and causal inference where one file comprises all the covariate and observed outcome information, and the second file consists of a list of all individuals who receive the active treatment. Under certain ignorability assumptions, the procedure properly propagates the error in the record linkage process, resulting in valid statistical inferences. To estimate the causal effects, we devise a two-stage procedure. The first stage of the procedure performs Bayesian record linkage to multiply-impute the treatment assignment for all individuals in the first file, while adjustments for covariates’ imbalance and imputation of missing potential outcomes are performed in the second stage. This procedure is used to evaluate the effect of Meals on Wheels services on mortality and healthcare utilization among homebound older adults in Rhode Island. In addition, an interpretable sensitivity analysis is developed to assess potential violations of the ignorability assumptions.

Smoking is one of the main risk factors that has affected human mortality and life expectancy over the past century. Smoking accounts for a large part of the nonlinearities in the growth of life expectancy and of the geographic and gender differences in mortality. As Bongaarts (Popul. Dev. Rev. 32 (2006) 605–628) and Janssen (Genus 74 (2018) 21) suggested, accounting for smoking could improve the quality of mortality forecasts due to the predictable nature of the smoking epidemic. We propose a new Bayesian hierarchical model to forecast life expectancy at birth for both genders and for 69 countries/regions with good data on smoking-related mortality. The main idea is to convert the forecast of the nonsmoking life expectancy at birth (i.e., life expectancy at birth removing the smoking effect) into life expectancy forecast through the use of the age-specific smoking attributable fraction (ASSAF). We introduce a new age-cohort model for the ASSAF and a Bayesian hierarchical model for nonsmoking life expectancy at birth. The forecast performance of the proposed method is evaluated by out-of-sample validation compared with four other commonly used methods for life expectancy forecasting. Improvements in forecast accuracy and model calibration based on the new method are observed.

Probabilistic hurricane storm surge forecasting using a high-fidelity model has been considered impractical due to the overwhelming computational expense to run thousands of simulations. This article demonstrates that modern statistical tools enable good forecasting performance using a small number of carefully chosen simulations. This article offers algorithms that quickly handle the massive output of a surge model while addressing the missing data at unsubmerged locations. Also included is a new optimal design criterion for selecting simulations that accounts for the log transform required to statistically model surge data. Hurricane Michael (2018) is used as a testbed for this investigation and provides evidence for the approach’s efficacy in comparison to the existing probabilistic surge forecast method.

Clustering is a form of unsupervised learning that aims to uncover latent groups within data based on similarity across a set of features. A common application of this in biomedical research is in delineating novel cancer subtypes from patient gene expression data, given a set of informative genes. However, it is typically unknown a priori what genes may be informative in discriminating between clusters and what the optimal number of clusters are. Few methods exist for performing unsupervised clustering of RNA-seq samples, and none currently adjust for between-sample global normalization factors, select cluster-discriminatory genes or account for potential confounding variables during clustering. To address these issues, we propose the feature selection and clustering of RNA-seq (FSCseq): a model-based clustering algorithm that utilizes a finite mixture of regression (FMR) model and the quadratic penalty method with a smoothly clipped absolute deviation (SCAD) penalty. The maximization is done by a penalized Classification EM algorithm, allowing us to include normalization factors and confounders in our modeling framework. Given the fitted model, our framework allows for subtype prediction in new patients via posterior probabilities of cluster membership, even in the presence of batch effects. Based on simulations and real data analysis, we show the advantages of our method relative to competing approaches.

Time-course experiments are commonly conducted to capture temporal changes. It is generally of interest to detect if any changes happen over time, which we define as a detection problem. If there is a change, it is informative to know when the change is, which we define as an identification problem. It is often desired to control Type I error rate at a nominal level while applying a testing procedure to detect or identify these changes. Quite a few analytic methods have been proposed. Most existing methods aim to solve either the detection problem or, more recently, the identification problem. Here, we propose to solve these two problems using a unified multiple-testing framework built upon an empirical Bayes change-point model. Our model provides a flexible framework that can account for sophisticated temporal gene expression patterns. We show that our testing procedure is valid and asymptotically optimal in the sense of rejecting the maximum number of null hypotheses, while the Bayesian false discovery rate (FDR) can be controlled at a predefined nominal level. Simulation studies and application to real transcriptome time-course data illustrate that our proposed model is a flexible and powerful method to capture various temporal patterns in analysis of time-course data.