We present numerical scaling results for tile energy level statistics in orthogonal and symplectic tight-binding Hamiltonian random matrix ensembles defined on disordered two and three-dimensional electronic systems with and without spinorbit coupling (SOC), respectively. In the metallic phase for weak disorder the nearest level spacing distribution function P(S), the number variance <(6N)^2>, and the two-point correlation function K_2(e) are shown to be described by the Gaussian random matrix theories. In the insulating phase, for strong disorder, the correlations vanish for large scales and the ordinary Poisson statistics is asymptotically recovered, which is consistent with localization of the corresponding eigenstates. At the Anderson metal-insulator transition we obtain new universal scale-invariant distribution functions describing the critical spectral density fluctuations.