You are probably really asking about spherical harmonics, ideally specified in a fine text, such as Modern Quantum Mechanics by Sakurai and Napolitano, for instance. You don't treat these angles as operators, as you normally work in their coordinate representation, where they are pure numerical angle variables (periodic ones).
To avoid confusion with the $\hat \bullet $ notation, let us keep using that for unit vectors and use Capitals for operators, instead, and lower case for numerical variables, such as $r,\theta,\phi$.
In the spherical coordinate representation,
\begin{align}
\mathbf L &= i \hbar \left(\frac{\hat{\boldsymbol{\theta}}}{\sin(\theta)} \frac{\partial}{\partial\phi} - \hat{\boldsymbol{\phi}} \frac{\partial}{\partial\theta}\right), \\
L^2 &= -\hbar^2 \left(\frac{1}{\sin(\theta)} \frac{\partial}{\partial\theta} \left(\sin(\theta) \frac{\partial}{\partial\theta}\right) + \frac{1}{\sin^2(\theta)}\frac{\partial^2}{\partial\phi^2}\right), \\
L_z &= -i \hbar \frac{\partial}{\partial\phi}~~~.
\end{align}
Atomic wavefunctions are of the form $\langle x,y,z|n,l,m\rangle=R_{nl}(r) Y^m_l(\theta, \phi) $, and we normally use the normalization $\langle \hat r|l,m\rangle=\langle \theta,\phi|l,m\rangle =Y^m_l(\theta,\phi)$, periodic functions. Consequently, using the standard "abuse" of notation for operator action on coordinate representations,
$$
\langle \hat r|L^2|l,m\rangle= \hbar^2 l(l+1) Y^m_l(\theta,\phi) \\
\\= -\hbar^2 \left(\frac{1}{\sin(\theta)} \frac{\partial}{\partial\theta} \left(\sin(\theta) \frac{\partial}{\partial\theta}\right) + \frac{1}{\sin^2(\theta)}\frac{\partial^2}{\partial\phi^2}\right) Y^m_l(\theta,\phi),\\
\langle \hat r|L_z|l,m\rangle=\hbar m Y^m_l(\theta,\phi)= -i \hbar \frac{\partial}{\partial\phi} Y^m_l(\theta,\phi).
$$
To address your questions then, the notional spherical operators you are talking about, $R,\Theta,\Phi$, have continuous c-number eigenvalues $r,\theta,\phi$, in principle, but they need not be used explicitly. When acting on the spherical coordinate representation wavefunctions we considered, they commute with each other and produce $\Theta Y^m_l(\theta,\phi)= \theta Y^m_l(\theta,\phi)$, $\Phi Y^m_l(\theta,\phi)= \phi Y^m_l(\theta,\phi)$, so they are out of the picture---taking away troubled formal underpinnings with them.
- They most emphatically do not commute with $L^2$ and (for $\Phi$) $L_z$, in general (except for $Y^0_0$), as evident above!
I also wouldn't spend any time on the gradients $-i\hbar(\partial_r,\partial_\theta, \partial_\phi)$ which underlie the respective canonical momentum operators, as you again never need use them, in practice.
The spherical harmonics $Y^m_l$ obey well-known orthogonality relations, where you must recall the angular measure involves a factor of $\sin\theta$. What expectation value do you have in mind? Typically, $\langle l,m| \Theta|l,m\rangle= 2\pi \int_0^\pi \!\!d\theta ~(\sin\theta) ~~\theta ~~|Y^m_l|^2$.
NB Response to comment question.
The $Y^m_l$ are eigenfunctions of $L^2$ and $\Theta $, but $\theta Y^m_l$ is an eigenfunction of $\Theta$ but not $L^2$, in general, so $[L^2,\Theta]\neq 0$. This is because θ is both an eigenvalue and a variable, so on the $Y^m_l$ basis where $L^2$ is represented by the above double gradient,
$$
[L^2,\theta] Y^m_l= L^2 (\theta Y^m_l) - \theta \hbar^2 l(l+1) Y^m_l\neq 0.
$$
To bypass unfamiliar features of this formalism, you may revert to the harmonic oscillator, with rationalized unit hamiltonian H, quadratic, $2H= P^2+X^2$, so that $[H,X]\neq 0$. Now, the Hermite functions $\psi_n(x)=\langle x|\psi_n\rangle$ are eigenfunctions of the Hamiltonian, $H\psi_n(x)= E_n \psi_n(x)$, and, of course, X. But the functions $X\psi_n(x)=x\psi_n(x)$ are not eigenfunctions of the Hamiltonian, anymore, and so
$$
\tfrac{1}{2} [-\hbar^2\partial_x^2 + x^2,x] \psi_n(x) =-\hbar^2\partial_x \psi_n(x) \neq 0~.
$$
You see then that the set of all functions of x, all eigenfunctions of X, is much larger than the set of eigenfunctions of H, and the two operators fail to commute.
