Classification with positive and unlabeled (PU) data frequently arises in bioinformatics, clinical data, and ecological studies, where collecting negative samples can be prohibitively expensive. While prior works on PU data focus on binary classification, in this paper we consider multiple positive labels, a practically important and common setting. We introduce a multinomial-PU model and an ordinal-PU model, suited to unordered and ordered labels respectively. We propose proximal gradient descent-based algorithms to minimize the l_{1,2}-penalized log-likelihood losses, with convergence guarantees to stationary points of the non-convex objective. Despite the challenging non-convexity induced by the presence-only data and multi-class labels, we prove statistical error bounds for the stationary points within a neighborhood around the true parameters under the high-dimensional regime. This is made possible through a careful characterization of the landscape of the log-likelihood loss in the neighborhood. In addition, simulations and two real data experiments demonstrate the empirical benefits of our algorithms compared to the baseline methods.