Well, you got it wrong. You rather need something like this, and possibly permute things to get the orderings right.
x = rand(16,872);
y = rand(4,872);
z = rand(18,872);
r = reshape(x, :, 1, 1, 872)
θ = reshape(y, 1, :, 1, 872)
ϕ = reshape(z, 1, 1, :, 872)
X = vec(r .* sin.(θ) .* cos.(ϕ))
Y = vec(r .* sin.(θ) .* sin.(ϕ))
Z = vec(r .* cos.(θ))
Why do you want to? Vectorizing this kind of code most of the time makes it harder to reason about and there should be no fundamental differences in speed. If you want to run it on GPU I could see a point. But as you can already see, it’s easy to get things wrong.