I'll try to clear some of the confusion, first: $$ 1= ⟨ψ|ψ⟩ = \int ⟨ψ|x⟩⟨x|ψ⟩dx=\int|\psi(x)|^2dx $$ and $$\int |x⟩⟨x|dx=I$$
In your first question you compute $|\psi(r,\theta,\phi)|^2 $, which is fine, it's the density that you're looking for. It will only equal 1 after integration over $r, \theta\ and\ \phi$.
The mistake in the second derivation is that you assumed $⟨χ_+|+⟨χ_−|=I$, while in reality $⟨χ_+|+⟨χ_−|=\sqrt2⟨1|$ (Positive x-direction spin state). You were looking for: $|χ_+⟩⟨χ_+|+|χ_-⟩⟨χ_−|=I$.
Calculating $⟨ψ|(|r,\theta,\phi⟩⟨r,\theta,\phi|⊗I)|ψ⟩$ should do the trick.
bra-ket and projectors
#bra-ket and projectors# I'llI'll just add a really short summary here, to clarify further:
⟨ψ| is a bra (dual vector), |ψ⟩ is a ket (vector in Hilbert space).
$⟨ψ|\phi⟩$ is a complex number.
$|\phi⟩⟨ψ|$ is an operator, and $|ψ⟩⟨ψ|$ is a projector on state ψ.
Operator $P$ is a projector if $P^2=P$.
If you want the probability density at position x, you project on the eigenspace of the position operator, using $|x⟩⟨x|$.
If you want the probability density at $r,θ\ and\ ϕ$ you use the projection on the eigenspace of these operators. You don't care which spin you get, so you don't project on the spin space and just use the identity.
$⟨ψ|(|r,\theta,\phi⟩⟨r,\theta,\phi|⊗I)|ψ⟩$, or $⟨ψ|(|r,\theta,\phi⟩⟨r,\theta,\phi|⊗(|χ_+⟩⟨χ_+|+|χ_-⟩⟨χ_−|))|ψ⟩=⟨ψ|(|r,\theta,\phi⟩⟨r,\theta,\phi|⊗|χ_+⟩⟨χ_+|+|r,\theta,\phi⟩⟨r,\theta,\phi|⊗|χ_-⟩⟨χ_−|)|ψ⟩$, if you prefer, will give you what you're looking for. Notice that every open bra has its closing ket, which means you will get a number, and it is symmetric, so the the number is real and positive.