Small-sample confidence regions in exponential families

This article presents an algorithm for small-sample conditional confidence regions for two or more parameters for any discrete regression model in the generalized linear interactive model family. Regions are constructed by careful inversion of conditional hypothesis tests. This method presupposes the use of approximate or exact techniques for enumerating the sample space for some components of the vector of sufficient statistics conditional on other components. Such enumeration may be performed exactly or by exact or approximate Monte Carlo, including the algorithms of Kolassa and Tanner (1994, Journal of the American Statistical Association 89, 697-702; 1999, Biometrics 55, 246-251). This method also assumes that one can compute certain conditional probabilities for a fixed value of the parameter vector. Because of a property of exponential families, one can use this set of conditional probabilities to directly compute the conditional probabilities associated with any other value of the vector of the parameters of interest. This observation dramatically reduces the computational effort required to invert the hypothesis test to obtain the confidence region. To construct a region with confidence level 1 - alpha, the algorithm begins with a grid of values for the parameters of interest. For each parameter vector on the grid (corresponding to the current null hypothesis), one transforms the initial set of conditional probabilities using exponential tilting and then calculates the p value for this current null hypothesis. The confidence region is the set of parameter values for which the p value is at least alpha.     

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.     

Monte Carlo EM for missing covariates in parametric regression models

We propose a method for estimating parameters for general parametric regression models with an arbitrary number of missing covariates. We allow any pattern of missing data and assume that the missing data mechanism is ignorable throughout. When the missing covariates are categorical, a useful technique for obtaining parameter estimates is the EM algorithm by the method of weights proposed in Ibrahim (1990, Journal of the American Statistical Association 85, 765-769). We extend this method to continuous or mixed categorical and continuous covariates, and for arbitrary parametric regression models, by adapting a Monte Carlo version of the EM algorithm as discussed by Wei and Tanner (1990, Journal of the American Statistical Association 85, 699-704). In addition, we discuss the Gibbs sampler for sampling from the conditional distribution of the missing covariates given the observed data and show that the appropriate complete conditionals are log-concave. The log-concavity property of the conditional distributions will facilitate a straightforward implementation of the Gibbs sampler via the adaptive rejection algorithm of Gilks and Wild (1992, Applied Statistics 41, 337-348). We assume the model for the response given the covariates is an arbitrary parametric regression model, such as a generalized linear model, a parametric survival model, or a nonlinear model. We model the marginal distribution of the covariates as a product of one-dimensional conditional distributions. This allows us a great deal of flexibility in modeling the distribution of the covariates and reduces the number of nuisance parameters that are introduced in the E-step. We present examples involving both simulated and real data.     

Robustified maximum likelihood estimation in generalized partial linear mixed model for longitudinal data

Summary .  In this article, we study the robust estimation of both mean and variance components in generalized partial linear mixed models based on the construction of robustified likelihood function. Under some regularity conditions, the asymptotic properties of the proposed robust estimators are shown. Some simulations are carried out to investigate the performance of the proposed robust estimators. Just as expected, the proposed robust estimators perform better than those resulting from robust estimating equations involving conditional expectation like Sinha (2004, Journal of the American Statistical Association 99, 451–460) and Qin and Zhu (2007, Journal of Multivariate Analysis 98, 1658–1683). In the end, the proposed robust method is illustrated by the analysis of a real data set.     

Multivariate continuation ratio models: connections and caveats

We develop semiparametric estimation methods for a pair of regressions that characterize the first and second moments of clustered discrete survival times. In the first regression, we represent discrete survival times through univariate continuation indicators whose expectations are modeled using a generalized linear model. In the second regression, we model the marginal pairwise association of survival times using the Clayton-Oakes cross-product ratio (Clayton, 1978, Biometrika 65, 141-151; Oakes, 1989, Journal of the American Statistical Association 84, 487-493). These models have recently been proposed by Shih (1998, Biometrics 54, 1115-1128). We relate the discrete survival models to multivariate multinomial models presented in Heagerty and Zeger (1996, Journal of the American Statistical Society 91, 1024-1036) and derive a paired estimating equations procedure that is computationally feasible for moderate and large clusters. We extend the work of Guo and Lin (1994, Biometrics 50, 632-639) and Shih (1998) to allow covariance weighted estimating equations and investigate the impact of weighting in terms of asymptotic relative efficiency. We demonstrate that the multinomial structure must be acknowledged when adopting weighted estimating equations and show that a naive use of GEE methods can lead to inconsistent parameter estimates. Finally, we illustrate the proposed methodology by analyzing psychological testing data previously summarized by TenHave and Uttal (1994, Applied Statistics 43, 371-384) and Guo and Lin (1994).     

Confidence bands for logistic regression with restricted predictor variables

Confidence bands are constructed for the logistic response function when there is an interval restriction on each of the predictor variables. The construction involves application of a general fitting procedure using Scheffé's S-method, described by Casella and Strawderman (1980, Journal of the American Statistical Association 75, 862-868). Specific details are given for the case of one predictor variable, along with details for a fixed-width alternative to the S-method bands. In the one-predictor case, Monte Carlo results suggest that both bands are conservative for small sample sizes, such as N = 25. By N = 200 the S-method's coverage probabilities are seen to attain their nominal levels while the fixed-width bands remain conservative. The procedures are exemplified with data from a genetic toxicology experiment.     

Response dynamics and receptive-field organization of catfish ganglion cells [published erratum appears in J Gen Physiol 1995 Aug;106(2):following 388]

Responses from catfish retinal ganglion cells were evoked by a spot or an annulus of light and were analyzed by a procedure identical to the one used previously to study catfish amacrine cells (Sakai H. M., and K.-I. Naka, 1992. Journal of Neurophysiology. 67:430-442.). In two- input white-noise experiments, a response evoked by simultaneous stimulation of the center and surround was decomposed into the components generated by the center and surround through a process of cross-correlation. The center and surround responses were also decomposed into their linear and nonlinear components so that the response dynamics of the linear and nonlinear components could be measured. We found that the concentric organization of the receptive field was determined by linear components, i.e., the first-order kernels generated by the center and surround were of opposite polarity. Both the center and surround generated second-order kernels with similar signatures, i.e., the second-order components formed a monotonic receptive field. The peak response time of the first- and second-order kernels from the surround was longer by approximately 20 ms than that of the center. Except for the DC potential present in the intracellular responses, almost identical first- and second-order kernels for the center and surround were obtained from both the intracellular response and spike discharges. Thus, information on concentric organization of a receptive field is translated into spike discharges with little loss of information. A train of spike discharges carries, simultaneously, at least four kinds of information: two linear and two nonlinear components, which originate in the receptive field center and the surround. A spike train is not a simple signaling device but is a carrier of complex and multiple signals. Victor, J. D., and R. M. Shapley (1979. Journal of General Physiology. 74:671-687.) discovered similarly that, in the cat retina, static second-order nonlinearity is encoded into spike trains. Results obtained in this study support the thesis that signals generated by the preganglionic cells are translated into spike discharges without major modification and that those signals can be recovered from the spike trains (Sakuranaga, M., Y. Ando, and K.-I. Naka. 1987. Journal of General Physiology. 90:229-259.; Korenberg, M. J., H. M. Sakai, and K.-I. Naka. 1989. Journal of Neurophysiology. 61:1110-1120.). Current injection studies have shown that such signal transmission is possible (Sakai, H. M., and K.-I. Naka, 1988a. Journal of Neurophysiology. 60:1549-1567.; 1990. Journal of Neurophysiology. 63:105-119.).     

Sequential ordinal modeling with applications to survival data

This paper considers the class of sequential ordinal models in relation to other models for ordinal response data. Markov chain Monte Carlo (MCMC) algorithms, based on the approach of Albert and Chib (1993, Journal of the American Statistical Association 88, 669-679), are developed for the fitting of these models. The ideas and methods are illustrated in detail with a real data example on the length of hospital stay for patients undergoing heart surgery. A notable aspect of this analysis is the comparison, based on marginal likelihoods and training sample priors, of several nonnested models, such as the sequential model, the cumulative ordinal model, and Weibull and log-logistic models.     

Adaptive regression splines in the Cox model

We develop a method for constructing adaptive regression spline models for the exploration of survival data. The method combines Cox's (1972, Journal of the Royal Statistical Society, Series B 34, 187-200) regression model with a weighted least-squares version of the multivariate adaptive regressi on spline (MARS) technique of Friedman (1991, Annals of Statistics 19, 1-141) to adaptively select the knots and covariates. The new technique can automatically fit models with terms that represent nonlinear effects and interactions among covariates. Applications based on simulated data and data from a clinical trial for myeloma are presented. Results from the myeloma application identified several important prognostic variables, including a possible nonmonotone relationship with survival in one laboratory variable. Results are compared to those from the adaptive hazard regression (HARE) method of Kooperberg, Stone, and Truong (1995, Journal of the American Statistical Association 90, 78-94).     

A Bayesian approach to nonlinear random effects models

Nonlinear random effects models are considered from the Bayesian point of view. The method of analysis follows closely that of Lindley and Smith (1972, Journal of the Royal Statistical Society, Series B 34, 1-42). The numerical method is related to the EM algorithm.     

A local influence approach applied to binary data from a psychiatric study

Recently, a lot of concern has been raised about assumptions needed in order to fit statistical models to incomplete multivariate and longitudinal data. In response, research efforts are being devoted to the development of tools that assess the sensitivity of such models to often strong but always, at least in part, unverifiable assumptions. Many efforts have been devoted to longitudinal data, primarily in the selection model context, although some researchers have expressed interest in the pattern-mixture setting as well. A promising tool, proposed by Verbeke et al. (2001, Biometrics 57, 43-50), is based on local influence (Cook, 1986, Journal of the Royal Statistical Society, Series B 48, 133-169). These authors considered the Diggle and Kenward (1994, Applied Statistics 43, 49-93) model, which is based on a selection model, integrating a linear mixed model for continuous outcomes with logistic regression for dropout. In this article, we show that a similar idea can be developed for multivariate and longitudinal binary data, subject to nonmonotone missingness. We focus on the model proposed by Baker, Rosenberger, and DerSimonian (1992, Statistics in Medicine 11, 643-657). The original model is first extended to allow for (possibly continuous) covariates, whereafter a local influence strategy is developed to support the model-building process. The model is able to deal with nonmonotone missingness but has some limitations as well, stemming from the conditional nature of the model parameters. Some analytical insight is provided into the behavior of the local influence graphs.     

Causal mediation analyses with rank preserving models

We present a linear rank preserving model (RPM) approach for analyzing mediation of a randomized baseline intervention's effect on a univariate follow-up outcome. Unlike standard mediation analyses, our approach does not assume that the mediating factor is also randomly assigned to individuals in addition to the randomized baseline intervention (i.e., sequential ignorability), but does make several structural interaction assumptions that currently are untestable. The G-estimation procedure for the proposed RPM represents an extension of the work on direct effects of randomized intervention effects for survival outcomes by Robins and Greenland (1994, Journal of the American Statistical Association 89, 737-749) and on intervention non-adherence by Ten Have et al. (2004, Journal of the American Statistical Association 99, 8-16). Simulations show good estimation and confidence interval performance by the proposed RPM approach under unmeasured confounding relative to the standard mediation approach, but poor performance under departures from the structural interaction assumptions. The trade-off between these assumptions is evaluated in the context of two suicide/depression intervention studies.     

Normal approximation diagnostics for the Cox model

We discuss two diagnostic methods for assessing the accuracy of the normal approximated confidence region to the likelihood-based confidence region for the Cox proportional hazards model with censored data. The proposed diagnostic methods are extensions of the contour measures of Hodges (1987, Journal of the American Statistical Association 82, 149-154) and Cook and Tsai (1990, Journal of the American Statistical Association 85, 770-777) and the curvature measures of Jennings (1986, Journal of the American Statistical Association 81, 471-476) and Cook and Tsai (1990). These methods are also illustrated in a study of hepatocyte growth factor in patients with lung cancer and a Mayo Clinic randomized study of participants with primary biliary cirrhosis.     

Multiple comparisons for analyzing dichotomous response

Dichotomous response models are common in many experimental settings. Statistical parameters of interest are typically the probabilities, pi, that an experimental unit will respond at the various treatment levels. Herein, simultaneous procedures are considered for multiple comparisons among these probabilities, with attention directed at construction of simultaneous confidence intervals for various functions of the pi. The inferences are based on the asymptotic normality of the maximum likelihood estimator of pi. Specific applications include all pairwise comparisons and comparisons with a fixed (control) treatment. Monte Carlo evaluations are undertaken to examine the small-sample properties of the various procedures. It is seen that use of the usual estimates of variance consistently leads to less-than-nominal empirical coverage for most sample sizes examined. For very large samples (total size greater than about 300), nominal coverage is achieved. A reformulation of the pairwise comparisons using a construction noted by Beal (1987, Biometrics 43, 941-950) is shown to exhibit generally nominal empirical coverage characteristics, and is recommended for use with small-to-moderate sample sizes.     

Maximum likelihood analysis of logistic regression models with incomplete covariate data and auxiliary information

This article presents a new method for maximum likelihood estimation of logistic regression models with incomplete covariate data where auxiliary information is available. This auxiliary information is extraneous to the regression model of interest but predictive of the covariate with missing data. Ibrahim (1990, Journal of the American Statistical Association 85, 765-769) provides a general method for estimating generalized linear regression models with missing covariates using the EM algorithm that is easily implemented when there is no auxiliary data. Vach (1997, Statistics in Medicine 16, 57-72) describes how the method can be extended when the outcome and auxiliary data are conditionally independent given the covariates in the model. The method allows the incorporation of auxiliary data without making the conditional independence assumption. We suggest tests of conditional independence and compare the performance of several estimators in an example concerning mental health service utilization in children. Using an artificial dataset, we compare the performance of several estimators when auxiliary data are available.     

Generalized hierarchical multivariate CAR models for areal data

In the fields of medicine and public health, a common application of areal data models is the study of geographical patterns of disease. When we have several measurements recorded at each spatial location (for example, information on p>/= 2 diseases from the same population groups or regions), we need to consider multivariate areal data models in order to handle the dependence among the multivariate components as well as the spatial dependence between sites. In this article, we propose a flexible new class of generalized multivariate conditionally autoregressive (GMCAR) models for areal data, and show how it enriches the MCAR class. Our approach differs from earlier ones in that it directly specifies the joint distribution for a multivariate Markov random field (MRF) through the specification of simpler conditional and marginal models. This in turn leads to a significant reduction in the computational burden in hierarchical spatial random effect modeling, where posterior summaries are computed using Markov chain Monte Carlo (MCMC). We compare our approach with existing MCAR models in the literature via simulation, using average mean square error (AMSE) and a convenient hierarchical model selection criterion, the deviance information criterion (DIC; Spiegelhalter et al., 2002, Journal of the Royal Statistical Society, Series B64, 583-639). Finally, we offer a real-data application of our proposed GMCAR approach that models lung and esophagus cancer death rates during 1991-1998 in Minnesota counties.     

Generalized linear mixed models with varying coefficients for longitudinal data

The routinely assumed parametric functional form in the linear predictor of a generalized linear mixed model for longitudinal data may be too restrictive to represent true underlying covariate effects. We relax this assumption by representing these covariate effects by smooth but otherwise arbitrary functions of time, with random effects used to model the correlation induced by among-subject and within-subject variation. Due to the usually intractable integration involved in evaluating the quasi-likelihood function, the double penalized quasi-likelihood (DPQL) approach of Lin and Zhang (1999, Journal of the Royal Statistical Society, Series B61, 381-400) is used to estimate the varying coefficients and the variance components simultaneously by representing a nonparametric function by a linear combination of fixed effects and random effects. A scaled chi-squared test based on the mixed model representation of the proposed model is developed to test whether an underlying varying coefficient is a polynomial of certain degree. We evaluate the performance of the procedures through simulation studies and illustrate their application with Indonesian children infectious disease data.     

Interval estimation of the LD50 based on an up-and-down experiment

It is well known that an up-and-down method can be more efficient than fixed-sample methods in estimating the LD50 of a quantal response curve. A problem that has not been addressed by many is that of obtaining a confidence interval for the LD50 from the up-and-down method. Dixon and Mood (1948, Journal of the American Statistical Association 43, 109-126) proposed a confidence interval using a maximum likelihood approach, but not much is known about its properties. In this paper, a new confidence interval for the LD50 based on turning points is obtained, which uses the concept of phi-mixing. Simulation results indicate that the coverage probabilities of both methods tend to be less than the nominal level unless the sample size is large. Even so, when the tolerance distribution is normal, the proposed confidence interval is found to be superior to Dixon's interval in terms of the coverage, the width, and stability. The advantages of the method do not appear to hold in the presence of nonnormal tolerance distribution.     

Survival estimation using splines

A nonparametric maximum likelihood procedure is given for estimating the survivor function from right-censored data. It approximates the hazard rate by a simple function such as a spline, with different approximations yielding different estimators. A special case is that proposed by Nelson (1969, Journal of Quality Technology 1, 27-52) and Altshuler (1970, Mathematical Biosciences 6, 1-11). The estimators are uniformly consistent and have the same asymptotic weak convergence properties as the Kaplan-Meier (1958, Journal of the American Statistical Association 53, 457-481) estimator. However, in small and in heavily censored samples, the simplest spline estimators have uniformly smaller mean squared error than do the Kaplan-Meier and Nelson-Altshuler estimators. The procedure is extended to estimate the baseline hazard rate and regression coefficients in the Cox (1972, Journal of the Royal Statistical Society, Series B 34, 187-220) proportional hazards model and is illustrated using experimental carcinogenesis data.