In this talk, we consider an elliptic partial differential equation with random diffusion coefficient, defined on a bounded two-dimensional domain with homogeneous Dirichlet boundary conditions. We are interested in deriving reliable approximations for the expected value of the solution and of expected values of linear quantities of interest, defined as weighted integrals of the solution over the domain. The deterministic PDE obtained for a fixed realization of the random diffusion coefficient is discretized by means of the conforming hp Virtual Element (VE) method. Novel error estimates for the VE approximation of quantities of interest. To discretize the integral in probability, we employ the Multi Level Monte Carlo (MLMC) method that considers multiple levels of spatial resolution: the estimator is first computed on a coarse resolution level, and successive correction terms are then added to improve accuracy. First, we analyze the h-version of the MLMC-VE method, which relies on a sequence of nested polygonal meshes of the domain. In practice, such sequence is constructed using a mesh-agglomeration strategy. Computable estimators for both the expected value of the solution and the expected value of quantities of interest are designed, and error estimates are derived. Then, a p-version of the MLMC-VE method is studied. Numerical tests that verify the theoretical error estimates and validate the proposed methodology will be presented.