We present a variational method for approximating the ground state of spin models close to (Richardson-Gaudin) integrability. This is done by variationally optimizing eigenstates of integrable Richardson-Gaudin models, where the toolbox of integrability allows for an efficient evaluation and minimization of the energy functional. The method is shown to return exact results for integrable models and improve substantially on perturbation theory for models close to integrability. For large integrability-breaking interactions, it is shown how (avoided) level crossings necessitate the use of excited states of integrable Hamiltonians in order to accurately describe the ground states of general non-integrable models.