The problem with the divergence of the fields you wrote is that it is ill-defined in the origin. So whatever you find is valid only if $r\neq0$ and we need to manually add the value of the divergence in the origin. How do we do it? We use the fact that the integral of the divergence in the volume is equal to the flux of the vector field on a surface which encloses that volume.
We start by computing the flux $\Phi$ of your vector fields on a spherical shell $S$ of radius $R$ i.e
$$\Phi = \int dS {1\over R^n}$$
where I used the fact that the vector field is always perpendicular to the surface so we can just integrate its value at $r=R$ on the surface $S$. Of course, in spherical coordinates, $dS=R^2\sin(\theta)d\theta d\phi$ hence
$$\Phi = R^{2-n}\int \sin(\theta)d\theta d\phi = 4\pi R^{2-n} $$
where the integral I did is just the solid angle $4\pi$.
This must be correspond to the integral of the divergence inside the volume.
As you can see, the flux on the surface is not always the same and can depend on $R$. This is because, except the $n=2$ case, the other fields decrease too fast / not fast enough for the flux to always be the same at whatever distance (the surface scales as $\sim R^2$ so you need to compensate for that or you will get an $R$-dependency)
For n=2 the flux is $4\pi$. But if we integrate the divergence, we get 0. To make things work out we need to define a Dirac's delta function in the origin that fixes thing. A bit tricky but works. See comments to your questions.
For n=1 the flux is $4\pi R$ and the divergence is $1/r^2$. If we integrate the divergence over the volume we get
$$4\pi\int^R r^2dr\; {1\over r^2} = 4\pi R$$ so all good. The flux (and hence also the integral of the divergence over the volume) are $R$ dependent but both quantities are the same. All good. We don't even have to care about the behavior in the origin, it just works..
For $n=3$ (and $n>2$ in general) things get tricky. We know, as you mentioned, that the flux (and hence the integral of the div) must be positive and indeed we get $\Phi = 4\pi /R>0$ (nut decreasing with $R$ as the field goes to 0 very quick). However if we try to do the integral in the volume we find
$$-4\pi\int ^R r^2dr {1\over r^4}=\lim_{\epsilon\to 0} 4\pi\left({1\over R} -{1\over \epsilon}\right)$$
where I used $\epsilon$ to exclude the origin. This is negative and does not converge. Again, the problem is that we don't know what is happening in the origin. We just know that there must be something so big (so infinite!) that it compensates the infinity of the above integral and also its being negative, it must be something like $\sim +4\pi /\epsilon$ with $\epsilon\to 0$ (I don't actually know, but I assume in the origin there is something like $4\pi\delta(r)/r$)
More in general, the problem is that the vector field is decreasing so fast that the only relevant contribution is in the origin and we don't know it. That will fix the "minus sign" paradox. This results in the divergence pointing towards the origin because in the origin there is something huge that compensates.
Your confusion comes from the fact that you consider the divergence a measure of the direction of the field (they all point radially outwards) where instead is a measure of what behavior the flux has as you move away from the origin (bigger values of $R$): constant if $n=0$, increasing if $n<2$, decreasing if $n>2$ and that is why the divergence is 0 [constant flux], positive [increasing flux] or negative [decreasing flux]. Roughly.
The divergence is also not a measure of how much vector field is generating but rather of how much flux the source in the origin provides.