In Eulcidean 3-space with coordinates $(r, \theta, \phi)$ where $\theta$ is the polar angle and $\phi$ the azimuthal angle, we may write the covariant divergence of a vector $E = E^\mu e\_\mu$ as 
$$E^{i}\_{||i} = \frac{1}{\sqrt{g}} \frac{\partial(\sqrt{g}E^i)}{\partial x^i}$$ where $g$ is the metric determinant: $r^4\sin^2(\theta)$. Often in physics one prefers to work with unit vectors instead of the natural basis for spherical coordinates, and so we make a change of frame to $$\hat{e\_r} = e\_r,\text{ } \hat{e\_\theta} = \frac{e\_\theta}{r}, \text{ } \hat{e\_\phi} = \frac{e\_\phi}{r\sin(\theta)}$$


Which is accompanied by the change of coordinates: $$\hat{E^r} = E^r,\text{ } \hat{E^\theta} = rE^\theta, \text{ } \hat{E^\phi} = r\sin(\theta)E^\phi$$


However it is common to want to work with components in the hatted system but take derivatives with respect to the unhatted coordinates, and so one finds an expression in terms of those quantities: 
$$\vec{\nabla} \cdot \vec{E} = E^{i}\_{||i} = E^{i}\_{;i} \\= \frac{1}{\sqrt{g}} \frac{\partial(\sqrt{g}E^i)}{\partial x^i} = \frac{1}{r^{2}} \frac{\partial\left(r^{2} E^{r}\right)}{\partial r}+\frac{1}{ \sin \theta} \frac{\partial}{\partial \theta}\left(E^{\theta} \sin \theta\right)+\frac{\partial E^{\phi}}{\partial \phi} \\ =\frac{1}{r^{2}} \frac{\partial\left(r^{2} \hat{E^{r}}\right)}{\partial r}+\frac{1}{r \sin \theta} \frac{\partial}{\partial \theta}\left(\hat{E^{\theta}} \sin \theta\right)+\frac{1}{r \sin \theta} \frac{\partial \hat{E^{\phi}}}{\partial \phi}$$


However one can then analyse $$\hat{E^{i}}\_{||i} = \frac{1}{\sqrt{\hat{g}}} \frac{\partial(\sqrt{\hat{g}}\hat{E^i})}{\partial \hat{x^i}}$$ and a quick analysis of our new frame's metric shows $\hat{g} = 1$ and thus: 
$$\hat{\nabla}\cdot\hat{E} = \frac{\partial \hat{E^{r}}}{\partial \hat{r}}+\frac{\partial \hat{E^{\theta}}}{\partial \hat{\theta}}+ \frac{\partial \hat{E^{\phi}}}{\partial \hat{\phi}}$$
$$
=\frac{\partial \hat{E^{r}}}{\partial r}+\frac{1}{r}\frac{\partial \hat{E^{\theta}}}{\partial \theta}+ \frac{1}{r\sin(\theta)}\frac{\partial \hat{E^{\phi}}}{\partial \phi} \neq \vec{\nabla} \cdot \vec{E}
$$


Why are these divergences not equal? I am taking the divergence of the same vector in both coordinates and then transforming such that there derivatives are with respect to the same coordinates of the same components. What is going wrong here?



Keeping the covariant divergence intact under changes of frame
The problem is that you made a change of frame but assumed that there is a corresponding coordinate system behind the new frame. The hatted vector fields are not holonomic, and hence there does not exists a "hatted coordinate system $\hat{x}$", and the assumption that you can use the standard expression for the divergence in a coordinate system is therefore invalid. 


Another way to think about this: geometrically, the divergence of a vector field is


$$ \nabla\cdot( \sum E^i e\_i ) $$


where $E^i$ are the coefficients. By Leibniz rule you get


$$ = \sum \nabla\_{e\_i} E^i + \sum E^i \nabla\cdot(e\_i) $$


when you performed the change of variables and wrote the expression $\hat{\nabla}\cdot \hat{E}$, you only accounted for the terms in the first sum and not the terms in the second sum (which vanish for $e\_i$ being the Euclidean coordinate vector fields). 


