Generalized estimating equations: A pragmatic and flexible approach to the marginal GLM modelling of correlated data in the behavioural sciences

Within behavioural research, non‐normally distributed data with a complicated structure are common. For instance, data can represent repeated observations of quantities on the same individual. The regression analysis of such data is complicated both by the interdependency of the observations (response variables) and by their non‐normal distribution. Over the last decade, such data have been more and more frequently analysed using generalized mixed‐effect models. Some researchers invoke the heavy machinery of mixed‐effect modelling to obtain the desired population‐level (marginal) inference, which can be achieved by using simpler tools—namely by marginal models. This paper highlights marginal modelling (using generalized estimating equations [GEE]) as an alternative method. In various situations, GEE can be based on fewer assumptions and directly generate estimates (population‐level parameters) which are of immediate interest to the behavioural researcher (such as population means). Using four examples from behavioural research, we demonstrate the use, advantages, and limits of the GEE approach as implemented within the functions of the ‘geepack’ package in R.     

Marginal hazard ratio estimates in joint frailty models for heart failure trials

This work is motivated by clinical trials in chronic heart failure disease, where treatment has effects both on morbidity (assessed as recurrent non‐fatal hospitalisations) and on mortality (assessed as cardiovascular death, CV death). Recently, a joint frailty proportional hazards model has been proposed for these kind of efficacy outcomes to account for a potential association between the risk rates for hospital admissions and CV death. However, more often clinical trial results are presented by treatment effect estimates that have been derived from marginal proportional hazards models, that is, a Cox model for mortality and an Andersen–Gill model for recurrent hospitalisations. We show how these marginal hazard ratios and their estimates depend on the association between the risk processes, when these are actually linked by shared or dependent frailty terms. First we derive the marginal hazard ratios as a function of time. Then, applying least false parameter theory, we show that the marginal hazard ratio estimate for the hospitalisation rate depends on study duration and on parameters of the underlying joint frailty model. In particular, we identify parameters, for example the treatment effect on mortality, that determine if the marginal hazard ratio estimate for hospitalisations is smaller, equal or larger than the conditional one. How this affects rejection probabilities is further investigated in simulation studies. Our findings can be used to interpret marginal hazard ratio estimates in heart failure trials and are illustrated by the results of the CHARM‐Preserved trial (where CHARM is the ‘Candesartan in Heart failure Assessment of Reduction in Mortality and morbidity’ programme).     

Unbiased Estimation of Odds Ratio by Using Negative Binomial Plans

The use of the negative binomial distribution in both the numerator and denominator in prospective studies leads to an unbiased estimate of the odds ratio and an exact expression for its variance. Sample sizes that minimize the variance of odds ratio estimates are specified. The variance of the odds ratio estimate is shown to be close to the Cramér-Rao lower bound.     

Estimates of Survival from Radiotelemetry: a Response to Guthery and Lusk

Abstract: Radiotelemetry has become an important and frequently used tool in wildlife research. Inferences drawn from radiotelemetry data depend on the assumption that the radiotransmitters are not influencing parameter(s) of interest. An article by Guthery and Lusk (2004) in the Wildlife Society Bulletin questioned the validity of this assumption for estimating survival rates of northern bobwhites (Colinus virginianus) using radiotelemetry data. In this evaluation, we address technical and philosophical flaws in Guthery and Lusk's (2004) critique of northern bobwhite studies utilizing radiotelemetry. They concluded that biologists should be skeptical of radiotelemetry studies and they advised researchers to design studies to address potential biases caused by radiotransmitters using independent data. Although we agree that researchers are responsible for testing key assumptions of their techniques, we believe Guthery and Lusk's (2004) conclusions were not well supported and were based on tenuous assumptions. Guthery and Lusk (2004) calculated the level of productivity (given as a fall age ratio) required to balance a simple population model that contained published estimates of annual survival and assumed an annual finite population growth rate of 1.0. We review their population model and show that the relationship between an annual survival rate and fall age ratio is nonlinear. This nonlinearity can lead to biased estimates of a fall age ratio, especially at lower values of annual survival. We also question the validity of using fall age ratios as an estimator of productivity. Further, we suggest that this assessment of a radiotransmitter effect from a survival rate itself is not appropriate. This rate can be depressed (or elevated) for a variety of reasons not related to the influence of radiotransmitters. In addition, Guthery and Lusk (2004) assumed that daily survival rates (as calculated from both annual and seasonal published estimates) were constant throughout the year; thus, they scaled daily survival rates from seasonal to annual estimates. Further, their meta-analysis was hindered by temporal pseudoreplication and a lack of independence among the observations used in the analysis. We conclude the weight of the evidence presented by Guthery and Lusk (2004) is not as strong as they claim because it fails to meet the test of sufficient causation. While scientists should always be skeptical and critical of assumptions of all methods employed in wildlife research, more rigorous tests are necessary before we discredit a valuable technique without sufficient empirical evidence.     

Latent Class Modeling Approaches for Assessing Diagnostic Error without a Gold Standard: With Applications to p53 Immunohistochemical Assays in Bladder Tumors

Improved characterization of tumors for purposes of guiding treatment decisions for cancer patients will require that accurate and reproducible assays be developed for a variety of tumor markers. No gold standards exist for most tumor marker assays. Therefore, estimates of assay sensitivity and specificity cannot be obtained unless a latent class model-based approach is used. Our goal in this article is to estimate sensitivity and specificity for p53 immunohistochemical assays of bladder tumors using data from a reproducibility study conducted by the National Cancer Institute Bladder Tumor Marker Network. We review latent class modeling approaches proposed by previous authors, and we find that many of these approaches impose assumptions about specimen heterogeneity that are not consistent with the biology of bladder tumors. We present flexible mixture model alternatives that are biologically plausible for our example, and we use them to estimate sensitivity and specificity for our p53 assay example. These mixture models are shown to offer an improvement over other methods in a variety of settings, but we caution that, in general, care must be taken in applying latent class models.     

Skew‐normal/independent linear mixed models for censored responses with applications to HIV viral loads

Often in biomedical studies, the routine use of linear mixed‐effects models (based on Gaussian assumptions) can be questionable when the longitudinal responses are skewed in nature. Skew‐normal/elliptical models are widely used in those situations. Often, those skewed responses might also be subjected to some upper and lower quantification limits (QLs; viz., longitudinal viral‐load measures in HIV studies), beyond which they are not measurable. In this paper, we develop a Bayesian analysis of censored linear mixed models replacing the Gaussian assumptions with skew‐normal/independent (SNI) distributions. The SNI is an attractive class of asymmetric heavy‐tailed distributions that includes the skew‐normal, skew‐t, skew‐slash, and skew‐contaminated normal distributions as special cases. The proposed model provides flexibility in capturing the effects of skewness and heavy tail for responses that are either left‐ or right‐censored. For our analysis, we adopt a Bayesian framework and develop a Markov chain Monte Carlo algorithm to carry out the posterior analyses. The marginal likelihood is tractable, and utilized to compute not only some Bayesian model selection measures but also case‐deletion influence diagnostics based on the Kullback–Leibler divergence. The newly developed procedures are illustrated with a simulation study as well as an HIV case study involving analysis of longitudinal viral loads.     

Cutting edge: selection of B lymphocyte subsets is regulated by natural IgM

Natural IgM has a wide range of actions in the immune system. Here we demonstrate that mice lacking serum IgM have an expansion in splenic marginal zone B cells with a proportionately smaller reduction in follicular B cells. The increase in the marginal zone-follicular B cell ratio (and an expansion in peritoneal B1a cells) is fully reversed by administration of polyclonal IgM, but not by two IgM monoclonals. Mice engineered to have a secreted oligoclonal IgM repertoire with an endogenous membrane IgM also exhibited a similar expansion of marginal zone B cells. We propose that natural IgM, by virtue of its polyreactivity, enhances Ag-driven signaling through the B cell receptor and promotes the formation of follicular B cells. These results demonstrate that natural IgM regulates the selection of B lymphocyte subsets.     

A review of global potentially available cropland estimates and their consequences for model‐based assessments

The world's population is growing and demand for food, feed, fiber, and fuel is increasing, placing greater demand on land and its resources for crop production. We review previously published estimates of global scale cropland availability, discuss the underlying assumptions that lead to differences between estimates, and illustrate the consequences of applying different estimates in model‐based assessments of land‐use change. The review estimates a range from 1552 to 5131 Mha, which includes 1550 Mha that is already cropland. Hence, the lowest estimates indicate that there is almost no room for cropland expansion, while the highest estimates indicate that cropland could potentially expand to over three times its current area. Differences can largely be attributed to institutional assumptions, i.e. which land covers/uses (e.g. forests or grasslands) are societally or governmentally allowed to convert to cropland, while there was little variation in biophysical assumptions. Estimates based on comparable assumptions showed a variation of up to 84%, which originated mainly from different underlying data sources. On the basis of this synthesis of the assumptions underlying these estimates, we constructed a high, a medium, and a low estimate of cropland availability that are representative of the range of estimates in the reviewed studies. We apply these estimates in a land‐change model to illustrate the consequences on cropland expansion and intensification as well as deforestation. While uncertainty in cropland availability is hardly addressed in global land‐use change assessments, the results indicate a large range of estimates with important consequences for model‐based assessments.     

The aggregate income losses from childhood stunting and the returns to a nutrition intervention aimed at reducing stunting

We undertake two calculations, one for all developing countries, the other for 34 developing countries that together account for 90% of the world’s stunted children. The first asks how much lower a country’s per capita income is today as a result of having a fraction of its workforce been stunted in childhood. We use a development accounting framework, relying on micro-econometric estimates of the effects of childhood stunting on adult wages through their effects on years of schooling, cognitive skills, and height, parsing out the relative contribution of each set of returns to avoid double counting. We estimate that, on average, the per capita income penalty from stunting is between 5–7%, depending on the assumption. In our second calculation we estimate the economic value and the costs associates with scaling up a package of nutrition interventions using the same methodology and set of assumptions used in the first calculation. We take a package of 10 nutrition interventions that has data on both effects and costs, and we estimate the rate-of-return to gradually introducing this program over a period of 10 years in 34 countries that together account for 90% of the world’s stunted children. We estimate a rate-of-return of 12%, and a benefit-cost ratio of 5:1-6:1.     

Regression survival analysis with an assumed copula for dependent censoring: a sensitivity analysis approach

SUMMARY: In clinical studies, when censoring is caused by competing risks or patient withdrawal, there is always a concern about the validity of treatment effect estimates that are obtained under the assumption of independent censoring. Because dependent censoring is nonidentifiable without additional information, the best we can do is a sensitivity analysis to assess the changes of parameter estimates under different assumptions about the association between failure and censoring. This analysis is especially useful when knowledge about such association is available through literature review or expert opinions. In a regression analysis setting, the consequences of falsely assuming independent censoring on parameter estimates are not clear. Neither the direction nor the magnitude of the potential bias can be easily predicted. We provide an approach to do sensitivity analysis for the widely used Cox proportional hazards models. The joint distribution of the failure and censoring times is assumed to be a function of their marginal distributions. This function is called a copula. Under this assumption, we propose an iteration algorithm to estimate the regression parameters and marginal survival functions. Simulation studies show that this algorithm works well. We apply the proposed sensitivity analysis approach to the data from an AIDS clinical trial in which 27% of the patients withdrew due to toxicity or at the request of the patient or investigator.     

Marginal Models Via GLS: A Convenient Yet Neglected Tool for the Analysis of Correlated Data in the Behavioural Sciences

Behavioural research often produces data that have a complicated structure. For instance, data can represent repeated observations of the same individual and suffer from heteroscedasticity as well as other technical snags. The regression analysis of such data is often complicated by the fact that the observations (response variables) are mutually correlated. The correlation structure can be quite complex and might or might not be of direct interest to the user. In any case, one needs to take correlations into account (e.g. by means of random‐effect specification) in order to arrive at correct statistical inference (e.g. for construction of the appropriate test or confidence intervals). Over the last decade, such data have been more and more frequently analysed using repeated‐measures ANOVA and mixed‐effects models. Some researchers invoke the heavy machinery of mixed‐effects modelling to obtain the desired population‐level (marginal) inference, which can be achieved by using simpler tools – namely marginal models. This paper highlights marginal modelling (using generalized least squares [GLS] regression) as an alternative method. In various concrete situations, such marginal models can be based on fewer assumptions and directly generate estimates (population‐level parameters) which are of immediate interest to the behavioural researcher (such as population mean). Sometimes, they might be not only easier to interpret but also easier to specify than their competitors (e.g. mixed‐effects models). Using five examples from behavioural research, we demonstrate the use, advantages, limits and pitfalls of marginal and mixed‐effects models implemented within the functions of the ‘nlme’ package in R.     

Adjusting for selection bias in retrospective, case-control studies

Retrospective case–control studies are more susceptibleto selection bias than other epidemiologic studies as by designthey require that both cases and controls are representativeof the same population. However, as cases and control recruitmentprocesses are often different, it is not always obvious thatthe necessary exchangeability conditions hold. Selection biastypically arises when the selection criteria are associatedwith the risk factor under investigation. We develop a methodwhich produces bias-adjusted estimates for the odds ratio. Ourmethod hinges on 2 conditions. The first is that a variablethat separates the risk factor from the selection criteria canbe identified. This is termed the "bias breaking" variable.The second condition is that data can be found such that a bias-correctedestimate of the distribution of the bias breaking variable canbe obtained. We show by means of a set of examples that suchbias breaking variables are not uncommon in epidemiologic settings.We demonstrate using simulations that the estimates of the oddsratios produced by our method are consistently closer to thetrue odds ratio than standard odds ratio estimates using logisticregression. Further, by applying it to a case–controlstudy, we show that our method can help to determine whetherselection bias is present and thus confirm the validity of studyconclusions when no evidence of selection bias can be found.     

Computing Bayes factors using thermodynamic integration

In the Bayesian paradigm, a common method for comparing two models is to compute the Bayes factor, defined as the ratio of their respective marginal likelihoods. In recent phylogenetic works, the numerical evaluation of marginal likelihoods has often been performed using the harmonic mean estimation procedure. In the present article, we propose to employ another method, based on an analogy with statistical physics, called thermodynamic integration. We describe the method, propose an implementation, and show on two analytical examples that this numerical method yields reliable estimates. In contrast, the harmonic mean estimator leads to a strong overestimation of the marginal likelihood, which is all the more pronounced as the model is higher dimensional. As a result, the harmonic mean estimator systematically favors more parameter-rich models, an artefact that might explain some recent puzzling observations, based on harmonic mean estimates, suggesting that Bayes factors tend to overscore complex models. Finally, we apply our method to the comparison of several alternative models of amino-acid replacement. We confirm our previous observations, indicating that modeling pattern heterogeneity across sites tends to yield better models than standard empirical matrices.     

Deletion diagnostics for alternating logistic regressions

Deletion diagnostics are introduced for the regression analysis of clustered binary outcomes estimated with alternating logistic regressions, an implementation of generalized estimating equations (GEE) that estimates regression coefficients in a marginal mean model and in a model for the intracluster association given by the log odds ratio. The diagnostics are developed within an estimating equations framework that recasts the estimating functions for association parameters based upon conditional residuals into equivalent functions based upon marginal residuals. Extensions of earlier work on GEE diagnostics follow directly, including computational formulae for one‐step deletion diagnostics that measure the influence of a cluster of observations on the estimated regression parameters and on the overall marginal mean or association model fit. The diagnostic formulae are evaluated with simulations studies and with an application concerning an assessment of factors associated with health maintenance visits in primary care medical practices. The application and the simulations demonstrate that the proposed cluster‐deletion diagnostics for alternating logistic regressions are good approximations of their exact fully iterated counterparts.     

Chaos in three species food chains

We study the dynamics of a three species food chain using bifurcation theory to demonstrate the existence of chaotic dynamics in the neighborhood of the equilibrium where the top species in the food chain is absent. The goal of our study is to demonstrate the presence of chaos in a class of ecological models, rather than just in a specific model. This work extends earlier numerical studies of a particular system by Hastings and Powell (1991) by showing that chaos occurs in a class of ecological models. The mathematical techniques we use are based on work by Guckenheimer and Holmes (1983) on co-dimension two bifurcations. However, restrictions on the equations we study imposed by ecological assumptions require a new and somewhat different analysis.     

Robust estimation of microbial diversity in theory and in practice

Quantifying diversity is of central importance for the study of structure, function and evolution of microbial communities. The estimation of microbial diversity has received renewed attention with the advent of large-scale metagenomic studies. Here, we consider what the diversity observed in a sample tells us about the diversity of the community being sampled. First, we argue that one cannot reliably estimate the absolute and relative number of microbial species present in a community without making unsupported assumptions about species abundance distributions. The reason for this is that sample data do not contain information about the number of rare species in the tail of species abundance distributions. We illustrate the difficulty in comparing species richness estimates by applying Chao''s estimator of species richness to a set of in silico communities: they are ranked incorrectly in the presence of large numbers of rare species. Next, we extend our analysis to a general family of diversity metrics (‘Hill diversities''), and construct lower and upper estimates of diversity values consistent with the sample data. The theory generalizes Chao''s estimator, which we retrieve as the lower estimate of species richness. We show that Shannon and Simpson diversity can be robustly estimated for the in silico communities. We analyze nine metagenomic data sets from a wide range of environments, and show that our findings are relevant for empirically-sampled communities. Hence, we recommend the use of Shannon and Simpson diversity rather than species richness in efforts to quantify and compare microbial diversity.     

Parentage versus two-generation analyses for estimating pollen-mediated gene flow in plant populations

Assessment of contemporary pollen-mediated gene flow in plants is important for various aspects of plant population biology, genetic conservation and breeding. Here, through simulations we compare the two alternative approaches for measuring pollen-mediated gene flow: (i) the NEIGHBORHOOD model--a representative of parentage analyses, and (ii) the recently developed TWOGENER analysis of pollen pool structure. We investigate their properties in estimating the effective number of pollen parents (N(ep)) and the mean pollen dispersal distance (delta). We demonstrate that both methods provide very congruent estimates of N(ep) and delta, when the methods' assumptions considering the shape of pollen dispersal curve and the mating system follow those used in data simulations, although the NEIGHBORHOOD model exhibits generally lower variances of the estimates. The violations of the assumptions, especially increased selfing or long-distance pollen dispersal, affect the two methods to a different degree; however, they are still capable to provide comparable estimates of N(ep). The NEIGHBORHOOD model inherently allows to estimate both self-fertilization and outcrossing due to the long-distance pollen dispersal; however, the TWOGENER method is particularly sensitive to inflated selfing levels, which in turn may confound and suppress the effects of distant pollen movement. As a solution we demonstrate that in case of TWOGENER it is possible to extract the fraction of intraclass correlation that results from outcrossing only, which seems to be very relevant for measuring pollen-mediated gene flow. The two approaches differ in estimation precision and experimental efforts but they seem to be complementary depending on the main research focus and type of a population studied.     

Marginal false discovery rate control for likelihood‐based penalized regression models

The popularity of penalized regression in high‐dimensional data analysis has led to a demand for new inferential tools for these models. False discovery rate control is widely used in high‐dimensional hypothesis testing, but has only recently been considered in the context of penalized regression. Almost all of this work, however, has focused on lasso‐penalized linear regression. In this paper, we derive a general method for controlling the marginal false discovery rate that can be applied to any penalized likelihood‐based model, such as logistic regression and Cox regression. Our approach is fast, flexible and can be used with a variety of penalty functions including lasso, elastic net, MCP, and MNet. We derive theoretical results under which the proposed method is valid, and use simulation studies to demonstrate that the approach is reasonably robust, albeit slightly conservative, when these assumptions are violated. Despite being conservative, we show that our method often offers more power to select causally important features than existing approaches. Finally, the practical utility of the method is demonstrated on gene expression datasets with binary and time‐to‐event outcomes.     

Estimating treatment effects with partially observed covariates using outcome regression with missing indicators

Missing data is a common issue in research using observational studies to investigate the effect of treatments on health outcomes. When missingness occurs only in the covariates, a simple approach is to use missing indicators to handle the partially observed covariates. The missing indicator approach has been criticized for giving biased results in outcome regression. However, recent papers have suggested that the missing indicator approach can provide unbiased results in propensity score analysis under certain assumptions. We consider assumptions under which the missing indicator approach can provide valid inferences, namely, (1) no unmeasured confounding within missingness patterns; either (2a) covariate values of patients with missing data were conditionally independent of treatment or (2b) these values were conditionally independent of outcome; and (3) the outcome model is correctly specified: specifically, the true outcome model does not include interactions between missing indicators and fully observed covariates. We prove that, under the assumptions above, the missing indicator approach with outcome regression can provide unbiased estimates of the average treatment effect. We use a simulation study to investigate the extent of bias in estimates of the treatment effect when the assumptions are violated and we illustrate our findings using data from electronic health records. In conclusion, the missing indicator approach can provide valid inferences for outcome regression, but the plausibility of its assumptions must first be considered carefully.