The $k$-th ($k = 3, 4, 5$) order backward differential formula (${\rm BDF}k$) is applied to develop the high order energy stable schemes for the molecular beam epitaxial model with slope selection. The numerical schemes are established by combining the convex splitting technique with the $k$-th order accurate Douglas-Dupont stabilization term in the form of $S\tau^{k−1}∆_h(\phi^n−\phi^{n−1}).$ With the help of the new constructed discrete gradient structure of the $k$-th order explicit extrapolation formula, the stabilized ${\rm BDF}k$ scheme is proved to preserve energy dissipation law at the discrete levels and unconditionally stable in the energy norm. By using the discrete orthogonal convolution kernels and the associated convolution embedding inequalities, the $L^2$ norm error estimate is established under a weak constraint of time-step size. Numerical simulations are presented to demonstrate the accuracy and efficiency of the proposed numerical schemes.