Statistical significance analysis of longitudinal gene expression data

MOTIVATION: Time-course microarray experiments are designed to study biological processes in a temporal fashion. Longitudinal gene expression data arise when biological samples taken from the same subject at different time points are used to measure the gene expression levels. It has been observed that the gene expression patterns of samples of a given tumor measured at different time points are likely to be much more similar to each other than are the expression patterns of tumor samples of the same type taken from different subjects. In statistics, this phenomenon is called the within-subject correlation of repeated measurements on the same subject, and the resulting data are called longitudinal data. It is well known in other applications that valid statistical analyses have to appropriately take account of the possible within-subject correlation in longitudinal data. RESULTS: We apply estimating equation techniques to construct a robust statistic, which is a variant of the robust Wald statistic and accounts for the potential within-subject correlation of longitudinal gene expression data, to detect genes with temporal changes in expression. We associate significance levels to the proposed statistic by either incorporating the idea of the significance analysis of microarrays method or using the mixture model method to identify significant genes. The utility of the statistic is demonstrated by applying it to an important study of osteoblast lineage-specific differentiation. Using simulated data, we also show pitfalls in drawing statistical inference when the within-subject correlation in longitudinal gene expression data is ignored.     

Genome-wide linkage analysis of systolic blood pressure: a comparison of two approaches to phenotype definition

Problem 1 of the Genetic Analysis Workshop 13(GAW13) contains longitudinal data of cardiovascular measurements from 330 pedigrees. The longitudinal data complicates the phenotype definition because multiple measurements are taken on each individual. To address this complication, we propose an approach that uses generalized estimating equations to obtain residuals for each time point for each person. The mean residual is then taken as the new phenotype with which to use in a variance components linkage analysis. We compare our phenotype definition approach to an approach that first reduces the multiple measurements to a single measurement and then models these summary statistics as regression terms in a variance components analysis. For each approach, multipoint linkage analysis was performed using the residuals and the SOLAR computer program. Our results show little difference between the methods based on the LOD scores.     

On Efficiency of Constrained Longitudinal Data Analysis Versus Longitudinal Analysis of Covariance

Summary : In randomized clinical trials, measurements are often collected on each subject at a baseline visit and several post‐randomization time points. The longitudinal analysis of covariance in which the postbaseline values form the response vector and the baseline value is treated as a covariate can be used to evaluate the treatment differences at the postbaseline time points. Liang and Zeger (2000, Sankhyā: The Indian Journal of Statistics, Series B 62, 134–148) propose a constrained longitudinal data analysis in which the baseline value is included in the response vector together with the postbaseline values and a constraint of a common baseline mean across treatment groups is imposed on the model as a result of randomization. If the baseline value is subject to missingness, the constrained longitudinal data analysis is shown to be more efficient for estimating the treatment differences at postbaseline time points than the longitudinal analysis of covariance. The efficiency gain increases with the number of subjects missing baseline and the number of subjects missing all postbaseline values, and, for the pre–post design, decreases with the absolute correlation between baseline and postbaseline values.     

Multikink quantile regression for longitudinal data with application to progesterone data analysis

Motivated by investigating the relationship between progesterone and the days in a menstrual cycle in a longitudinal study, we propose a multikink quantile regression model for longitudinal data analysis. It relaxes the linearity condition and assumes different regression forms in different regions of the domain of the threshold covariate. In this paper, we first propose a multikink quantile regression for longitudinal data. Two estimation procedures are proposed to estimate the regression coefficients and the kink points locations: one is a computationally efficient profile estimator under the working independence framework while the other one considers the within-subject correlations by using the unbiased generalized estimation equation approach. The selection consistency of the number of kink points and the asymptotic normality of two proposed estimators are established. Second, we construct a rank score test based on partial subgradients for the existence of the kink effect in longitudinal studies. Both the null distribution and the local alternative distribution of the test statistic have been derived. Simulation studies show that the proposed methods have excellent finite sample performance. In the application to the longitudinal progesterone data, we identify two kink points in the progesterone curves over different quantiles and observe that the progesterone level remains stable before the day of ovulation, then increases quickly in 5 to 6 days after ovulation and then changes to stable again or drops slightly.     

Quantile regression for longitudinal data using the asymmetric Laplace distribution

In longitudinal studies, measurements of the same individuals are taken repeatedly through time. Often, the primary goal is to characterize the change in response over time and the factors that influence change. Factors can affect not only the location but also more generally the shape of the distribution of the response over time. To make inference about the shape of a population distribution, the widely popular mixed-effects regression, for example, would be inadequate, if the distribution is not approximately Gaussian. We propose a novel linear model for quantile regression (QR) that includes random effects in order to account for the dependence between serial observations on the same subject. The notion of QR is synonymous with robust analysis of the conditional distribution of the response variable. We present a likelihood-based approach to the estimation of the regression quantiles that uses the asymmetric Laplace density. In a simulation study, the proposed method had an advantage in terms of mean squared error of the QR estimator, when compared with the approach that considers penalized fixed effects. Following our strategy, a nearly optimal degree of shrinkage of the individual effects is automatically selected by the data and their likelihood. Also, our model appears to be a robust alternative to the mean regression with random effects when the location parameter of the conditional distribution of the response is of interest. We apply our model to a real data set which consists of self-reported amount of labor pain measurements taken on women repeatedly over time, whose distribution is characterized by skewness, and the significance of the parameters is evaluated by the likelihood ratio statistic.     

Morphometry and allometry of the primate humerus

The primate distal humerus has been used both in phylogenetic reconstruction and in assessing locomotor and postural adaptations. This study uses an allometric approach to predict locomotor patterns of extant primates regardless of phylogenetic position. By showing the relationship between form and function in living primate taxa it will be possible to use this data set to predict locomotor behavior of extinct primates. Several linear measurements were taken from the distal humerus of 71 extant primate species (anthropoids and prosimians). Allometric regressions of each measurement were performed with mandibular M₂ area as a surrogate for body size. These measurements were used to determine if significant differences in distal humerus morphology exist among locomotor groups. The results were then used to test several hypotheses about the relationship between humeral form and function. For example, the hypothesis that suspensory primates have a large medial epicondyle is confirmed; the hypothesis that terrestrial quadrupeds have a deep olecranon fossa could not be confirmed with quantitative data. In addition to this hypothesis testing, the residuals from the allometric regressions of the humeral measurements were used in a discriminant functions analysis to estimate locomotor behavior from distal humerus morphology. The discriminant functions analysis correctly reclassified 64/71 (90%) species.     

A generalized nonparametric test for lattice-ordered means

Treatment means in factorial experiments are lattice ordered when there is an increase in mean response as the level of any factor is increased while holding the other factors fixed. Such means occur naturally in many experiments. A nonparametric test for lattice-ordered means involving a Kendall-type statistic will be summarized for k-factor factorial experiments. Specifically, the form of the test statistic and variance under the null hypothesis will be presented. In addition, a normalized version of the test statistic will be discussed and applied to relevant data.     

A morphometric analysis of the fetal craniofacies by ultrasound: fetal cephalometry

We present here a set of 24 standardized linear measurements that describe the growth of different craniofacial structures in the normal fetus from 16 to 36 weeks of gestation. These measurements were taken from 89 pregnant women, who had from 1 to 3 ultrasonographic evaluations during the pregnancy (16, 26, and 36 weeks of gestation). All the values presented here were obtained using the technique described by Escobar et al. The mean and standard deviation was calculated for each measurement and was used to estimate the normal growth pattern of each variable. Approximate confidence intervals for the mean of each variable were constructed for use in identifying unusually low or high values. The confidence intervals are available in graphic form by request. These data will not only contribute to an understanding of fetal craniofacial growth and development in utero, but in addition, it will help to make the diagnoses of mild craniofacial anomalies that would not be detected by the routine ultrasonographic examination. We suggest that this procedure should be included if not in all routine obstetrical ultrasound evaluations, then at least in the more extensive level II obstetrical ultrasound.     

Bayesian approaches to joint cure-rate and longitudinal models with applications to cancer vaccine trials

Complex issues arise when investigating the association between longitudinal immunologic measures and time to an event, such as time to relapse, in cancer vaccine trials. Unlike many clinical trials, we may encounter patients who are cured and no longer susceptible to the time-to-event endpoint. If there are cured patients in the population, there is a plateau in the survival function, S(t), after sufficient follow-up. If we want to determine the association between the longitudinal measure and the time-to-event in the presence of cure, existing methods for jointly modeling longitudinal and survival data would be inappropriate, since they do not account for the plateau in the survival function. The nature of the longitudinal data in cancer vaccine trials is also unique, as many patients may not exhibit an immune response to vaccination at varying time points throughout the trial. We present a new joint model for longitudinal and survival data that accounts both for the possibility that a subject is cured and for the unique nature of the longitudinal data. An example is presented from a cancer vaccine clinical trial.     

Day-to-day variation in nitrogenase activity of Alnus incana explained by weather variables: a multivariate time series analysis

A modelling system is described that indicates the extent to which day-to-day variations in nitrogenase activity in young Alnus incana (L.) Moench, grown in defined conditions in the field, may be affected by weather conditions both during and prior to the day of measurement. Nitrogenase activity (acetylene reduction activity, ARA) was measured weekly on intact field-grown grey alder (A. incana) plants, 0.15–0.42 m tall at planting, nodulated with Frankia. The measurements were done at noon on two groups of plants in 1987 and on two other groups in 1988. Each group was made up of five or six plants. Seven weather variables: daily sunshine hours, daily mean, maximum and minimum air temperature, daily mean and 1300 h relative humidity, and daily rainfall were used. The relation between log(ARA/leaf area) and the weather variables were analysed using a PLS model (partial least squares projection to latent structures). The advantage of PLS is that it can handle x-variables that are correlated. Data from 1987 were chosen as a training set. Multivariate PLS time series analysis was made by adding, in a stepwise manner, the weather data up to 5 d before the day of measurement. This procedure gave six models with n * 7 x-variables (n= 1–6). With the models from the time series analysis of 1987 data, true predictions of ARA per leaf area were made from weather data 1988 (test set 1) and from ‘early-season’ weather data from 1987 and 1988 (test set 2). The variation in ARA/leaf area could be predicted from the weather conditions. The predictions of the two test sets improved when the weather conditions one and two days before the day of measurements were added to the model. The further addition of weather data from 3 to 5 d before the day of measurement did not improve the model. The good predictions of ARA/leaf area show that the alders responded to the variable weather conditions in the same way in 1988 as in 1987, despite the ten-fold difference in size (leaf area) at the end of the growing season. Among the weather variables, air temperature and the daily sunshine hours were positively correlated to ARA, while relative air humidity and rainfall were negatively correlated to ARA. The daily minimum temperature and rainfall appeared to have least impact on ARA. By use of PLS, we could extract information out of a data set containing highly correlated x-variables, information that is non-accessible with conventional statistical tools such as multiple regression. When making measurements of nitrogenase activities under field conditions, we propose that attention should be paid to the weather conditions on the days preceding the day of measurement. The day-to-day variation in nitrogenase activity is discussed with reference to known effects of stress factors under controlled conditions.     

Methods for the statistical analysis of binary data in split-cluster designs

Split-cluster designs are frequently used in the health sciences when naturally occurring clusters such as multiple sites or organs in the same subject are assigned to different treatments. However, statistical methods for the analysis of binary data arising from such designs are not well developed. The purpose of this article is to propose and evaluate a new procedure for testing the equality of event rates in a design dividing each of k clusters into two segments having multiple sites (e.g., teeth, lesions). The test statistic proposed is a generalization of a previously published procedure based on adjusting the standard Pearson chi-square statistic, but can also be derived as a score test using the approach of generalized estimating equations.     

Joint analysis of longitudinal data with informative right censoring

Longitudinal data arise when subjects are followed over a period of time. A commonly encountered complication in the analysis of such data is the variable length of follow-up due to right censorship. This can be further exacerbated by the possible dependency between the censoring time and the longitudinal measurements. This article proposes a combination of a semiparametric transformation model for the censoring time and a linear mixed effects model for the longitudinal measurements. The dependency is handled via latent variables which are naturally incorporated. We show that the likelihood function has an explicit form and develops a two-stage estimation procedure to avoid direct maximization over a high-dimensional parameter space. The resulting estimators are shown to be consistent and asymptotically normal, with a closed form for the variance-covariance matrix that can be used to obtain a plug-in estimator. Finite sample performance of the proposed approach is assessed through extensive simulations. The method is applied to renal disease data.     

Regression analysis of mean quality-adjusted lifetime with censored data

In clinical trials of chronic diseases such as acquired immunodeficiency syndrome, cancer, or cardiovascular diseases, the concept of quality-adjusted lifetime (QAL) has received more and more attention. In this paper, we consider the problem of how the covariates affect the mean QAL when the data are subject to right censoring. We allow a very general form for the mean model as a function of covariates. Using the idea of inverse probability weighting, we first construct a simple weighted estimating equation for the parameters in our mean model. We then find the form of the most efficient estimating equation, which yields the most efficient estimator for the regression parameters. Since the most efficient estimator depends on the distribution of the health history processes, and thus cannot be estimated nonparametrically, we consider different approaches for improving the efficiency of the simple weighted estimating equation using observed data. The applicability of these methods is demonstrated by both simulation experiments and a data example from a breast cancer clinical trial study.     

Time‐Varying Latent Effect Model for Longitudinal Data with Informative Observation Times

Summary In analysis of longitudinal data, it is not uncommon that observation times of repeated measurements are subject‐specific and correlated with underlying longitudinal outcomes. Taking account of the dependence between observation times and longitudinal outcomes is critical under these situations to assure the validity of statistical inference. In this article, we propose a flexible joint model for longitudinal data analysis in the presence of informative observation times. In particular, the new procedure considers the shared random‐effect model and assumes a time‐varying coefficient for the latent variable, allowing a flexible way of modeling longitudinal outcomes while adjusting their association with observation times. Estimating equations are developed for parameter estimation. We show that the resulting estimators are consistent and asymptotically normal, with variance–covariance matrix that has a closed form and can be consistently estimated by the usual plug‐in method. One additional advantage of the procedure is that it provides a unified framework to test whether the effect of the latent variable is zero, constant, or time‐varying. Simulation studies show that the proposed approach is appropriate for practical use. An application to a bladder cancer data is also given to illustrate the methodology.     

A New Nonparametric Approach to the Comparison of K Independent Samples of Response Curves II: A K Sample Generalization of the FRIEDMAN Test

A K sample generalization of the FRIEDMAN test (1937) is introduced which can be used as a nonparametric procedure for testing the homogeneity of the profiles of K independent samples of response curves measured at T identical points of time. While a similar procedure in LEHMACHER & WALL (1978), section 3, is based on T combined tests, each of them at level a/T, here a finite and asymptotic test is presented which is based on a single test statistic. The application of the new multivariate test is illustrated by the same numerical example as in LEHMACHER & WALL (1978). The properties of this test are discussed and compared with the combined test mentioned above.     

Survival analysis of longitudinal microarrays

MOTIVATION: The development of methods for linking gene expressions to various clinical and phenotypic characteristics is an active area of genomic research. Scientists hope that such analysis may, for example, describe relationships between gene function and clinical events such as death or recovery. Methods are available for relating gene expression to measurements that are categorized or continuous, but there is less work in relating expressions to an observed event time such as time to death, response or relapse. When gene expressions are measured over time, there are methods for differentiating temporal patterns. However, methods have not yet been proposed for the survival analysis of longitudinally collected microarrays. RESULTS: We describe an approach for the survival analysis of longitudinal gene expression data. We construct a measure of association between the time to an event and gene expressions collected over time. Statistical significance is addressed using permutations and control of the false discovery rate. Our proposed method is illustrated on a dataset from a multi-center research study of inflammation and response to injury that aims to uncover the biological reasons why patients can have dramatically different outcomes after suffering a traumatic injury (www.gluegrant.org).     

The Wilcoxon signed rank test for paired comparisons of clustered data

The Wilcoxon signed rank test is a frequently used nonparametric test for paired data (e.g., consisting of pre- and posttreatment measurements) based on independent units of analysis. This test cannot be used for paired comparisons arising from clustered data (e.g., if paired comparisons are available for each of two eyes of an individual). To incorporate clustering, a generalization of the randomization test formulation for the signed rank test is proposed, where the unit of randomization is at the cluster level (e.g., person), while the individual paired units of analysis are at the subunit within cluster level (e.g., eye within person). An adjusted variance estimate of the signed rank test statistic is then derived, which can be used for either balanced (same number of subunits per cluster) or unbalanced (different number of subunits per cluster) data, with an exchangeable correlation structure, with or without tied values. The resulting test statistic is shown to be asymptotically normal as the number of clusters becomes large, if the cluster size is bounded. Simulation studies are performed based on simulating correlated ranked data from a signed log-normal distribution. These studies indicate appropriate type I error for data sets with > or =20 clusters and a superior power profile compared with either the ordinary signed rank test based on the average cluster difference score or the multivariate signed rank test of Puri and Sen. Finally, the methods are illustrated with two data sets, (i) an ophthalmologic data set involving a comparison of electroretinogram (ERG) data in retinitis pigmentosa (RP) patients before and after undergoing an experimental surgical procedure, and (ii) a nutritional data set based on a randomized prospective study of nutritional supplements in RP patients where vitamin E intake outside of study capsules is compared before and after randomization to monitor compliance with nutritional protocols.     

A Simple Two‐Sample Rank Test for Multivariate Survival Outcomes with Left Truncation and Right Censoring

A distribution‐free two‐sample rank test is proposed for testing for differences between survival distributions in the analysis of biomedical studies in which two groups of subjects are followed over time for a particular outcome, which may recur. This method is motivated by an observational HIV (human immunodeficiency virus) study in which a group of HIV‐seropositive women and a comparable group of HIV‐seronegative women were examined every 6 months for the presence of cervical intraepithelial neoplasia (CIN), the cervical cancer precursor. Women entered the study serially and were subject to random loss to follow‐up. Only women free of CIN at study entry were followed resulting in left‐truncated survival times. If a woman is found to be CIN infected at a later examination, she is treated and then followed until CIN recurs. The two groups of women were compared at both occurrences of CIN on the basis of rank statistics. For the first occurrence of CIN, survival times since the beginning of the study (based on calendar time) are compared. For a recurrence of CIN, survival times since the first development of CIN are compared. The proposed test statistic for an overall difference between the two groups follows a chi‐square distribution with two degrees of freedom. Simulation results demonstrate the usefulness of the proposed test proposed test statistic, which reduces to the Gehan statistic if each person is followed only to the first failure and there is no serial enrollment.     

From SNPs to genes: disease association at the gene level

Interpreting Genome-Wide Association Studies (GWAS) at a gene level is an important step towards understanding the molecular processes that lead to disease. In order to incorporate prior biological knowledge such as pathways and protein interactions in the analysis of GWAS data it is necessary to derive one measure of association for each gene. We compare three different methods to obtain gene-wide test statistics from Single Nucleotide Polymorphism (SNP) based association data: choosing the test statistic from the most significant SNP; the mean test statistics of all SNPs; and the mean of the top quartile of all test statistics. We demonstrate that the gene-wide test statistics can be controlled for the number of SNPs within each gene and show that all three methods perform considerably better than expected by chance at identifying genes with confirmed associations. By applying each method to GWAS data for Crohn's Disease and Type 1 Diabetes we identified new potential disease genes.     

Marginal analysis of longitudinal ordinal data with misclassification in both response and covariates

Marginal methods have been widely used for the analysis of longitudinal ordinal and categorical data. These models do not require full parametric assumptions on the joint distribution of repeated response measurements but only specify the marginal or even association structures. However, inference results obtained from these methods often incur serious bias when variables are subject to error. In this paper, we tackle the problem that misclassification exists in both response and categorical covariate variables. We develop a marginal method for misclassification adjustment, which utilizes second‐order estimating functions and a functional modeling approach, and can yield consistent estimates and valid inference for mean and association parameters. We propose a two‐stage estimation approach for cases in which validation data are available. Our simulation studies show good performance of the proposed method under a variety of settings. Although the proposed method is phrased to data with a longitudinal design, it also applies to correlated data arising from clustered and family studies, in which association parameters may be of scientific interest. The proposed method is applied to analyze a dataset from the Framingham Heart Study as an illustration.