Ensemble density functional theory (DFT) furnishes a rigorous theoretical framework for describing the non-dynamic electron correlation arising from (near) degeneracy of several electronic configurations. Ensemble DFT naturally leads to fractional occupation numbers (FONs) for several Kohn-Sham (KS) orbitals, which thereby become variational parameters of the methodology. The currently available implementation of ensemble DFT in the form of the spin-restricted ensemble-referenced KS (REKS) method was originally designed for systems with only two fractionally occupied KS orbitals, which was sufficient to accurately describe dissociation of a single chemical bond or the singlet ground state of biradicaloid species. To extend applicability of the method to systems with several dissociating bonds or to polyradical species, more fractionally occupied orbitals must be included in the ensemble description. Here we investigate a possibility of developing the extended REKS methodology with the help of the generalized valence bond (GVB) wavefunction theory. The use of GVB enables one to derive a simple and physically transparent energy expression depending explicitly on the FONs of several KS orbitals. In this way, a version of the REKS method with four electrons in four fractionally occupied orbitals is derived and its accuracy in the calculation of various types of strongly correlated molecules is investigated. We propose a possible scheme to ameliorate the partial size-inconsistency that results from perfect spin-pairing. We conjecture that perfect pairing natural orbital (NO) functionals of reduced density matrix functional theory (RDMFT) should also display partial size-inconsistency.