In non-proportional hazards models, the hazard ratio for a unit increase in covariate value is not constant but varies over time. Existing approaches to estimating time-varying log hazard ratio include various spline approximations and maximum penalised partial likelihood. We consider improvements to these methods under the plausible assumption that the hazard ratio changes from its short-term to long-term value in a monotonic fashion. A monotone B-spline estimate based on equidistant knots with the last few coefficients constrained to be equal works reasonably well. We also propose a constrained maximum penalised partial likelihood approach with the constraints removed through re-parameterisation. A novel feature of the proposed method is that it is based on a log product of spacings penalty rather than the usual roughness penalty, which makes selection of smoothing parameter easier. The utility of the proposed methods is demonstrated using a real data example and simulations.

SummaryA brief survey on methods to handle non‐proportional hazards in survival analysis is given with emphasis on short‐term and long‐term hazard ratio modelling. A drawback of the existing model of this nature is that except at time zero or infinity, the hazard ratio for a unit increase in the value of a covariate depends on the starting value. With two or more covariates, the hazard ratio for a unit increase in one covariate with other covariates held fixed depends in an unintended way on the values of the other covariates. We propose an alternative way to model short‐term and long‐term hazard ratios without the above drawbacks through a judicious choice of covariate‐time interactions. Under the new model, it is easier to describe the time‐varying effect of each covariate on the hazard. Nonparametric maximum likelihood estimation for the new model can be carried out in the same way as for the existing model. We also propose a product version of the existing model, which overcomes its second drawback but not the first. The advocated covariate–time interaction model provides a better fit to the Veterans Administration lung cancer data set than the original and product versions of the existing model.

In survival analysis, one way to deal with non-proportional hazards is to model short-term and long-term hazard ratios. The existing model of this nature has no control over how fast the hazard ratio is changing over time. We add a parameter to the existing model to allow the hazard ratio to change over time at different speed. A nonparametric maximum likelihood approach is used to estimate the model parameters. The existing model is a special case of the extended model when the speed parameter is 0, which leads naturally to a way of testing the adequacy of the existing model. Simulation results show that there can be substantial bias in the estimation of the short-term and long-term hazard ratio if the speed parameter is fixed incorrectly at 0 rather than estimated. The extended model is fitted to three real data sets to shed new insights, including the observation that converging hazards does not necessarily imply the odds are proportional.

In linear quantile regression, the regression coefficients for different quantiles are typically estimated separately. Efforts to improve the efficiency of estimators are often based on assumptions of commonality among the slope coefficients. We propose instead a two-stage procedure whereby the regression coefficients are first estimated separately and then smoothed over quantile level. Due to the strong correlation between coefficient estimates at nearby quantile levels, existing bandwidth selectors will pick bandwidths that are too small. To remedy this, we use 10-fold cross-validation to determine a common bandwidth inflation factor for smoothing the intercept as well as slope estimates. Simulation results suggest that the proposed method is effective in pooling information across quantile levels, resulting in estimates that are typically more efficient than the separately obtained estimates and the interquantile shrinkage estimates derived using a fused penalty function. The usefulness of the proposed method is demonstrated in a real data example.

To adjust the quantile function estimated using a parametric model, the parametric function is composed with the quantile function of the probability integral transformed data. One round of bandwidth selection suffices as adjustments at all quantile levels can be obtained by smoothing the same set of probability integral transformed data. This is in contrast to the customary additive adjustment which requires the user to transform the data differently for estimating different quantiles. Another advantage of the proposed method is that it yields a diagnostic plot useful in assessing the goodness of fit of the assumed model. Compared with the additive approach, function compositional adjustment pays more attention to the fact that the quantity being estimated is a quantile function. As a result, it enjoys intrinsically some desirable properties such as range preservation and invariance to increasing transformation. It is also more amenable to the study of monotonicity as the composition of two quantile functions is a quantile function. In further support of the new adjustment method, Taylor series approximation and results from three simulation studies suggest that the adjusted estimator is robust to model misspecification, and can be more efficient than direct nonparametric estimation. We illustrate the proposed adjustment method using two examples from water resources and human biomonitoring studies.

The article develops a hybrid variational Bayes (VB) algorithm that combines the mean-field and stochastic linear regression fixed-form VB methods. The new estimation algorithm can be used to approximate any posterior without relying on conjugate priors. We propose a divide and recombine strategy for the analysis of large datasets, which partitions a large dataset into smaller subsets and then combines the variational distributions that have been learned in parallel on each separate subset using the hybrid VB algorithm. We also describe an efficient model selection strategy using cross-validation, which is straightforward to implement as a by-product of the parallel run. The proposed method is applied to fitting generalized linear mixed models. The computational efficiency of the parallel and hybrid VB algorithm is demonstrated on several simulated and real datasets. Supplementary material for this article is available online.

SummaryPooling of data is often carried out to protect privacy or to save cost, with the claimed advantage that it does not lead to much loss of efficiency. We argue that this does not give the complete picture as the estimation of different parameters is affected to different degrees by pooling. We establish a ladder of efficiency loss for estimating the mean, variance, skewness and kurtosis, and more generally multivariate joint cumulants, in powers of the pool size. The asymptotic efficiency of the pooled data non‐parametric/parametric maximum likelihood estimator relative to the corresponding unpooled data estimator is reduced by a factor equal to the pool size whenever the order of the cumulant to be estimated is increased by one. The implications of this result are demonstrated in case–control genetic association studies with interactions between genes. Our findings provide a guideline for the discriminate use of data pooling in practice and the assessment of its relative efficiency. As exact maximum likelihood estimates are difficult to obtain if the pool size is large, we address briefly how to obtain computationally efficient estimates from pooled data and suggest Gaussian estimation and non‐parametric maximum likelihood as two feasible methods.

SummaryUsing data collected from the ‘Sequenced treatment alternatives to relieve depression’ study, we use logistic regression to predict whether a patient will respond to treatment on the basis of early symptom change and patient characteristics. Model selection criteria such as the Akaike information criterion AIC and mean-squared-error of prediction MSEP may not be appropriate if the aim is to predict with a high degree of certainty who will respond or not respond to treatment. Towards this aim, we generalize the definition of the positive and negative predictive value curves to the case of multiple predictors. We point out that it is the ordering rather than the precise values of the response probabilities which is important, and we arrive at a unified approach to model selection via two-sample rank tests. To avoid overfitting, we define a cross-validated version of the positive and negative predictive value curves and compare these curves after smoothing for various models. When applied to the study data, we obtain a ranking of models that differs from those based on AIC and MSEP, as well as a tree-based method and regularized logistic regression using a lasso penalty. Our selected model performs consistently well for both 4-week-ahead and 7-week-ahead predictions.

Human biomonitoring of exposure to environmental chemicals is important. Individual monitoring is not viable because of low individual exposure level or insufficient volume of materials and the prohibitive cost of taking measurements from many subjects. Pooling of samples is an efficient and cost‐effective way to collect data. Estimation is, however, complicated as individual values within each pool are not observed but are only known up to their average or weighted average. The distribution of such averages is intractable when the individual measurements are lognormally distributed, which is a common assumption. We propose to replace the intractable distribution of the pool averages by a Gaussian likelihood to obtain parameter estimates. If the pool size is large, this method produces statistically efficient estimates, but regardless of pool size, the method yields consistent estimates as the number of pools increases. An empirical Bayes (EB) Gaussian likelihood approach, as well as its Bayesian analog, is developed to pool information from various demographic groups by using a mixed‐effect formulation. We also discuss methods to estimate the underlying mean–variance relationship and to select a good model for the means, which can be incorporated into the proposed EB or Bayes framework. By borrowing strength across groups, the EB estimator is more efficient than the individual group‐specific estimator. Simulation results show that the EB Gaussian likelihood estimates outperform a previous method proposed for the National Health and Nutrition Examination Surveys with much smaller bias and better coverage in interval estimation, especially after correction of bias. Copyright © 2014 John Wiley & Sons, Ltd.

Using data collected from the ‘Sequenced treatment alternatives to relieve depression’ study, we use logistic regression to predict whether a patient will respond to treatment on the basis of early symptom change and patient characteristics. Model selection criteria such as the Akaike information criterion AIC and mean-squared-error of prediction MSEP may not be appropriate if the aim is to predict with a high degree of certainty who will respond or not respond to treatment. Towards this aim, we generalize the definition of the positive and negative predictive value curves to the case of multiple predictors. We point out that it is the ordering rather than the precise values of the response probabilities which is important, and we arrive at a unified approach to model selection via two-sample rank tests. To avoid overfitting, we define a cross-validated version of the positive and negative predictive value curves and compare these curves after smoothing for various models. When applied to the study data, we obtain a ranking of models that differs from those based on AIC and MSEP, as well as a tree-based method and regularized logistic regression using a lasso penalty. Our selected model performs consistently well for both 4-week-ahead and 7-week-ahead predictions.

There is much recent interest in finding rare genetic variants associated with various diseases. Owing to the scarcity of rare mutations, single-variant analyses often lack power. To enable pooling of information across variants, we use a random effect formulation within a retrospective modeling framework that respects the retrospective data collecting mechanism of case–control studies. More concretely, we model the control allele frequencies of the variants as random effects, and the systematic differences between the case and control frequencies as fixed effects, resulting in a mixed model. The use of Poisson approximation and gamma-distributed random effects results in a generalized negative binomial distribution for the joint distribution of the control and case frequencies. Variants are selected by conducting stepwise likelihood ratio tests. The superiority of the proposed method over two existing variant selection methods is demonstrated in a simulation study. The effects of non-gamma random effects and correlated variants are also found to be not too detrimental in the simulation study. When the proposed procedure is applied to identify rare variants associated with obesity, it identifies one additional variant not picked up by existing methods.

Haplotype information could lead to more powerful tests of genetic association than single-locus analyses but it is not easy to estimate haplotype frequencies from genotype data due to phase ambiguity. The challenge is compounded when individuals are pooled together to save costs or to increase sample size, which is crucial in the study of rare variants. Existing expectationmaximization type algorithms are slow and cannot cope with large pool size or long haplotypes. We show that by collapsing the total allele frequencies of each pool suitably, the maximum likelihood estimates of haplotype frequencies based on the collapsed data can be calculated very quickly regardless of pool size and haplotype length. We provide a running time analysis to demonstrate the considerable savings in time that the collapsed data method can bring. The method is particularly well suited to estimating certain union probabilities useful in the study of rare variants. We provide theoretical and empirical evidence to suggest that the proposed estimation method will not suffer much loss in efficiency if the variants are rare. We use the method to analyze re-sequencing data collected from a case control study involving 148 obese persons and 150 controls. Focusing on a region containing 25 rare variants around theMGLL gene, our method selects three rare variants as potentially causal. This is more parsimonious than the 12 variants selected by a recently proposed covering method. From another set of 32 rare variants aroundthe FAAH gene, we discover an interesting potential interaction between two of them. Copyright (c) 2012 John Wiley & Sons, Ltd.

The article develops a hybrid Variational Bayes algorithm that combines the mean-field and fixed-form Variational Bayes methods. The new estimation algorithm can be used to approximate any posterior without relying on conjugate priors. We propose a divide and recombine strategy for the analysis of large datasets, which partitions a large dataset into smaller pieces and then combines the variational distributions that have been learnt in parallel on each separate piece using the hybrid Variational Bayes algorithm. The proposed method is applied to fitting generalized linear mixed models. The computational efficiency of the parallel and hybrid Variational Bayes algorithm is demonstrated on several simulated and real datasets.

Background Pooling is a cost effective way to collect data for genetic association studies, particularly for rare genetic variants. It is of interest to estimate the haplotype frequencies, which contain more information than single locus statistics. By viewing the pooled genotype data as incomplete data, the expectation-maximization (EM) algorithm is the natural algorithm to use, but it is computationally intensive. A recent proposal to reduce the computational burden is to make use of database information to form a list of frequently occurring haplotypes, and to restrict the haplotypes to come from this list only in implementing the EM algorithm. There is, however, the danger of using an incorrect list, and there may not be enough database information to form a list externally in some applications. Results We investigate the possibility of creating an internal list from the data at hand. One way to form such a list is to collapse the observed total minor allele frequencies to “zero” or “at least one”, which is shown to have the desirable effect of amplifying the haplotype frequencies. To improve coverage, we propose ways to add and remove haplotypes from the list, and a benchmarking method to determine the frequency threshold for removing haplotypes. Simulation results show that the EM estimates based on a suitably augmented and trimmed collapsed data list (ATCDL) perform satisfactorily. In two scenarios involving 25 and 32 loci respectively, the EM-ATCDL estimates outperform the EM estimates based on other lists as well as the collapsed data maximum likelihood estimates. Conclusions The proposed augmented and trimmed CD list is a useful list for the EM algorithm to base upon in estimating the haplotype distributions of rare variants. It can handle more markers and larger pool size than existing methods, and the resulting EM-ATCDL estimates are more efficient than the EM estimates based on other lists.

MOTIVATION:It has been claimed in the literature that pooling DNA samples is efficient in estimating haplotype frequencies. There is, however, no theoretical justification based on calculation of statistical efficiency. In fact, the limited evidence given so far is based on simulation studies with small numbers of loci. With rapid advance in technology, it is of interest to see if pooling is still efficient when the number of loci increases.METHODS:Instead of resorting to simulation studies, we make use of asymptotic statistical theory to perform exact calculation of the efficiency of pooling relative to no pooling in the estimation of haplotype frequencies. As an intermediate step, we use the log-linear formulation of the haplotype probabilities and derive the asymptotic variance-covariance matrix of the maximum likelihood estimators of the canonical parameters of the log-linear model.RESULTS:Based on our calculations under linkage equilibrium, pooling can suffer huge loss in efficiency relative to no pooling when there are more than three independent loci and the alleles are not rare. Pooling works better for rare alleles. In particular, if all the minor allele frequencies are 0.05, pooling maintains an advantage over no pooling until the number of independent loci reaches 6. High linkage disequilibrium effectively reduces the number of independent loci by ruling out certain haplotypes from occurring. Similar calculations of efficiency for the case of no pooling justify the common belief that it is not worthwhile to use molecular methods to resolve the phase ambiguity of individual genotype data.AVAILABILITY:The R codes for the calculation are available at http://www.stat.nus.edu.sg/∼staxj/poolingCONTACT:stakuka@nus.edu.sg.

MOTIVATION:The multifactor-dimensionality reduction (MDR) method has been widely used in multi-locus interaction analysis. It reduces dimensionality by partitioning the multi-locus genotypes into a high-risk group and a low-risk group according to whether the genotype-specific risk ratio exceeds a fixed threshold or not. Alternatively, one can maximize the chi(2) value exhaustively over all possible ways of partitioning the multi-locus genotypes into two groups, and we aim to show that this is computationally feasible.METHODS:We advocate finding the optimal MDR (OMDR) that would have resulted from an exhaustive search over all possible ways of partitioning the multi-locus genotypes into two groups. It is shown that this optimal MDR can be obtained efficiently using an ordered combinatorial partitioning (OCP) method, which differs from the existing MDR method in the use of a data-driven rather than fixed threshold. The generalized extreme value distribution (GEVD) theory is applied to find the optimal order of gene combination and assess statistical significance of interactions.RESULTS:The computational complexity of OCP strategy is linear in the number of multi-locus genotypes in contrast with an exponential order for the naive exhaustive search strategy. Simulation studies show that OMDR can be more powerful than MDR with substantial power gain possible when the partitioning of OMDR is different from that of MDR. The analysis results of a breast cancer dataset show that the use of GEVD accelerates the determination of interaction order and reduces the time cost for P-value calculation by more than 10-fold.AVAILABILITY:C++ program is available at http://home.ustc.edu.cn/~zhanghan/ocp/ocp.html

Objective: There are currently no clinically useful assessments that can reliably predict-early in treatment-whether a particular depressed patient will respond to a particular antidepressant. We explored the possibility of using baseline features and early symptom change to predict which patients will and which patients will not respond to treatment.Method: Participants were 2,280 outpatients enrolled in the Sequenced Treatment Alternatives to Relieve Depression (STAR*D) study who had complete 16-item Quick Inventory of Depressive Symptomatology-self-report (QIDS-SR16) records at baseline, week 2, and week 6 (primary outcome) of treatment with citalopram. Response was defined as a >= 50% reduction in QIDS-SR16 score by week 6. By developing a recursive subsetting algorithm, we used both baseline variables and change in QIDS-SR16 scores from baseline to week 2 to predict response/nonresponse to treatment for as many patients as possible with controlled accuracy, while reserving judgment for the rest.Results: Baseline variables by themselves were not clinically useful predictors, whereas symptom change from baseline to week 2 identified 280 nonresponders, of which 227 were true nonresponders. By subsetting recursively according to both baseline features and symptom change, we were able to identify 505 nonresponders, of which 403 were true nonresponders, to achieve a clinically meaningful negative predictive value of 0.8, which was upheld in cross-validation analyses.Conclusions: Recursive subsetting based on baseline features and early symptom change allows predictions of nonresponse that are sufficiently certain for clinicians to spare identified patients from prolonged exposure to ineffective treatment, thereby personalizing depression management and saving time and cost. Trial Registration: clinicaltrials.gov Identifier: NCT00021528 J Clin Psychiatry 2010;71(11):1502 1508 (C) Copyright 2010 Physicians Postgraduate Press, Inc.

The analysis of clustered binary data is a common task in many areas of application. Parametric approaches to the analysis of such data are numerous, but there has been much recent interest in nonparametric and semiparametric approaches. When cluster sizes are unequal, an assumption is often made of compatibility of marginal distributions in order for semiparametric approaches to be developed when there is little replication for different cluster sizes. Here, we use the marginal compatibility assumption to extend flexible semiparametric Bayesian methods able to shrink towards a “parametric backbone” to the situation where there are few replicated observations for distinct cluster sizes and each distinct value of a covariate. A motivating application is the analysis of developmental toxicology data where pregnant laboratory animals are exposed to a dose of some potentially toxic compound and interest lies in describing the distribution, as a function of the dose level, of the number of fetuses exhibiting some characteristic abnormality. Flexible semiparametric methods are required here, as the data typically exhibit overdispersion and complex structure. We also consider a further extension appropriate to the analysis of clustered binary data in the situation where there is little or no replication for distinct covariate values.

>The survival analysis literature has always lagged behind the categorical data literature in developing methods to analyze clustered or multivariate data. While estimators based on working

SummaryIn a regression model, the joint distribution for each finite sample of units is determined by a function px(y) depending only on the list of covariate values x=(x(u1),…,x(un)) on the sampled units. No random sampling of units is involved. In biological work, random sampling is frequently unavoidable, in which case the joint distribution p(y,x) depends on the sampling scheme. Regression models can be used for the study of dependence provided that the conditional distribution p(y|x) for random samples agrees with px(y) as determined by the regression model for a fixed sample having a non-random configuration x. The paper develops a model that avoids the concept of a fixed population of units, thereby forcing the sampling plan to be incorporated into the sampling distribution. For a quota sample having a predetermined covariate configuration x, the sampling distribution agrees with the standard logistic regression model with correlated components. For most natural sampling plans such as sequential or simple random sampling, the conditional distribution p(y|x) is not the same as the regression distribution unless px(y) has independent components. In this sense, most natural sampling schemes involving binary random-effects models are biased. The implications of this formulation for subject-specific and population-averaged procedures are explored.