You're confusing the expectation values of the vector $\mathbf{r}=(x,y,z)$ and its magnitude $r=\sqrt{x^2+y^2+z^2}$. Because the hydrogen hamiltonian is parity invariant, all its eigenfunctions are chosen to have a definite parity. This means that the expectation value $\langle\mathbf{r}\rangle$ will always be zero because each component is the integral of an odd function times an even probability distribution function.
The second expectation value you mention,
$$\langle n,l,m|r|n,l,m\rangle=\int \mathrm{d}x \mathrm{d}y \mathrm{d}z \psi_{n,l,m}^\ast(x,y,z)\sqrt{x^2+y^2+z^2}\psi_{n,l,m}(x,y,z),$$
is quite different, since you're not counting the direction of the distances you add. This is the same as comparing the averages of $x$ and $|x|$ over some 1D probability distribution function $p(x)$. The average of $x$ will be zero if $p(x)$ is even, while the average of $|x|$ will only vanish is $p$ is a delta function, concentrated at the origin. Thus $\langle|x|\rangle$ can be quite large while $\langle x\rangle$ is zero because of cancellations.
This is not to say that expectation values of $\mathbf{r}$ are not interesting, but one must simply be more careful. The fact that the diagonal matrix elements vanish says that the eigenstates have no permanent dipole moment - which of course they can't as they are eigenstates of an isotropic system. What doesn't vanish, however, are the transition matrix elements,
$$\langle n',l',m'|\mathbf{r}|n,l,m\rangle,$$
which play an important role in the hydrogen atom's interaction with radiation. (Specifically, they control the leading-order interaction energy, which is the dipole coupling $-q\mathbf{r}\cdot\mathbf{E}$.) Because the components of $\mathbf{r}$ are themselves spherical harmonics of degree 1 (or linear combinations of them), the matrix element above will involve the spherical integral
$$\int\mathrm{d}\Omega \, Y_{l'm'}^\ast Y_{1,\mu} Y_{lm}$$
for $\mu=-1,0,1$.
Because the spherical harmonics have a rich orthogonality structure, only a few of these integrals will survive. Specifically, you need $\Delta l=|l-l'|=1$ and $m'=m+\mu$, which has the physical content of the dipole selection rules: for a dipole coupling to radiation, only S-P, P-D, D-F, ... transitions are allowed. Further, only a few spatial orientations are allowed: $m$ can increase or decrease by one (which happens if $|\mu|=1$, and corresponds to circularly polarized light in either direction, with the electric field rotating in the $x,y$ plane) or stay the same (when $m=0$ and the light is linearly polarized along $z$).
Best Answer
You have to compute $\int d\Omega\ Y_{lm}^*(\theta,\phi)\hat{J}_iY_{lm}(\theta,\phi)$ where $d\Omega=\sin\theta d\theta d\phi$, $J_i$ are the angular momenta operators represented in position space and $Y_{lm}$ are the wavefunctions for the state $|lm\rangle$, i.e. spherical harmonics.