The present work reports the theoretical considerations of the magnetocrystalline anisotropy of ferrous ions induced by tetravalent dopants in yttrium iron garnet. Using the spin Hamiltonian developed earlier by us and the molecular field (h) approximation we derive the cubic anisotropy constants K1 and K2 at zero temperature. We adopt the Alben’s et al. model of twelve inequivalent Fe2+ sites in silicon-substituted yttrium iron garnet. Results are given for h = 400, 300, 200 and the spin Hamiltonian parameters with the trigonal Δ = 300, 400, 500, 600, 700 and the nontrigonal crystal field parameter Γ = 200, 300 cm−1. The agreement with the experimental K1 and K2 is quite good. The discussion reveals that the properties of the far and near sites in the two-center model can now be theoretically explained. The theoretical ratios of K1(far) to K1(near) agree well with experiment. Thus our results speak in favor of the orbital singlet rather than the doublet model assumed previously for Fe2+ in silicon- or germanium-substituted yttrium iron garnets.