We consider analysis of sparsely sampled multilevel functional data, where the basic observational unit is a function and data have a natural hierarchy of basic units. An example is when functions are recorded at multiple visits for each subject. Multilevel functional principal component analysis was proposed recently for such data when functions are densely recorded. Here, we consider the case when functions are sparsely sampled and may contain only a few observations per function. We exploit the multilevel structure of covariance operators and achieve data reduction by principal component decompositions at both between‐subject and within‐subject levels. We address inherent methodological differences in the sparse sampling context to: (i) estimate the covariance operators; (ii) estimate the functional principal component scores; and (iii) predict the underlying curves. Through simulations, the proposed method is able to discover dominating modes of variations and reconstruct underlying curves well even in sparse settings. Our approach is illustrated by two applications, the Sleep Heart Health Study and eBay auctions. Copyright © 2014 John Wiley & Sons, Ltd