We present an approach to improve the convergence of a multilevel Monte Carlo method using a high-order discretization scheme for elliptic partial differential equations with random coefficients. We describe in detail a multilevel estimator which utilizes a fourth-order accurate solution of the elliptic PDEs. The superiority of this fourth-order estimator is reflected in terms of fewer multilevel Monte Carlo levels and samples compared to the existing method. Numerical experiments with small correlation length and high variance are reported.