There is extensive literature on methods, such as propensity scoring, for estimating the causal effects for two treatments using real world data. Much less work has been done for the more general setting with three or more treatments.  Whereas the literature has suggested that these propensity-based methods do not naturally extend to the multi-level treatment case, we show, using the concept of weak unconfoundedness, that adjusting for or matching on a scalar function of the covariates removes biases associated with observed covariates. We focused on subclassification and matching approaches as these have found to be effective for two treatments and are among the most popular methods in that setting. We apply the proposed methods to an analysis of the effectiveness of treatments for fibromyalgia from a prospective observational study. We also carried out a simulation study to assess the performance of those new methods relative to such approaches like: pairwise propensity score matching; matching on the Mahalanobis distance of all covariates; matching on the set of propensity scores (with the number of scores equal to the number of distinct treatment levels minus one (Rassen, 2013)); weighting on the inverse of the binary treatment propensity scores (McCaffrey, 2013). The simulations suggest that the proposed methods are simple and viable options for comparing the effectiveness of three or more treatments. RASSEN et al.: Matching by propensity score in cohort studies with three treatment groups. Epidemiology 24, 401–9. MCCAFFREY et al.: A tutorial on propensity score estimation for multiple treatments using generalized boosted models. Stat. Med. 32, 3388–414.