A comparison of several methods for analyzing censored data

The purpose of this study was to compare the performance of several methods for statistically analyzing censored datasets [i.e. datasets that contain measurements that are less than the field limit-of-detection (LOD)] when estimating the 95th percentile and the mean of right-skewed occupational exposure data. The methods examined were several variations on the maximum likelihood estimation (MLE) and log-probit regression (LPR) methods, the common substitution methods, several non-parametric (NP) quantile methods for the 95th percentile and the NP Kaplan-Meier (KM) method. Each method was challenged with computer-generated censored datasets for a variety of plausible scenarios where the following factors were allowed to vary randomly within fairly wide ranges: the true geometric standard deviation, the censoring point or LOD and the sample size. This was repeated for both a single-laboratory scenario (i.e. single LOD) and a multiple-laboratory scenario (i.e. three LODs) as well as a single lognormal distribution scenario and a contaminated lognormal distribution scenario. Each method was used to estimate the 95th percentile and mean for the censored datasets (the NP quantile methods estimated only the 95th percentile). For each scenario, the method bias and overall imprecision (as indicated by the root mean square error or rMSE) were calculated for the 95th percentile and mean. No single method was unequivocally superior across all scenarios, although nearly all of the methods excelled in one or more scenarios. Overall, only the MLE- and LPR-based methods performed well across all scenarios, with the robust versions generally showing less bias than the standard versions when challenged with a contaminated lognormal distribution and multiple LODs. All of the MLE- and LPR-based methods were remarkably robust to departures from the lognormal assumption, nearly always having lower rMSE values than the NP methods for the exposure scenarios postulated. In general, the MLE methods tended to have smaller rMSE values than the LPR methods, particularly for the small sample size scenarios. The substitution methods tended to be strongly biased, but in some scenarios had the smaller rMSE values, especially for sample sizes <20. Surprisingly, the various NP methods were not as robust as expected, performing poorly in the contaminated distribution scenarios for both the 95th percentile and the mean. In conclusion, when using the rMSE rather than bias as the preferred comparison metric, the standard MLE method consistently outperformed the so-called robust variations of the MLE-based and LPR-based methods, as well as the various NP methods, for both the 95th percentile and the mean. When estimating the mean, the standard LPR method tended to outperform the robust LPR-based methods. Whenever bias is the main consideration, the robust MLE-based methods should be considered. The KM method, currently hailed by some as the preferred method for estimating the mean when the lognormal distribution assumption is questioned, did not perform well for either the 95th percentile or mean and is not recommended.     

有序多分类重复测量资料的广义线性混合效应模型分析

目的 探讨广义线性混合效应模型在有序多分类重复测量资料分析中的应用及SAS9.1的GLIMMIX和NLMIXED过程实现.方法 为了评价某新药治疗糖尿病神经病变的临床疗效,采用以安慰剂为对照的随机双盲临床试验.在各个随访时间记录各受试者的神经病变主觉症状总分,并根据减分率评定疗效.建立广义线性混合效应模型,并分别用线性化法和数值法积分近似法进行参数估计,利用SAS中的GLIMMIX和NLMIXED过程得以实现.结果 2种参数估计方法 结果 很接近.疗效的组间差别有统计学意义（P〈0.000 1）,试验组疗效优于安慰剂组;各个疗程间的疗效差别有统计学意义（P〈0.000 1）,且疗程越大疗效越好; 治疗前神经病变主觉症状总分对疗效有影响（P=0.061 3,接近显著性水平）,其值越高,越容易治愈,提示病情严重的患者相比病情轻微的患者治愈效果更好.另外用数值法积分近似法还给出了随机截距和随机斜率的统计显著性检验.结论 采用广义线性混合效应模型对有序多分类重复测量临床资料进行统计分析,可以更客观的进行药物疗效评价.     

重复测量资料非线性混合效应模型应用与实现

目的探讨重复测量资料非线性分析技术、SAS软件NLMIXED过程及在群体药动学的应用.方法结合重复测量数据特点,采用最大似然原理进行参数估计,建立非线性混合效应参数模型.结果该模型不仅考虑了个体内和个体间变异,而且也考虑了参数间的非线性,允许固定效应和随机效应进入模型的非线性部分;可方便地分析随机缺失等非均衡数据;有助于引入其他解释变量时最佳模型的选择,更客观地解释其对代谢过程的影响.结论当重复测量资料不满足线性条件时,使用非线性混合效应模型能更客观地反映原数据特征,挖掘资料蕴藏的信息,弥补线性理论分析非线性重复测量资料之不足.     

Penalized likelihood estimation for semiparametric mixed models,with application to alcohol treatment research

In this article, we implement a practical computational method for various semiparametric mixed effects models, estimating nonlinear functions by penalized splines. We approximate the integration of the penalized likelihood with respect to random effects with the use of adaptive Gaussian quadrature, which we can conveniently implement in SAS procedure NLMIXED. We carry out the selection of smoothing parameters through approximated generalized cross‐validation scores. Our method has two advantages: (1) the estimation is more accurate than the current available quasi‐likelihood method for sparse data, for example, binary data; and (2) it can be used in fitting more sophisticated models. We show the performance of our approach in simulation studies with longitudinal outcomes from three settings: binary, normal data after Box–Cox transformation, and count data with log‐Gamma random effects. We also develop an estimation method for a longitudinal two‐part nonparametric random effects model and apply it to analyze repeated measures of semicontinuous daily drinking records in a randomized controlled trial of topiramate. Copyright © 2012 John Wiley & Sons, Ltd.     

A likelihood reformulation method in non-normal random effects models

In this paper, we propose a practical computational method to obtain the maximum likelihood estimates (MLE) for mixed models with non-normal random effects. By simply multiplying and dividing a standard normal density, we reformulate the likelihood conditional on the non-normal random effects to that conditional on the normal random effects. Gaussian quadrature technique, conveniently implemented in SAS Proc NLMIXED, can then be used to carry out the estimation process. Our method substantially reduces computational time, while yielding similar estimates to the probability integral transformation method (J. Comput. Graphical Stat. 2006; 15:39-57). Furthermore, our method can be applied to more general situations, e.g. finite mixture random effects or correlated random effects from Clayton copula. Simulations and applications are presented to illustrate our method.     

Joint mixed‐effects models for causal inference with longitudinal data

Causal inference with observational longitudinal data and time‐varying exposures is complicated due to the potential for time‐dependent confounding and unmeasured confounding. Most causal inference methods that handle time‐dependent confounding rely on either the assumption of no unmeasured confounders or the availability of an unconfounded variable that is associated with the exposure (eg, an instrumental variable). Furthermore, when data are incomplete, validity of many methods often depends on the assumption of missing at random. We propose an approach that combines a parametric joint mixed‐effects model for the study outcome and the exposure with g‐computation to identify and estimate causal effects in the presence of time‐dependent confounding and unmeasured confounding. G‐computation can estimate participant‐specific or population‐average causal effects using parameters of the joint model. The joint model is a type of shared parameter model where the outcome and exposure‐selection models share common random effect(s). We also extend the joint model to handle missing data and truncation by death when missingness is possibly not at random. We evaluate the performance of the proposed method using simulation studies and compare the method to both linear mixed‐ and fixed‐effects models combined with g‐computation as well as to targeted maximum likelihood estimation. We apply the method to an epidemiologic study of vitamin D and depressive symptoms in older adults and include code using SAS PROC NLMIXED software to enhance the accessibility of the method to applied researchers.     

分类重复测量数据的非线性混合效应模型

[目的]探讨分类重复测量数据的非线性混合效应模型及SAS8.0软件NLMIXED过程实现。[方法]直接拟合分类反应变量的非线性概率模型,结合重复测量资料的特点,采用附加高斯积分来获得最大似然的参数估计。[结果]非线性混合效应模型能很好地拟合分类反应变量的重复测量资料,它允许固定效应和随机效应进入模型的非线性部分,可方便地分析随机缺失等非均衡数据。[结论]分类反应变量重复测量资料的非线性混合效应模型分析结果合理、容易解释,为分类重复测量资料提供一种新的分析思路。     

A note on point estimation of the hazard ratio in exponential distributions

The maximum likelihood estimator (MLE) of the ratio of the hazard rates in two exponential distributions is biased. This bias can be important when sample sizes are small or the ratio of these two hazard rates is large. When there is either no censoring or type II censoring, we propose using the uniformly minimum variance unbiased estimator (UMVUE). We show that using the UMVUE instead of the MLE reduces the mean-squared error (MSE). We have found that the UMVUE always has MSEs smaller than the MLE. We have also found that the UMVUE leads to important reductions in the MSE when the sample size used to calculate the hazard rate in the numerator of the hazard ratio is small (say less than or equal to 20) regardless of the sample size in the denominator. In the presence of type I censoring, the proposed estimator and the MLE are both biased. On the basis of a Monte Carlo study, however, we obtain similar reductions in the MSE using the UMVUE, as for no censoring or type II censoring.     

A distribution-based multiple imputation method for handling bivariate pesticide data with values below the limit of detection

Background

Environmental and biomedical researchers frequently encounter laboratory data constrained by a lower limit of detection (LOD). Commonly used methods to address these left-censored data, such as simple substitution of a constant for all values < LOD, may bias parameter estimation. In contrast, multiple imputation (MI) methods yield valid and robust parameter estimates and explicit imputed values for variables that can be analyzed as outcomes or predictors.

Objective

In this article we expand distribution-based MI methods for left-censored data to a bivariate setting, specifically, a longitudinal study with biological measures at two points in time.

Methods

We have presented the likelihood function for a bivariate normal distribution taking into account values < LOD as well as missing data assumed missing at random, and we use the estimated distributional parameters to impute values < LOD and to generate multiple plausible data sets for analysis by standard statistical methods. We conducted a simulation study to evaluate the sampling properties of the estimators, and we illustrate a practical application using data from the Community Participatory Approach to Measuring Farmworker Pesticide Exposure (PACE3) study to estimate associations between urinary acephate (APE) concentrations (indicating pesticide exposure) at two points in time and self-reported symptoms.

Results

Simulation study results demonstrated that imputed and observed values together were consistent with the assumed and estimated underlying distribution. Our analysis of PACE3 data using MI to impute APE values < LOD showed that urinary APE concentration was significantly associated with potential pesticide poisoning symptoms. Results based on simple substitution methods were substantially different from those based on the MI method.

Conclusions

The distribution-based MI method is a valid and feasible approach to analyze bivariate data with values < LOD, especially when explicit values for the nondetections are needed. We recommend the use of this approach in environmental and biomedical research.     

Joint models for efficient estimation in proportional hazards regression models

In survival studies, information lost through censoring can be partially recaptured through repeated measures data which are predictive of survival. In addition, such data may be useful in removing bias in survival estimates, due to censoring which depends upon the repeated measures. Here we investigate joint models for survival T and repeated measurements Y, given a vector of covariates Z. Mixture models indexed as f (T/Z) f (Y/T,Z) are well suited for assessing covariate effects on survival time. Our objective is efficiency gains, using non-parametric models for Y in order to avoid introducing bias by misspecification of the distribution for Y. We model (T/Z) as a piecewise exponential distribution with proportional hazards covariate effect. The component (Y/T,Z) has a multinomial model. The joint likelihood for survival and longitudinal data is maximized, using the EM algorithm. The estimate of covariate effect is compared to the estimate based on the standard proportional hazards model and an alternative joint model based estimate. We demonstrate modest gains in efficiency when using the joint piecewise exponential joint model. In a simulation, the estimated efficiency gain over the standard proportional hazards model is 6.4 per cent. In clinical trial data, the estimated efficiency gain over the standard proportional hazards model is 10.2 per cent.     

Bias due to missing exposure data using complete-case analysis in the proportional hazards regression model

We studied bias due to missing exposure data in the proportional hazards regression model when using complete-case analysis (CCA). Eleven missing data scenarios were considered: one with missing completely at random (MCAR), four missing at random (MAR), and six non-ignorable missingness scenarios, with a variety of hazard ratios, censoring fractions, missingness fractions and sample sizes. When missingness was MCAR or dependent only on the exposure, there was negligible bias (2-3 per cent) that was similar to the difference between the estimate in the full data set with no missing data and the true parameter. In contrast, substantial bias occurred when missingness was dependent on outcome or both outcome and exposure. For models with hazard ratio of 3.5, a sample size of 400, 20 per cent censoring and 40 per cent missing data, the relative bias for the hazard ratio ranged between 7 per cent and 64 per cent. We observed important differences in the direction and magnitude of biases under the various missing data mechanisms. For example, in scenarios where missingness was associated with longer or shorter follow-up, the biases were notably different, although both mechanisms are MAR. The hazard ratio was underestimated (with larger bias) when missingness was associated with longer follow-up and overestimated (with smaller bias) when associated with shorter follow-up. If it is known that missingness is associated with a less frequently observed outcome or with both the outcome and exposure, CCA may result in an invalid inference and other methods for handling missing data should be considered.     

Joint modeling longitudinal semi‐continuous data and survival,with application to longitudinal medical cost data

It has been increasingly common to analyze simultaneously repeated measures and time to failure data. In this paper we propose a joint model when the repeated measures are semi‐continuous, characterized by the presence of a large portion of zero values, as well as right skewness of non zero (positive) values. Examples include monthly medical costs, car insurance annual claims, or annual number of hospitalization days. A random effects two‐part model is used to describe respectively the odds of being positive and the level of positive values. The random effects from the two‐part model are then incorporated in the hazard of the failure time to form the joint model. The estimation can be carried out by Gaussian quadrature techniques conveniently implemented in SAS Proc NLMIXED. Our model is applied to longitudinal (monthly) medical costs of 1455 chronic heart‐failure patients from the clinical data repository at the University of Virginia. Copyright © 2008 John Wiley & Sons, Ltd.     

Empirical Bayes Gaussian likelihood estimation of exposure distributions from pooled samples in human biomonitoring

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.     

Mixed treatment comparison of repeated measurements of a continuous endpoint: an example using topical treatments for primary open‐angle glaucoma and ocular hypertension

Mixed treatment comparison (MTC) meta‐analyses estimate relative treatment effects from networks of evidence while preserving randomisation. We extend the MTC framework to allow for repeated measurements of a continuous endpoint that varies over time. We used, as a case study, a systematic review and meta‐analysis of intraocular pressure (IOP) measurements from randomised controlled trials evaluating topical ocular hypotensives in primary open‐angle glaucoma or ocular hypertension because IOP varies over the day and over the treatment course, and repeated measurements are frequently reported. We adopted models for conducting MTC in W inBUGS (The BUGS Project, Cambridge, UK) to allow for repeated IOP measurements and to impute missing standard deviations of the raw data using the predictive distribution from observations with standard deviations. A flexible model with an unconstrained baseline for IOP variations over time and time‐invariant random treatment effects fitted the data well. We also adopted repeated measures models to allow for class effects; assuming treatment effects to be exchangeable within classes slightly improved model fit but could bias estimated treatment effects if exchangeability assumptions were not valid. We enabled all timepoints to be included in the analysis, allowing for repeated measures to increase precision around treatment effects and avoid bias associated with selecting timepoints for meta‐analysis.The methods we developed for modelling repeated measures and allowing for missing data may be adapted for use in other MTC meta‐analyses. Copyright © 2011 John Wiley & Sons, Ltd.     

Empirical Bayes estimates generated in a hierarchical summary ROC analysis agreed closely with those of a full Bayesian analysis

BACKGROUND AND OBJECTIVE: A range of fixed-effect and random-effects meta-analytic methods are available to obtain summary estimates of measures of diagnostic test accuracy. The hierarchical summary receiver operating characteristic (HSROC) model proposed by Rutter and Gatsonis in 2001 represents a general framework for the meta-analysis of diagnostic test studies that allows different parameters to be defined as a fixed effect or random effects within the same model. The Bayesian method used for fitting the model is complex, however, and the model is not widely used. The objective of this report is to show how the model may be fitted using the SAS procedure NLMIXED and to compare the results to the fully Bayesian analysis using an example. METHODS: The HSROC model, its assumptions, and its interpretation are described. The advantages of this model over the usual summary ROC (SROC) regression model are outlined. A complex example is used to compare the estimated SROC curves, expected operating points, and confidence intervals using the alternative approaches to fitting the model. RESULTS: The empirical Bayes estimates obtained using NLMIXED agree closely with those obtained using the fully Bayesian analysis. CONCLUSION: This alternative and more straightforward method for fitting the HSROC model makes the model more accessible to meta-analysts.     

暴露评估中处理未检出值的多检出限回归方法及其应用

目的探索暴露评估中多检出限回归在存在多个检出限的情况下未检出值数据填补中的应用。方法介绍多检出限回归方法中次序统计量的计算方法,并将方法应用与于处理2004年全国淡水鱼重金属镉残留量数据,对未检出值进行填补,并将填补后的参数估计与简单替换的方法进行比较分析。结果多检出限回归填补数据后估计的2004年全国镉残留平均值(mg/kg)为0.01509,3种不同替换值(0、1/2和1倍检出限)的估计结果分别为0.014759、0.015270和0.015781。结论多检出限回归填补估计方法利用检出值的信息对未检出值进行填补,较简单替换法更为合理。通过SAS宏程序可以方便得以实现,值得在食品暴露评估中推广应用。     

重复测量数据的混合模型及其MIXED过程实现—混合线性模型及其SAS软件实现（二）

目的：重复测量数据存在自相关及随机误差分布于不同层次,不宜使用常规分析方法,本文研讨使用混合线性模型及SAS软件实现的分析方法;方法：利用MIXED对多个处理组的重复测量数据进行混合模型分析。结果：通过固定效应与随机效应及对协方差矩阵的估计,使重复测量数据得以合理的分析。结论：MIXED可以有效地,全面地分析重复测量数据。     

A toolkit for measurement error correction,with a focus on nutritional epidemiology

Exposure measurement error is a problem in many epidemiological studies, including those using biomarkers and measures of dietary intake. Measurement error typically results in biased estimates of exposure‐disease associations, the severity and nature of the bias depending on the form of the error. To correct for the effects of measurement error, information additional to the main study data is required. Ideally, this is a validation sample in which the true exposure is observed. However, in many situations, it is not feasible to observe the true exposure, but there may be available one or more repeated exposure measurements, for example, blood pressure or dietary intake recorded at two time points. The aim of this paper is to provide a toolkit for measurement error correction using repeated measurements. We bring together methods covering classical measurement error and several departures from classical error: systematic, heteroscedastic and differential error. The correction methods considered are regression calibration, which is already widely used in the classical error setting, and moment reconstruction and multiple imputation, which are newer approaches with the ability to handle differential error. We emphasize practical application of the methods in nutritional epidemiology and other fields. We primarily consider continuous exposures in the exposure‐outcome model, but we also outline methods for use when continuous exposures are categorized. The methods are illustrated using the data from a study of the association between fibre intake and colorectal cancer, where fibre intake is measured using a diet diary and repeated measures are available for a subset. © 2014 The Authors. Statistics in Medicine Published by John Wiley & Sons, Ltd.     

双反应变量重复测量资料分析及MIXED过程实现

目的探讨双反应变量重复测量资料的分析原理与方法及SAS软件PROCMIXED过程的应用。方法结合双反应变量重复测量数据的特点,采用SAS软件的MIXED过程对其进行分析,建立线性混合效应模型。结果该模型不仅考虑了每个变量多次重复测量结果之间的相关性,也考虑了两个变量之间的相关性,同时还引入固定效应和随机效应,结合数据特征分析,结果更为可信。结论对双反应变量非独立重复测量资料,可以把数据之间的相关性分解为重复测量间相关性和变量间相关性两部分,采用MIXED过程不仅可对其相关性做出明晰深入的分析,且可保证数据分析结果解释更符合实际。     

Censored correlated cytokine concentrations: multivariate Tobit regression using clustered variance estimation

Interest in cytokines as markers for the function of the immune system is increasing. Methods quantifying cytokine concentrations are often subject to detection limits, which lead to non‐detectable observations and censored distributions. When distributions are skewed, geometric mean ratios (GMRs) can be used to describe the relative concentration between two cytokines, and the GMR ratio (GMRR) can be used to compare two groups. The problem is how to estimate GMRRs from censored distributions.We evaluated methods, including simple deletion and substitution, in simulated and real data. One method applies Tobit directly to the censored difference between the two cytokine log‐concentrations (Diff). However, censoring is correlated to the outcome and is therefore not independent. The correlation increases as the correlation between the two log‐concentrations decreases. We propose a Tobit stacking method that uses clustered variance–covariance estimation allowing homogeneous (Stackc) or inhomogeneous (Stackh) variances. We compare it with direct estimation of the bivariate Tobit likelihood function (Bitobit) and multiple imputation. We assess sensitivity to inhomogeneity and non‐normality. Simulations show that deletion and substitution are empirically biased and that Diff has an empirical bias, which increases as the correlation between the log‐concentrations decreases. Estimates from multiple imputation, Stackh and Bitobit are almost identical. The estimates exhibit small empirical bias for both homogeneous and inhomogeneous normal distributions. For skewed mixture and heavy‐tailed distributions, they perform reasonably well if censoring is less than 30%. We recommend these methods to estimate GMRRs. At least one of the methods is available in Stata, R or SAS. Copyright © 2012 John Wiley & Sons, Ltd.