The PX-EM algorithm for fast stable fitting of Henderson's mixed model

This paper presents procedures for implementing the PX-EM algorithm of Liu, Rubin and Wu to compute REML estimates of variance covariance components in Henderson''s linear mixed models. The class of models considered encompasses several correlated random factors having the same vector length e.g., as in random regression models for longitudinal data analysis and in sire-maternal grandsire models for genetic evaluation. Numerical examples are presented to illustrate the procedures. Much better results in terms of convergence characteristics (number of iterations and time required for convergence) are obtained for PX-EM relative to the basic EM algorithm in the random regression.     

Estimating Haplotype Effects for Survival Data

Summary Genetic association studies often investigate the effect of haplotypes on an outcome of interest. Haplotypes are not observed directly, and this complicates the inclusion of such effects in survival models. We describe a new estimating equations approach for Cox's regression model to assess haplotype effects for survival data. These estimating equations are simple to implement and avoid the use of the EM algorithm, which may be slow in the context of the semiparametric Cox model with incomplete covariate information. These estimating equations also lead to easily computable, direct estimators of standard errors, and thus overcome some of the difficulty in obtaining variance estimators based on the EM algorithm in this setting. We also develop an easily implemented goodness‐of‐fit procedure for Cox's regression model including haplotype effects. Finally, we apply the procedures presented in this article to investigate possible haplotype effects of the PAF‐receptor on cardiovascular events in patients with coronary artery disease, and compare our results to those based on the EM algorithm.     

Joint modeling of survival and longitudinal data: likelihood approach revisited

The maximum likelihood approach to jointly model the survival time and its longitudinal covariates has been successful to model both processes in longitudinal studies. Random effects in the longitudinal process are often used to model the survival times through a proportional hazards model, and this invokes an EM algorithm to search for the maximum likelihood estimates (MLEs). Several intriguing issues are examined here, including the robustness of the MLEs against departure from the normal random effects assumption, and difficulties with the profile likelihood approach to provide reliable estimates for the standard error of the MLEs. We provide insights into the robustness property and suggest to overcome the difficulty of reliable estimates for the standard errors by using bootstrap procedures. Numerical studies and data analysis illustrate our points.     

Standard errors for EM estimates in generalized linear models with random effects

A procedure is derived for computing standard errors of EM estimates in generalized linear models with random effects. Quadrature formulas are used to approximate the integrals in the EM algorithm, where two different approaches are pursued, i.e., Gauss-Hermite quadrature in the case of Gaussian random effects and nonparametric maximum likelihood estimation for an unspecified random effect distribution. An approximation of the expected Fisher information matrix is derived from an expansion of the EM estimating equations. This allows for inferential arguments based on EM estimates, as demonstrated by an example and simulations.     

A multilevel mixed effects varying coefficient model with multilevel predictors and random effects for modeling hospitalization risk in patients on dialysis

For patients on dialysis, hospitalizations remain a major risk factor for mortality and morbidity. We use data from a large national database, United States Renal Data System, to model time-varying effects of hospitalization risk factors as functions of time since initiation of dialysis. To account for the three-level hierarchical structure in the data where hospitalizations are nested in patients and patients are nested in dialysis facilities, we propose a multilevel mixed effects varying coefficient model (MME-VCM) where multilevel (patient- and facility-level) random effects are used to model the dependence structure of the data. The proposed MME-VCM also includes multilevel covariates, where baseline demographics and comorbidities are among the patient-level factors, and staffing composition and facility size are among the facility-level risk factors. To address the challenge of high-dimensional integrals due to the hierarchical structure of the random effects, we propose a novel two-step approximate EM algorithm based on the fully exponential Laplace approximation. Inference for the varying coefficient functions and variance components is achieved via derivation of the standard errors using score contributions. The finite sample performance of the proposed estimation procedure is studied through simulations.     

Employing a Monte Carlo Algorithm in Newton-Type Methods for Restricted Maximum Likelihood Estimation of Genetic Parameters

Estimation of variance components by Monte Carlo (MC) expectation maximization (EM) restricted maximum likelihood (REML) is computationally efficient for large data sets and complex linear mixed effects models. However, efficiency may be lost due to the need for a large number of iterations of the EM algorithm. To decrease the computing time we explored the use of faster converging Newton-type algorithms within MC REML implementations. The implemented algorithms were: MC Newton-Raphson (NR), where the information matrix was generated via sampling; MC average information(AI), where the information was computed as an average of observed and expected information; and MC Broyden''s method, where the zero of the gradient was searched using a quasi-Newton-type algorithm. Performance of these algorithms was evaluated using simulated data. The final estimates were in good agreement with corresponding analytical ones. MC NR REML and MC AI REML enhanced convergence compared to MC EM REML and gave standard errors for the estimates as a by-product. MC NR REML required a larger number of MC samples, while each MC AI REML iteration demanded extra solving of mixed model equations by the number of parameters to be estimated. MC Broyden''s method required the largest number of MC samples with our small data and did not give standard errors for the parameters directly. We studied the performance of three different convergence criteria for the MC AI REML algorithm. Our results indicate the importance of defining a suitable convergence criterion and critical value in order to obtain an efficient Newton-type method utilizing a MC algorithm. Overall, use of a MC algorithm with Newton-type methods proved feasible and the results encourage testing of these methods with different kinds of large-scale problem settings.     

Stochastic algorithms for Markov models estimation with intermittent missing data

Multistate Markov models are frequently used to characterize disease processes, but their estimation from longitudinal data is often hampered by complex patterns of incompleteness. Two algorithms for estimating Markov chain models in the case of intermittent missing data in longitudinal studies, a stochastic EM algorithm and the Gibbs sampler, are described. The first can be viewed as a random perturbation of the EM algorithm and is appropriate when the M step is straightforward but the E step is computationally burdensome. It leads to a good approximation of the maximum likelihood estimates. The Gibbs sampler is used for a full Bayesian inference. The performances of the two algorithms are illustrated on two simulated data sets. A motivating example concerned with the modelling of the evolution of parasitemia by Plasmodium falciparum (malaria) in a cohort of 105 young children in Cameroon is described and briefly analyzed.     

Multilevel Latent Class Models with Dirichlet Mixing Distribution

Summary Latent class analysis (LCA) and latent class regression (LCR) are widely used for modeling multivariate categorical outcomes in social science and biomedical studies. Standard analyses assume data of different respondents to be mutually independent, excluding application of the methods to familial and other designs in which participants are clustered. In this article, we consider multilevel latent class models, in which subpopulation mixing probabilities are treated as random effects that vary among clusters according to a common Dirichlet distribution. We apply the expectation‐maximization (EM) algorithm for model fitting by maximum likelihood (ML). This approach works well, but is computationally intensive when either the number of classes or the cluster size is large. We propose a maximum pairwise likelihood (MPL) approach via a modified EM algorithm for this case. We also show that a simple latent class analysis, combined with robust standard errors, provides another consistent, robust, but less‐efficient inferential procedure. Simulation studies suggest that the three methods work well in finite samples, and that the MPL estimates often enjoy comparable precision as the ML estimates. We apply our methods to the analysis of comorbid symptoms in the obsessive compulsive disorder study. Our models' random effects structure has more straightforward interpretation than those of competing methods, thus should usefully augment tools available for LCA of multilevel data.     

Application of dynamic metabolic flux analysis for process modeling: Robust flux estimation with regularization,confidence bounds,and selection of elementary modes

In macroscopic dynamic models of fermentation processes, elementary modes (EM) derived from metabolic networks are often used to describe the reaction stoichiometry in a simplified manner and to build predictive models by parameterizing kinetic rate equations for the EM. In this procedure, the selection of a set of EM is a key step which is followed by an estimation of their reaction rates and of the associated confidence bounds. In this paper, we present a method for the computation of reaction rates of cellular reactions and EM as well as an algorithm for the selection of EM for process modeling. The method is based on the dynamic metabolic flux analysis (DMFA) proposed by Leighty and Antoniewicz (2011, Metab Eng, 13(6), 745–755) with additional constraints, regularization and analysis of uncertainty. Instead of using estimated uptake or secretion rates, concentration measurements are used directly to avoid an amplification of measurement errors by numerical differentiation. It is shown that the regularized DMFA for EM method is significantly more robust against measurement noise than methods using estimated rates. The confidence intervals for the estimated reaction rates are obtained by bootstrapping. For the selection of a set of EM for a given st oichiometric model, the DMFA for EM method is combined with a multiobjective genetic algorithm. The method is applied to real data from a CHO fed-batch process. From measurements of six fed-batch experiments, 10 EM were identified as the smallest subset of EM based upon which the data can be described sufficiently accurately by a dynamic model. The estimated EM reaction rates and their confidence intervals at different process conditions provide useful information for the kinetic modeling and subsequent process optimization.     

Joint models for a primary endpoint and multiple longitudinal covariate processes

Summary .   Joint models are formulated to investigate the association between a primary endpoint and features of multiple longitudinal processes. In particular, the subject-specific random effects in a multivariate linear random-effects model for multiple longitudinal processes are predictors in a generalized linear model for primary endpoints. Li, Zhang, and Davidian (2004, Biometrics 60 , 1–7) proposed an estimation procedure that makes no distributional assumption on the random effects but assumes independent within-subject measurement errors in the longitudinal covariate process. Based on an asymptotic bias analysis, we found that their estimators can be biased when random effects do not fully explain the within-subject correlations among longitudinal covariate measurements. Specifically, the existing procedure is fairly sensitive to the independent measurement error assumption. To overcome this limitation, we propose new estimation procedures that require neither a distributional or covariance structural assumption on covariate random effects nor an independence assumption on within-subject measurement errors. These new procedures are more flexible, readily cover scenarios that have multivariate longitudinal covariate processes, and can be implemented using available software. Through simulations and an analysis of data from a hypertension study, we evaluate and illustrate the numerical performances of the new estimators.     

Extension of the SAEM algorithm for nonlinear mixed models with 2 levels of random effects

This article focuses on parameter estimation of multilevel nonlinearmixed-effects models (MNLMEMs). These models are used to analyzedata presenting multiple hierarchical levels of grouping (clusterdata, clinical trials with several observation periods, ...).The variability of the individual parameters of the regressionfunction is thus decomposed as a between-subject variabilityand higher levels of variability (e.g. within-subject variability).We propose maximum likelihood estimates of parameters of thoseMNLMEMs with 2 levels of random effects, using an extensionof the stochastic approximation version of expectation–maximization(SAEM)–Monte Carlo Markov chain algorithm. The extendedSAEM algorithm is split into an explicit direct expectation–maximization(EM) algorithm and a stochastic EM part. Compared to the originalalgorithm, additional sufficient statistics have to be approximatedby relying on the conditional distribution of the second levelof random effects. This estimation method is evaluated on pharmacokineticcrossover simulated trials, mimicking theophylline concentrationdata. Results obtained on those data sets with either the SAEMalgorithm or the first-order conditional estimates (FOCE) algorithm(implemented in the nlme function of R software) are compared:biases and root mean square errors of almost all the SAEM estimatesare smaller than the FOCE ones. Finally, we apply the extendedSAEM algorithm to analyze the pharmacokinetic interaction oftenofovir on atazanavir, a novel protease inhibitor, from theAgence Nationale de Recherche sur le Sida 107-Puzzle 2 study.A significant decrease of the area under the curve of atazanaviris found in patients receiving both treatments.     

A general maximum likelihood analysis of variance components in generalized linear models

This paper describes an EM algorithm for nonparametric maximum likelihood (ML) estimation in generalized linear models with variance component structure. The algorithm provides an alternative analysis to approximate MQL and PQL analyses (McGilchrist and Aisbett, 1991, Biometrical Journal 33, 131-141; Breslow and Clayton, 1993; Journal of the American Statistical Association 88, 9-25; McGilchrist, 1994, Journal of the Royal Statistical Society, Series B 56, 61-69; Goldstein, 1995, Multilevel Statistical Models) and to GEE analyses (Liang and Zeger, 1986, Biometrika 73, 13-22). The algorithm, first given by Hinde and Wood (1987, in Longitudinal Data Analysis, 110-126), is a generalization of that for random effect models for overdispersion in generalized linear models, described in Aitkin (1996, Statistics and Computing 6, 251-262). The algorithm is initially derived as a form of Gaussian quadrature assuming a normal mixing distribution, but with only slight variation it can be used for a completely unknown mixing distribution, giving a straightforward method for the fully nonparametric ML estimation of this distribution. This is of value because the ML estimates of the GLM parameters can be sensitive to the specification of a parametric form for the mixing distribution. The nonparametric analysis can be extended straightforwardly to general random parameter models, with full NPML estimation of the joint distribution of the random parameters. This can produce substantial computational saving compared with full numerical integration over a specified parametric distribution for the random parameters. A simple method is described for obtaining correct standard errors for parameter estimates when using the EM algorithm. Several examples are discussed involving simple variance component and longitudinal models, and small-area estimation.     

Haplotype Inference for Population Data with Genotyping Errors

Inference of haplotypes is important in genetic epidemiology studies. However, all large genotype data sets have errors due to the use of inexpensive genotyping machines that are fallible and shortcomings in genotyping scoring softwares, which can have an enormous impact on haplotype inference. In this article, we propose two novel strategies to reduce the impact induced by genotyping errors in haplotype inference. The first method makes use of double sampling. For each individual, the “GenoSpectrum” that consists of all possible genotypes and their corresponding likelihoods are computed. The second method is a genotype clustering algorithm based on multi‐genotyping data, which also assigns a “GenoSpectrum” for each individual. We then describe two hybrid EM algorithms (called DS‐EM and MG‐EM) that perform haplotype inference based on “GenoSpectrum” of each individual obtained by double sampling and multi‐genotyping data. Both simulated data sets and a quasi real‐data set demonstrate that our proposed methods perform well in different situations and outperform the conventional EM algorithm and the HMM algorithm proposed by Sun, Greenwood, and Neal (2007, Genetic Epidemiology 31 , 937–948) when the genotype data sets have errors.     

Estimating haplotype frequencies and standard errors for multiple single nucleotide polymorphisms

Estimating haplotype frequencies becomes increasingly important in the mapping of complex disease genes, as millions of single nucleotide polymorphisms (SNPs) are being identified and genotyped. When genotypes at multiple SNP loci are gathered from unrelated individuals, haplotype frequencies can be accurately estimated using expectation-maximization (EM) algorithms (Excoffier and Slatkin, 1995; Hawley and Kidd, 1995; Long et al., 1995), with standard errors estimated using bootstraps. However, because the number of possible haplotypes increases exponentially with the number of SNPs, handling data with a large number of SNPs poses a computational challenge for the EM methods and for other haplotype inference methods. To solve this problem, Niu and colleagues, in their Bayesian haplotype inference paper (Niu et al., 2002), introduced a computational algorithm called progressive ligation (PL). But their Bayesian method has a limitation on the number of subjects (no more than 100 subjects in the current implementation of the method). In this paper, we propose a new method in which we use the same likelihood formulation as in Excoffier and Slatkin's EM algorithm and apply the estimating equation idea and the PL computational algorithm with some modifications. Our proposed method can handle data sets with large number of SNPs as well as large numbers of subjects. Simultaneously, our method estimates standard errors efficiently, using the sandwich-estimate from the estimating equation, rather than the bootstrap method. Additionally, our method admits missing data and produces valid estimates of parameters and their standard errors under the assumption that the missing genotypes are missing at random in the sense defined by Rubin (1976).     

Maximum likelihood estimation in random effects cure rate models with nonignorable missing covariates

We introduce a method of parameter estimation for a random effects cure rate model. We also propose a methodology that allows us to account for nonignorable missing covariates in this class of models. The proposed method corrects for possible bias introduced by complete case analysis when missing data are not missing completely at random and is motivated by data from a pair of melanoma studies conducted by the Eastern Cooperative Oncology Group in which clustering by cohort or time of study entry was suspected. In addition, these models allow estimation of cure rates, which is desirable when we do not wish to assume that all subjects remain at risk of death or relapse from disease after sufficient follow-up. We develop an EM algorithm for the model and provide an efficient Gibbs sampling scheme for carrying out the E-step of the algorithm.     

非线性再生散度随机效应模型的极大似然估计及EM算法

非线性再生散度随机效应模型是指数族非线性随机效应模型和非线性再生散度模型的推广和发展．通过视模型中的随机效应为假想的缺失数据和应用Metropolis-Hastings（MH）算法,提出了模型参数极大似然估计的Monte-Carlo EM（MCEM）算法,并用模拟研究和实例分析说明了该算法的可行性．     

Small-sample inference for incomplete longitudinal data with truncation and censoring in tumor xenograft models

In cancer drug development, demonstrating activity in xenograft models, where mice are grafted with human cancer cells, is an important step in bringing a promising compound to humans. A key outcome variable is the tumor volume measured in a given period of time for groups of mice given different doses of a single or combination anticancer regimen. However, a mouse may die before the end of a study or may be sacrificed when its tumor volume quadruples, and its tumor may be suppressed for some time and then grow back. Thus, incomplete repeated measurements arise. The incompleteness or missingness is also caused by drastic tumor shrinkage (<0.01 cm3) or random truncation. Because of the small sample sizes in these models, asymptotic inferences are usually not appropriate. We propose two parametric test procedures based on the EM algorithm and the Bayesian method to compare treatment effects among different groups while accounting for informative censoring. A real xenograft study on a new antitumor agent, temozolomide, combined with irinotecan is analyzed using the proposed methods.     

Functional mixed effects models

In this article, a new class of functional models in which smoothing splines are used to model fixed effects as well as random effects is introduced. The linear mixed effects models are extended to nonparametric mixed effects models by introducing functional random effects, which are modeled as realizations of zero-mean stochastic processes. The fixed functional effects and the random functional effects are modeled in the same functional space, which guarantee the population-average and subject-specific curves have the same smoothness property. These models inherit the flexibility of the linear mixed effects models in handling complex designs and correlation structures, can include continuous covariates as well as dummy factors in both the fixed or random design matrices, and include the nested curves models as special cases. Two estimation procedures are proposed. The first estimation procedure exploits the connection between linear mixed effects models and smoothing splines and can be fitted using existing software. The second procedure is a sequential estimation procedure using Kalman filtering. This algorithm avoids inversion of large dimensional matrices and therefore can be applied to large data sets. A generalized maximum likelihood (GML) ratio test is proposed for inference and model selection. An application to comparison of cortisol profiles is used as an illustration.     

The mid-domain effect in ectomycorrhizal fungi: range overlap along an elevation gradient on Mount Fuji,Japan

Mid-domain effect (MDE) models predict that the random placement of species'' ranges within a bounded geographical area leads to increased range overlap and species richness in the center of the bounded area. These models are frequently applied to study species-richness patterns of macroorganisms, but the MDE in relation to microorganisms is poorly understood. In this study, we examined the characteristics of the MDE in richness patterns of ectomycorrhizal (EM) fungi, an ecologically important group of soil symbionts. We conducted intensive soil sampling to investigate overlap among species ranges and the applicability of the MDE to EM fungi in four temperate forest stands along an elevation gradient on Mount Fuji, Japan. Molecular analyses using direct sequencing revealed 302 EM fungal species. Of 73 EM fungal species found in multiple stands, 72 inhabited a continuous range along the elevation gradient. The maximum overlap in species range and the highest species richness occurred at elevations in the middle of the gradient. The observed richness pattern also fit within the 95% confidence interval of the mid-domain null model, supporting the role of the MDE in EM fungal richness. Deviation in observed richness from the mean of the mid-domain null estimation was negatively correlated with some environmental factors, including precipitation and soil C/N, indicating that unexplained richness patterns could be driven by these environmental factors. Our results clearly support the existence of microbial species'' ranges along environmental gradients and the potential applicability of the MDE to better understand microbial diversity patterns.     

Online nonlinear sequential Bayesian estimation of a biological wastewater treatment process

Online estimation of unknown state variables is a key component in the accurate modelling of biological wastewater treatment processes due to a lack of reliable online measurement systems. The extended Kalman filter (EKF) algorithm has been widely applied for wastewater treatment processes. However, the series approximations in the EKF algorithm are not valid, because biological wastewater treatment processes are highly nonlinear with a time-varying characteristic. This work proposes an alternative online estimation approach using the sequential Monte Carlo (SMC) methods for recursive online state estimation of a biological sequencing batch reactor for wastewater treatment. SMC is an algorithm that makes it possible to recursively construct the posterior probability density of the state variables, with respect to all available measurements, through a random exploration of the states by entities called ‘particle’. In this work, the simplified and modified Activated Sludge Model No. 3 with nonlinear biological kinetic models is used as a process model and formulated in a dynamic state-space model applied to the SMC method. The performance of the SMC method for online state estimation applied to a biological sequencing batch reactor with online and offline measured data is encouraging. The results indicate that the SMC method could emerge as a powerful tool for solving online state and parameter estimation problems without any model linearization or restrictive assumptions pertaining to the type of nonlinear models for biological wastewater treatment processes.