We develop a novel method for counterfactual analysis based on observational data using prediction intervals for units under different exposures. Unlike methods that target heterogeneous or conditional average treatment effects of an exposure, the proposed approach aims to take into account the irreducible dispersions of counterfactual outcomes so as to quantify the relative impact of different exposures. The prediction intervals are constructed in a distribution-free and model-robust manner based on the conformal prediction approach. The computational obstacles to this approach are circumvented by leveraging properties of a tuning-free method that learns sparse additive predictor models for counterfactual outcomes. The method is illustrated using both real and synthetic data.