The propensity score method is widely used to estimate the average treatment effect in observational studies. However, it is generally confined to binary treatment assignment. In an extension to the settings of a multi-level treatment, Imbens proposed a generalized propensity score which is the conditional probability of receiving a particular level of the treatment given pre-treatment variables. The average treatment effect can then be estimated by conditioning solely on the generalized propensity score under the assumption of weak unconfoundedness. In the present work, we adopted this approach and conducted extensive simulations to evaluate the performance of several methods using the generalized propensity score, including subclassification, matching, inverse probability of treatment weighting, and covariate adjustment. Compared with other methods, inverse probability of treatment weighting had the preferred overall performance. We then applied these methods to a retrospective cohort study of 200,722 pregnant women. The impact of exposure to different types of antidepressant medications (no exposure, SSRI only, non-SSRI only, and both) during pregnancy on two important infant outcomes (birth weight and respiratory distress) were assessed.