A mathematical model of gaseous slip flow in ultra-thin film gas bearings is numerically analyzed incorporating effects of surface roughness, which is characterized by fractal geometry. The Weierstrass-Mandelbrot (W-M) function is presented to represent the multiscale self-affine roughness of the surface. A modified Reynolds equation incorporating velocity slip boundary condition is applied for the arbitrary range of Knudsen numbers in the slip and transition regimes. The effects of bearing number, Knudsen number, geometry parameters of the bearing and roughness parameters on the complex flow behaviors of the gas bearing are investigated and discussed. Numerical solutions are obtained for various bearing configurations under the coupled effects of rarefaction and roughness. The results indicate that roughness has a more significant effect on higher Knudsen number (rarefaction effect) flows with higher relative roughness. Surface with larger fractal dimensions yield more frequency variations in the surface profile, which result in an obviously larger incremental pressure loss. The Poiseuille number increases not only with increasing of rarefaction effect but also with increasing the surface roughness. It can also be observed that the current study is in good agreement with solutions obtained from the linearized Boltzmann equation.