The paper introduces a fourth-order augmented matched interface and boundary (AMIB) method for solving elliptic interface problems with complex interfaces and piecewise smooth coefficients in two and three dimensions. To resolve the challenge posed by non-constant coefficients within the AMIB framework, the fast Fourier transform (FFT) Poisson solver of the existing AMIB methods is replaced by a geometric multigrid method to efficiently invert the Laplacian discretization matrix. In this work, a fourth order multigrid method will be employed in the framework of the AMIB method for elliptic interface problems with variable coefficients in two and three dimensions. Based on a Cartesian mesh, the standard fourth-order finite differences are employed to approximate the first and second derivatives involved in the Laplacian with variable coefficients. Near the interface, a fourth-order ray-casting matched interface and boundary (MIB) scheme is generalized to variable coefficient problems to enforce interface jump conditions in the corrected finite difference discretization. The augmented formulation of the AMIB allows us to decouple the interface treatments from the inversion of the Laplacian discretization matrix, so that one essentially solves an elliptic subproblem without interfaces. A fourth order geometric multigrid method is introduced to solve this subproblem with a Dirichlet boundary condition, where fourth order one-sided finite difference approximations are considered near the boundary in all grid levels. The proposed multigrid method significantly enhances the computational efficiency in solving variable coefficient problems, while achieving a fourth-order accuracy in accommodating complex interfaces and discontinuous solutions.