Hybrid density functional calculations are used to model substitutional lithium (LiZn) and interstitial lithium (Lii) defects in ZnO. The calculated transition levels and formation energies are compared to those obtained using a semilocal density functional scheme which lead to a severe band gap underestimation. It is found that the acceptor level of LiZn is deeper than has been predicted by local density functional calculations. Further, we compare the defect formation energies obtained using the two methodologies. It is shown that is the dominant form in Zn-rich and O-rich p-type material, but is still more stable than Lii in O-rich n-type material.