Evaluating the impact of non-randomized treatment on various health outcomes is difficult in observational studies because of the presence of covariates that may affect both the treatment or exposure received and the outcome of interest. In the present study, we develop a semiparametric multiply robust multiple imputation method for estimating average treatment effects in such studies. Our method combines information from multiple propensity score models and outcome regression models, and is multiply robust in that it produces consistent estimators for the average causal effects if at least one of the models is correctly specified. Our proposed estimators show promising performances even with incorrect models. Compared with existing fully parametric approaches, our proposed method is more robust against model misspecifications. Compared with fully non-parametric approaches, our proposed method does not have the problem of curse of dimensionality and achieves dimension reduction by combining information from multiple models. In addition, it is less sensitive to the extreme propensity score estimates compared with inverse propensity score weighted estimators and augmented estimators. The asymptotic properties of our method are developed and the simulation study shows the advantages of our proposed method compared with some existing methods in terms of balancing efficiency, bias, and coverage probability. Rubin's variance estimation formula can be used for estimating the variance of our proposed estimators. Finally, we apply our method to 2009-2010 National Health Nutrition and Examination Survey (NHANES) to examine the effect of exposure to perfluoroalkyl acids (PFAs) on kidney function.
