We quantify the effect of uncertainties on quantities of interest related to contact mechanics of rough surfaces. Specifically, we consider the problem of frictionless non adhesive normal contact between two semi infinite linear elastic solids subject to uncertainties. These uncertainties may for example originate from an incomplete surface description. To account for surface uncertainties, we model a rough surface as a suitable Gaussian random field whose covariance function encodes the surface's roughness, which is experimentally measurable. Within this stochastic framework, we first introduce the complete random contact model, which includes the precise definition of the considered class of rough random surfaces as well as the study of a practical random surface generator. Then, we introduce the multilevel Monte Carlo method which is a computationally efficient sampling method for the computation of statistical moments of uncertain system output's, such as those derived from contact simulations. In particular, we consider two different quantities of interest, namely the contact area and the number of contact clusters, and show via numerical experiments that the multilevel Monte Carlo method offers significant computational gains compared to an approximation via a classic Monte Carlo sampling.