I am convinced you would not have asked any of these questions had you gone through 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 angles. 

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$. 

As you may have learned, 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=Y^m_l(\theta,\phi)$, periodic functions; so that, using the standard "abuse" of notation for [operator action on coordinate representations](https://en.wikipedia.org/wiki/Bra–ket_notation#Spinless_position–space_wave_function), 
$$
\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 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](https://en.wikipedia.org/wiki/Spherical_harmonics#Orthogonality_and_normalization), 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 functions of *x*, eigenfunctions of *X* is much larger than the set of eigenfunctions of *H*, and the two operators fail to commute.