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$).
Hint: Rather than using spherical coordinates, which are singular where the 3D Dirac delta function has support, work instead in Cartesian coordinates $\vec{r}=(x,y,z)$ and use the defining property
$$ \iiint_{\mathbb{R}^3}\! d^3r ~f(\vec{r})~ \delta^3(\vec{r})~=~f(\vec{0}) $$
of the 3D Dirac delta function.
Best Answer
$\langle z\rangle=\int_0^\infty r^3dr[\int_0^{2\pi}d\phi\int_0^\pi \sin \theta \cos \theta d\theta]=0$ As the $\theta $ integration is zero.
But the symmetry argument is clear if the integration is written is Cartesian coordinates.In that case $$\langle z\rangle=\int_{-\infty}^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty z dz |\psi|^2 dx dy$$
As you can see the integration over $z$ is odd and therefore zero.