This dissertation presents a survey of ab initio models developed based on the merger of multireference (MR) methods and density functional theory (DFT) in order to provide an accurate and efficient description of electron correlation effects in strongly correlated systems. We first introduce a reduced density matrix (RDM)-based formulation of multiconfiguration pair-density functional theory (MC-PDFT) which addresses two of the three common problems in MR+DFT framework: double counting of the electron correlation and symmetry dilemma. MC-PDFT minimizes the double counting of electron correlation by computing the classical effects within MR part while accounting for quantum mechanical interactions via DFT part. Symmetry dilemma is also addressed through a change of variables from spin densities to on-top pair-density (OTPD). In order to resolve the third issue in MR+DFT framework, the computational cost barrier of the MR methods, we adopt the variational two-electron reduced density matrix (v2RDM)-driven complete active-space self-consistent field (CASSCF) approach. The favorable polynomial computational cost of v2RDM-CASSCF allows one to go beyond the active space size limitations of conventional configuration interaction (CI)-based MR methods. In order to reduce the delocalization error (DE) plaguing almost all density functionals, we extend MC-PDFT to its global and range-separated hybrid variants where a fraction of local exchange from OTPD functionals is replaced with its nonlocal counterpart computed by v2RDM-CASSCF reference RDMs. The efficiency and accuracy of our MC-PDFT-based models have also been demonstrated through their application to a wide variety of realistic and challenging molecular systems with dominant MR character such as the calculation of dissociation potential energy curves for di- and polyatomic molecules, reaction energy barriers of 1,3-dipolar cycloaddition reaction of ozone to ethylene and acetylene and singlet/triplet energy gaps of large members of oligocene molecular series. After providing numerical evidence for usefulness of our models for strongly correlated systems, we focus on the sources of errors and metrics for error quantification in DFT. Through introducing constrained search-Kohn-Sham density functional theory (CS-KSDFT), we have addressed two fundamental and controversial problems in KS-DFT: inaccessibility of the exact density within finite basis set and lack of a universal mathematical metric for the density error. Applying CS-KSDFT to strongly correlated systems, such as the triple-bond dissociation of N2, we have numerically showed that non-interacting KS-DFT electron densities can be far more accurate than those calculated by the conventional exchange-correlation (XC) functionals, regardless of the size of the basis set. By applying our basis set-independent metric to rank the performance of conventional density functionals for strongly correlated systems, we have numerically presented that the errors caused by the approximate forms of XC functionals (as opposed to density-driven errors) are the main source of error in KS-DFT. Lastly, we have analyzed the contentious onset of open-shell character in the singlet ground state zig-zag narrow graphene nanoribbons via a variety of metrics such as effectively unpaired electrons, natural orbital occupation numbers, singlet/triplet energy gaps and structural indicator of C-C bond length alternation.