MZINBVA: variational approximation for multilevel zero-inflated negative-binomial models for association analysis in microbiome surveys.

As our understanding of the microbiome has expanded, so has the recognition of its critical role in human health and disease, thereby emphasizing the importance of testing whether microbes are associated with environmental factors or clinical outcomes. However, many of the fundamental challenges that concern microbiome surveys arise from statistical and experimental design issues, such as the sparse and overdispersed nature of microbiome count data and the complex correlation structure among samples. For example, in the human microbiome project (HMP) dataset, the repeated observations across time points (level 1) are nested within body sites (level 2), which are further nested within subjects (level 3). Therefore, there is a great need for the development of specialized and sophisticated statistical tests. In this paper, we propose multilevel zero-inflated negative-binomial models for association analysis in microbiome surveys. We develop a variational approximation method for maximum likelihood estimation and inference. It uses optimization, rather than sampling, to approximate the log-likelihood and compute parameter estimates, provides a robust estimate of the covariance of parameter estimates and constructs a Wald-type test statistic for association testing. We evaluate and demonstrate the performance of our method using extensive simulation studies and an application to the HMP dataset. We have developed an R package MZINBVA to implement the proposed method, which is available from the GitHub repository https://github.com/liudoubletian/MZINBVA.

On variance estimation of the inverse probability-of-treatment weighting estimator: A tutorial for different types of propensity score weights.

Propensity score methods, such as inverse probability-of-treatment weighting (IPTW), have been increasingly used for covariate balancing in both observational studies and randomized trials, allowing the control of both systematic and chance imbalances. Approaches using IPTW are based on two steps: (i) estimation of the individual propensity scores (PS), and (ii) estimation of the treatment effect by applying PS weights. Thus, a variance estimator that accounts for both steps is crucial for correct inference. Using a variance estimator which ignores the first step leads to overestimated variance when the estimand is the average treatment effect (ATE), and to under or overestimated estimates when targeting the average treatment effect on the treated (ATT). In this article, we emphasize the importance of using an IPTW variance estimator that correctly considers the uncertainty in PS estimation. We present a comprehensive tutorial to obtain unbiased variance estimates, by proposing and applying a unifying formula for different types of PS weights (ATE, ATT, matching and overlap weights). This can be derived either via the linearization approach or M-estimation. Extensive R code is provided along with the corresponding large-sample theory. We perform simulation studies to illustrate the behavior of the estimators under different treatment and outcome prevalences and demonstrate appropriate behavior of the analytical variance estimator. We also use a reproducible analysis of observational lung cancer data as an illustrative example, estimating the effect of receiving a PET-CT scan on the receipt of surgery.

Designing individually randomized group treatment trials with repeated outcome measurements using generalized estimating equations.

Individually randomized group treatment (IRGT) trials, in which the clustering of outcome is induced by group-based treatment delivery, are increasingly popular in public health research. IRGT trials frequently incorporate longitudinal measurements, of which the proper sample size calculations should account for correlation structures reflecting both the treatment-induced clustering and repeated outcome measurements. Given the relatively sparse literature on designing longitudinal IRGT trials, we propose sample size procedures for continuous and binary outcomes based on the generalized estimating equations approach, employing the block exchangeable correlation structures with different correlation parameters for the treatment arm and for the control arm, and surveying five marginal mean models with different assumptions of time effect: no-time constant treatment effect, linear-time constant treatment effect, categorical-time constant treatment effect, linear time by treatment interaction, and categorical time by treatment interaction. Closed-form sample size formulas are derived for continuous outcomes, which depends on the eigenvalues of the correlation matrices; detailed numerical sample size procedures are proposed for binary outcomes. Through simulations, we demonstrate that the empirical power agrees well with the predicted power, for as few as eight groups formed in the treatment arm, when data are analyzed using the matrix-adjusted estimating equations for the correlation parameters with a bias-corrected sandwich variance estimator.

Leveraging baseline covariates to analyze small cluster-randomized trials with a rare binary outcome.

Cluster-randomized trials (CRTs) involve randomizing entire groups of participants-called clusters-to treatment arms but are often comprised of a limited or fixed number of available clusters. While covariate adjustment can account for chance imbalances between treatment arms and increase statistical efficiency in individually randomized trials, analytical methods for individual-level covariate adjustment in small CRTs have received little attention to date. In this paper, we systematically investigate, through extensive simulations, the operating characteristics of propensity score weighting and multivariable regression as two individual-level covariate adjustment strategies for estimating the participant-average causal effect in small CRTs with a rare binary outcome and identify scenarios where each adjustment strategy has a relative efficiency advantage over the other to make practical recommendations. We also examine the finite-sample performance of the bias-corrected sandwich variance estimators associated with propensity score weighting and multivariable regression for quantifying the uncertainty in estimating the participant-average treatment effect. To illustrate the methods for individual-level covariate adjustment, we reanalyze a recent CRT testing a sedation protocol in 31 pediatric intensive care units.

On Variance of the Treatment Effect in the Treated When Estimated by Inverse Probability Weighting.

In the analysis of observational studies, inverse probability weighting (IPW) is commonly used to consistently estimate the average treatment effect (ATE) or the average treatment effect in the treated (ATT). The variance of the IPW ATE estimator is often estimated by assuming that the weights are known and then using the so-called "robust" (Huber-White) sandwich estimator, which results in conservative standard errors (SEs). Here we show that using such an approach when estimating the variance of the IPW ATT estimator does not necessarily result in conservative SE estimates. That is, assuming the weights are known, the robust sandwich estimator may be either conservative or anticonservative. Thus, confidence intervals for the ATT using the robust SE estimate will not be valid, in general. Instead, stacked estimating equations which account for the weight estimation can be used to compute a consistent, closed-form variance estimator for the IPW ATT estimator. The 2 variance estimators are compared via simulation studies and in a data analysis of the association between smoking and gene expression.

Finite-sample adjustments in variance estimators for clustered competing risks regression.

The marginal Fine-Gray proportional subdistribution hazards model is a popular approach to directly study the association between covariates and the cumulative incidence function with clustered competing risks data, which often arise in multicenter randomized trials or multilevel observational studies. To account for the within-cluster correlations between failure times, the uncertainty of the regression parameters estimators is quantified by the robust sandwich variance estimator, which may have unsatisfactory performance with a limited number of clusters. To overcome this limitation, we propose four bias-corrected variance estimators to reduce the negative bias of the usual sandwich variance estimator, extending the bias-correction techniques from generalized estimating equations with noncensored exponential family outcomes to clustered competing risks outcomes. We further compare their finite-sample operating characteristics through simulations and two real data examples. In particular, we found the Mancl and DeRouen (MD) type sandwich variance estimator generally has the smallest bias. Furthermore, with a small number of clusters, the Wald t -confidence interval with the MD sandwich variance estimator carries close to nominal coverage for the cluster-level effect parameter. The t -confidence intervals based on the sandwich variance estimator with any one of the three types of multiplicative bias correction or the z -confidence interval with the Morel, Bokossa and Neerchal (MBN) type sandwich variance estimator have close to nominal coverage for the individual-level effect parameter. Finally, we develop a user-friendly R package crrcbcv implementing the proposed sandwich variance estimators to assist practical applications.

Variance estimators for weighted and stratified linear dose-response function estimators using generalized propensity score.

Propensity score methods are widely used in observational studies for evaluating marginal treatment effects. The generalized propensity score (GPS) is an extension of the propensity score framework, historically developed in the case of binary exposures, for use with quantitative or continuous exposures. In this paper, we proposed variance estimators for treatment effect estimators on continuous outcomes. Dose-response functions (DRFs) were estimated through weighting on the inverse of the GPS, or using stratification. Variance estimators were evaluated using Monte Carlo simulations. Despite the use of stabilized weights, the variability of the weighted estimator of the DRF was particularly high, and none of the variance estimators (a bootstrap-based estimator, a closed-form estimator especially developed to take into account the estimation step of the GPS, and a sandwich estimator) were able to adequately capture this variability, resulting in coverages below the nominal value, particularly when the proportion of the variation in the quantitative exposure explained by the covariates was large. The stratified estimator was more stable, and variance estimators (a bootstrap-based estimator, a pooled linearized estimator, and a pooled model-based estimator) more efficient at capturing the empirical variability of the parameters of the DRF. The pooled variance estimators tended to overestimate the variance, whereas the bootstrap estimator, which intrinsically takes into account the estimation step of the GPS, resulted in correct variance estimations and coverage rates. These methods were applied to a real data set with the aim of assessing the effect of maternal body mass index on newborn birth weight.

Variance estimation in inverse probability weighted Cox models.

Inverse probability weighted Cox models can be used to estimate marginal hazard ratios under different point treatments in observational studies. To obtain variance estimates, the robust sandwich variance estimator is often recommended to account for the induced correlation among weighted observations. However, this estimator does not incorporate the uncertainty in estimating the weights and tends to overestimate the variance, leading to inefficient inference. Here we propose a new variance estimator that combines the estimation procedures for the hazard ratio and weights using stacked estimating equations, with additional adjustments for the sum of terms that are not independently and identically distributed in a Cox partial likelihood score equation. We prove analytically that the robust sandwich variance estimator is conservative and establish the asymptotic equivalence between the proposed variance estimator and one obtained through linearization by Hajage et al. in 2018. In addition, we extend our proposed variance estimator to accommodate clustered data. We compare the finite sample performance of the proposed method with alternative methods through simulation studies. We illustrate these different variance methods in both independent and clustered data settings, using a bariatric surgery dataset and a multiple readmission dataset, respectively. To facilitate implementation of the proposed method, we have developed an R package ipwCoxCSV.

Analysis of recurrent hypoglycemic events.

Hypoglycemia is a major safety concern for diabetic patients. Hypoglycemic events can be modeled based on time to recurrent events or count data. In this article, we evaluated a gamma frailty model with variance estimated by the inverse of observed Fisher information matrix, a gamma frailty model with the sandwich variance estimator, and a piecewise negative binomial regression model. Simulations showed that the sandwich variance estimator performed better when the frailty model is mis-specified, and the piecewise negative binomial regression sometimes fails to converge. All three methods were applied to a dataset from a clinical trial evaluating insulin treatments.

Interaction analysis under misspecification of main effects: Some common mistakes and simple solutions.

The statistical practice of modeling interaction with two linear main effects and a product term is ubiquitous in the statistical and epidemiological literature. Most data modelers are aware that the misspecification of main effects can potentially cause severe type I error inflation in tests for interactions, leading to spurious detection of interactions. However, modeling practice has not changed. In this article, we focus on the specific situation where the main effects in the model are misspecified as linear terms and characterize its impact on common tests for statistical interaction. We then propose some simple alternatives that fix the issue of potential type I error inflation in testing interaction due to main effect misspecification. We show that when using the sandwich variance estimator for a linear regression model with a quantitative outcome and two independent factors, both the Wald and score tests asymptotically maintain the correct type I error rate. However, if the independence assumption does not hold or the outcome is binary, using the sandwich estimator does not fix the problem. We further demonstrate that flexibly modeling the main effect under a generalized additive model can largely reduce or often remove bias in the estimates and maintain the correct type I error rate for both quantitative and binary outcomes regardless of the independence assumption. We show, under the independence assumption and for a continuous outcome, overfitting and flexibly modeling the main effects does not lead to power loss asymptotically relative to a correctly specified main effect model. Our simulation study further demonstrates the empirical fact that using flexible models for the main effects does not result in a significant loss of power for testing interaction in general. Our results provide an improved understanding of the strengths and limitations for tests of interaction in the presence of main effect misspecification. Using data from a large biobank study "The Michigan Genomics Initiative", we present two examples of interaction analysis in support of our results.

The Calculus of M-Estimation in R with geex.

M-estimation, or estimating equation, methods are widely applicable for point estimation and asymptotic inference. In this paper, we present an R package that can find roots and compute the empirical sandwich variance estimator for any set of user-specified, unbiased estimating equations. Examples from the M-estimation primer by Stefanski and Boos (2002) demonstrate use of the software. The package also includes a framework for finite sample, heteroscedastic, and autocorrelation variance corrections, and a website with an extensive collection of tutorials.

Multilevel meta-analysis of multiple regression coefficients from single-case experimental studies.

The focus of the current study is on handling the dependence among multiple regression coefficients representing the treatment effects when meta-analyzing data from single-case experimental studies. We compare the results when applying three different multilevel meta-analytic models (i.e., a univariate multilevel model avoiding the dependence, a multivariate multilevel model ignoring covariance at higher levels, and a multivariate multilevel model modeling the existing covariance) to deal with the dependent effect sizes. The results indicate better estimates of the overall treatment effects and variance components when a multivariate multilevel model is applied, independent of modeling or ignoring the existing covariance. These findings confirm the robustness of multilevel modeling to misspecifying the existing covariance at the case and study level in terms of estimating the overall treatment effects and variance components. The results also show that the overall treatment effect estimates are unbiased regardless of the underlying model, but the between-case and between-study variance components are biased in certain conditions. In addition, the between-study variance estimates are particularly biased when the number of studies is smaller than 40 (i.e., 10 or 20) and the true value of the between-case variance is relatively large (i.e., 8). The observed bias is larger for the between-case variance estimates compared to the between-study variance estimates when the true between-case variance is relatively small (i.e., 0.5).

A critical issue of using the variance of the total in the linearization method - In the context of unequal probability sampling.

Publicly available national survey data are useful for the evidence-based research to advance our understanding of important questions in the health and biomedical sciences. Appropriate variance estimation is a crucial step to evaluate the strength of evidence in the data analysis. In survey data analysis, the conventional linearization method for estimating the variance of a statistic of interest uses the variance estimator of the total based on linearized variables. We warn that this common practice may result in undesirable consequences such as susceptibility to data shift and severely inflated variance estimates, when unequal weights are incorporated into variance estimation. We propose to use the variance estimator of the mean (mean-approach) instead of the variance estimator of the total (total-approach). We show a superiority of the mean-approach through analytical investigations. A real data example (the National Comorbidity Survey Replication) and simulation-based studies strongly support our conclusion.

A note on the empirical Bayes heterogeneity variance estimator in meta-analysis.

This paper focuses on the empirical Bayes (EB) or Mandel-Paule estimator of the heterogeneity variance in meta-analysis, which was discussed by Morris and proposed in earlier publications by Mandel and Paule in an inter-laboratory context. The relationship of the EB estimator to other heterogeneity variance estimators typically used in meta-analysis is explored, and approximate variance estimators for the EB estimate of the heterogeneity variance are proposed based on the M-estimation method. Statistical inference for the overall treatment effect using the EB estimator and the proposed standard errors is discussed using two example data sets from meta-analysis applications.

Power and sample size requirements for GEE analyses of cluster randomized crossover trials.

The cluster randomized crossover design has been proposed to improve efficiency over the traditional parallel cluster randomized design, which often involves a limited number of clusters. In recent years, the cluster randomized crossover design has been increasingly used to evaluate the effectiveness of health care policy or programs, and the interest often lies in quantifying the population-averaged intervention effect. In this paper, we consider the two-treatment two-period crossover design, and develop sample size procedures for continuous and binary outcomes corresponding to a population-averaged model estimated by generalized estimating equations, accounting for both within-period and interperiod correlations. In particular, we show that the required sample size depends on the correlation parameters through an eigenvalue of the within-cluster correlation matrix for continuous outcomes and through two distinct eigenvalues of the correlation matrix for binary outcomes. We demonstrate that the empirical power corresponds well with the predicted power by the proposed formulae for as few as eight clusters, when outcomes are analyzed using the matrix-adjusted estimating equations for the correlation parameters concurrently with a suitable bias-corrected sandwich variance estimator.

Power and sample size calculation incorporating misspecifications of the variance function in comparative clinical trials with over-dispersed count data.

Over-dispersed count data are frequently observed in clinical trials where the primary endpoint is occurrence of clinical events. Sample sizes of comparative clinical trials with these data are typically calculated under negative binomial models or quasi-Poisson models with specified variance functions, or under the assumption that the specified "working" variance functions are correctly specified. In this article, we propose a sample size formula anticipating misspecifications of the working variance function. We derived a method based on the asymptotic distribution of a Wald test statistic with a sandwich-type robust variance estimator under quasi-Poisson models. Under misspecifications of the working variance function, the asymptotic variance of the estimator of the treatment effect is expressed as a form involving both true and working variance functions. Our sample size formula includes several existing formulas as special cases when the working variance function is correctly specified as the true variance function. We also consider a sensitivity analysis for possible misspecifications of the "true" variance function when estimating sample sizes using our formula. A simulation study demonstrated the adequacy of our formulas in finite sample size settings. An application to a clinical trial to evaluate the treatment effect on prevention of COPD exacerbation is provided.

Closed-form variance estimator for weighted propensity score estimators with survival outcome.

Propensity score (PS) methods are widely used in observational studies for evaluating marginal treatment effects. PS-weighting is a popular PS-based method that allows for estimating both the average treatment effect on the overall population (ATE) and the average treatment effect on the treated population (ATT). Previous research has shown that the variance of the treatment effect is accurately estimated only if the variance estimator takes into account the fact that the propensity score is itself estimated from the available data in a first step of the analysis. In 2016, Austin showed that the bootstrap-based variance estimator was the only existing estimator resulting in approximately correct estimates of standard errors when evaluating a survival outcome and a Cox model was used to estimate a marginal hazard ratio (HR). This author stressed the need to develop a closed-form variance estimator of the marginal HR accounting for the estimation of the PS. In the present research, we developed such variance estimators both for the ATE and ATT. We evaluated their performance with an extensive simulation study and compared them to bootstrap-based variance estimators and to naive variance estimators that do not account for the estimation step. We found that the performance of the proposed variance estimators was similar to that of the bootstrap-based estimators. The proposed variance estimators provide an alternative to the bootstrap estimator, particularly interesting in situations in which time-consumption and/or reproducibility are an important issue. An implementation has been developed for the R software and is freely available (package hrIPW).

On the multiple imputation variance estimator for control-based and delta-adjusted pattern mixture models.

Control-based pattern mixture models (PMM) and delta-adjusted PMMs are commonly used as sensitivity analyses in clinical trials with non-ignorable dropout. These PMMs assume that the statistical behavior of outcomes varies by pattern in the experimental arm in the imputation procedure, but the imputed data are typically analyzed by a standard method such as the primary analysis model. In the multiple imputation (MI) inference, Rubin's variance estimator is generally biased when the imputation and analysis models are uncongenial. One objective of the article is to quantify the bias of Rubin's variance estimator in the control-based and delta-adjusted PMMs for longitudinal continuous outcomes. These PMMs assume the same observed data distribution as the mixed effects model for repeated measures (MMRM). We derive analytic expressions for the MI treatment effect estimator and the associated Rubin's variance in these PMMs and MMRM as functions of the maximum likelihood estimator from the MMRM analysis and the observed proportion of subjects in each dropout pattern when the number of imputations is infinite. The asymptotic bias is generally small or negligible in the delta-adjusted PMM, but can be sizable in the control-based PMM. This indicates that the inference based on Rubin's rule is approximately valid in the delta-adjusted PMM. A simple variance estimator is proposed to ensure asymptotically valid MI inferences in these PMMs, and compared with the bootstrap variance. The proposed method is illustrated by the analysis of an antidepressant trial, and its performance is further evaluated via a simulation study.

Closed-form REML estimators and sample size determination for mixed effects models for repeated measures under monotone missingness.

We derive the closed-form restricted maximum likelihood estimator and Kenward-Roger's variance estimator for fixed effects in the mixed effects model for repeated measures (MMRM) when the missing data pattern is monotone. As an important application of the analytic result, we present the formula for calculating the power of treatment comparison using the Wald t-test with the Kenward-Roger adjusted variance estimate in MMRM. It allows adjustment for baseline covariates without the need to specify the covariate distribution in randomized trials. A simple two-step procedure is proposed to determine the sample size needed to achieve the targeted power. The proposed method performs well for both normal and moderately non-normal data even in small samples (n=20) in simulations. An antidepressant trial is analyzed for illustrative purposes. Copyright © 2017 John Wiley & Sons, Ltd.

Covariance estimators for generalized estimating equations (GEE) in longitudinal analysis with small samples.

Generalized estimating equations (GEE) is a general statistical method to fit marginal models for longitudinal data in biomedical studies. The variance-covariance matrix of the regression parameter coefficients is usually estimated by a robust "sandwich" variance estimator, which does not perform satisfactorily when the sample size is small. To reduce the downward bias and improve the efficiency, several modified variance estimators have been proposed for bias-correction or efficiency improvement. In this paper, we provide a comprehensive review on recent developments of modified variance estimators and compare their small-sample performance theoretically and numerically through simulation and real data examples. In particular, Wald tests and t-tests based on different variance estimators are used for hypothesis testing, and the guideline on appropriate sample sizes for each estimator is provided for preserving type I error in general cases based on numerical results. Moreover, we develop a user-friendly R package "geesmv" incorporating all of these variance estimators for public usage in practice.