You asked for a proof from "first principles". So let's do it. I'll highlight the most common sources of errors and I'll show an alternative proof later that doesn't require any knowledge of tensor calculus or Einstein notation.
The hard way
First, the coordinates convention:
$$(r,\theta,\phi) \rightarrow (x,y,z) = (r\sin\theta\cos\phi,\;r\sin\theta\sin\phi,\;r\cos\theta)$$
The same way we can express $(x,y,z)$ as $x\,\mathbf{\hat e}_x + y\,\mathbf{\hat e}_y + z\,\mathbf{\hat e}_z$, we can also express $(r,\theta,\phi)$ as $r'\,\mathbf{\hat e}_r + \theta'\,\mathbf{\hat e}_\theta + \phi'\,\mathbf{\hat e}_\phi$, but now the coefficients are not the same: $(r',\theta',\phi') \neq (r,\theta,\phi)$, in general. This is because spherical coordinates are curvilinear, so the basis vectors are not the same at all points. For small variations, however, they are very similar. More precisely, relative to a point $\vec{\mathbf p}_0 = (x,y,z)$, a neighbor point $\vec{\mathbf p}_1 = (x+\Delta x,\;y+\Delta y,\;z+\Delta z)$ can be described by $\Delta \vec{\mathbf p} = (\Delta x,\Delta y,\Delta z)$ and, in spherical coordinates, if this variation is "infinitesimal", then $d\vec{\mathbf p} = (dr, d\theta, d\phi) = dr\,\mathbf{\hat e}_r + d\theta\,\mathbf{\hat e}_\theta + d\phi\,\mathbf{\hat e}_\phi$. This is basically the motivation for defining the (unnormalized) basis as:
$$
\vec{\mathbf e}_r = \frac{\partial\vec{\mathbf p}}{\partial r}, \quad
\vec{\mathbf e}_\theta = \frac{\partial\vec{\mathbf p}}{\partial \theta}, \quad
\vec{\mathbf e}_\phi = \frac{\partial\vec{\mathbf p}}{\partial \phi}
$$
But this is not normalized yet. Coincidentally, $||\partial\vec{\mathbf p}/\partial r||$ turns out to be $1$, but $||\partial\vec{\mathbf p}/\partial \theta|| = r$, as we'll see. So the actual basis should be defined as:
$$
\hat{\mathbf e}_r = \frac{\vec{\mathbf e}_r}{||\vec{\mathbf e}_r||}, \quad
\hat{\mathbf e}_\theta = \frac{\vec{\mathbf e}_\theta}{||\vec{\mathbf e}_\theta||}, \quad
\hat{\mathbf e}_\phi = \frac{\vec{\mathbf e}_\phi}{||\vec{\mathbf e}_\phi||}
$$
Explicitly:
$$\begin{align}
\vec{\mathbf e}_r &=\;( & \sin\theta\cos\phi&, & \sin\theta\sin\phi&, & \cos\theta) \\
\vec{\mathbf e}_\theta &=\;( & r\cos\theta\cos\phi&, & r\cos\theta\sin\phi&, & -r\sin\theta) \\
\vec{\mathbf e}_\phi &=\;( & -r\sin\theta\sin\phi&, & r\sin\theta\cos\phi&, & 0)
\end{align}$$
$$\begin{align}
||\vec{\mathbf e}_r||^2 &= \sin^2\theta(\cos^2\phi + \sin^2\phi) + \cos^2\theta & &= 1 \\
||\vec{\mathbf e}_\theta||^2 &= r^2\cos^2\theta(\cos^2\phi + \sin^2\phi) + r^2\sin\theta & &= r^2 \\
||\vec{\mathbf e}_\phi||^2 &= r^2\sin^2\theta(\sin^2\phi + \cos^2\phi) & &= r^2\sin^2\theta
\end{align}$$
$$\begin{align}
\hat{\mathbf e}_r &= \vec{\mathbf e}_r & &= & &( & \sin\theta\cos\phi&, & \sin\theta\sin\phi&, & \cos\theta) \\
\hat{\mathbf e}_\theta &= \vec{\mathbf e}_\theta/r & &= & &( & \cos\theta\cos\phi&, & \cos\theta\sin\phi&, & -\sin\theta) \\
\hat{\mathbf e}_\phi &= \vec{\mathbf e}_\phi/(r\sin\theta) & &= & &( & -\sin\phi&, & \cos\phi&, & 0)
\end{align}$$
You may verify that this also forms an orthogonal basis (hence orthonormal). For example:
$$\begin{align}
\hat{\mathbf e}_r \cdot \hat{\mathbf e}_\theta &= \sin\theta\cos\theta\cos^2\phi + \sin\theta\cos\theta\sin^2\phi - \sin\theta\cos\theta \\
&= 0
\end{align}$$
That does not need to happen in general.
To go from one set of coordinates to the other using the basis vectors, we solve:
$$
\begin{bmatrix}\hat{\mathbf e}_r \\ \hat{\mathbf e}_\theta \\ \hat{\mathbf e}_\phi\end{bmatrix} =
\begin{bmatrix}
\sin\theta\cos\phi & \sin\theta\sin\phi & \cos\theta \\
\cos\theta\cos\phi & \cos\theta\sin\phi & -\sin\theta \\
-\sin\phi & \cos\phi & 0
\end{bmatrix}
\begin{bmatrix}\hat{\mathbf e}_x \\ \hat{\mathbf e}_y \\ \hat{\mathbf e}_z\end{bmatrix}
$$
for $\hat{\mathbf e}_x$, $\hat{\mathbf e}_y$, and $\hat{\mathbf e}_z$ in terms of $\hat{\mathbf e}_r$, $\hat{\mathbf e}_\theta$, and $\hat{\mathbf e}_\phi$. Then any vector $\vec{\mathbf p} = x\,\mathbf{\hat e}_x + y\,\mathbf{\hat e}_y + z\,\mathbf{\hat e}_z$ can be written in the form $r'\,\mathbf{\hat e}_r + \theta'\,\mathbf{\hat e}_\theta + \phi'\,\mathbf{\hat e}_\phi$ by simple substitution. Since this particular basis is orthonormal, there's an alternative way: simply use the dot product. For example, to get $r'$:
$$\begin{align}
\vec{\mathbf p}\cdot\mathbf{\hat e}_r &= r'\,\mathbf{\hat e}_r\cdot\mathbf{\hat e}_r + \theta'\,\mathbf{\hat e}_\theta\cdot\mathbf{\hat e}_r + \phi'\,\mathbf{\hat e}_\phi\cdot\mathbf{\hat e}_r \\
&= r'
\end{align}$$
Now to the gradient. Using matrix notation, we can write the gradient as a row vector and the formula for the chain rule becomes:
$$\begin{align}
\vec\nabla f &= \begin{bmatrix}\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z}\end{bmatrix} \\
&= \begin{bmatrix}\frac{\partial f}{\partial r} & \frac{\partial f}{\partial \theta} & \frac{\partial f}{\partial \phi}\end{bmatrix}
\begin{bmatrix}
\frac{\partial r}{\partial x} & \frac{\partial r}{\partial y} & \frac{\partial
r}{\partial z} \\
\frac{\partial \theta}{\partial x} & \frac{\partial \theta}{\partial y} & \frac{\partial \theta}{\partial z} \\
\frac{\partial \phi}{\partial x} & \frac{\partial \phi}{\partial y} & \frac{\partial \phi}{\partial z}
\end{bmatrix}
\end{align}$$
Call the matrix on the right $J$ (it's the Jacobian matrix). Note that this also works the other way around too:
$$
\begin{bmatrix}\frac{\partial f}{\partial r} & \frac{\partial f}{\partial \theta} & \frac{\partial f}{\partial \phi}\end{bmatrix} \\
= \begin{bmatrix}\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z}\end{bmatrix}
\begin{bmatrix}
\frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} & \frac{\partial
z}{\partial r} \\
\frac{\partial x}{\partial \theta} & \frac{\partial y}{\partial \theta} & \frac{\partial z}{\partial \theta} \\
\frac{\partial x}{\partial \phi} & \frac{\partial y}{\partial \phi} & \frac{\partial z}{\partial \phi}
\end{bmatrix}
$$
And call this other matrix $J'$. We can invert the first equation to get $\vec\nabla f \,J^{-1} = \vec\nabla f \,J'$ $\Rightarrow$ $\vec\nabla f \,\left(J^{-1}-J'\right) = 0$. Since this works for an arbitrary $f$, we have $J^{-1} - J' = 0$ $\Rightarrow$ $J' = J^{-1}$. An important consequence is that, in general:
$$
\frac{\partial a}{\partial b} \neq \left(\frac{\partial b}{\partial a}\right)^{-1}
$$
It appears that the OP made this mistake in a comment, confusing $\partial r/\partial x$ with $(\partial x/\partial r)^{-1} = 1/(\sin\theta\cos\phi)$, as would be the case if we were using regular (instead of partial) derivatives.
Now we have two ways of calculating the matrix $J$. Directly or by calculating $J'$ first and then inverting it. Let's do it directly. We're gonna need the expressions for $r$, $\theta$, and $\phi$ in terms of $x$, $y$, and $z$ (for other coordinate systems this might be very difficult to obtain):
$$\begin{align}
r &= \sqrt{x^2+y^2+z^2} \\
\theta &= \arctan\left(\frac{\sqrt{x^2+y^2}}{z}\right) \\
\phi &= \arctan\left(\frac{y}{x}\right)
\end{align}$$
The partial derivatives are:
$$\begin{align}
\frac{\partial r}{\partial x} &= \frac{x}{\sqrt{x^2+y^2+z^2}} & &= \sin\theta\cos\phi \\
\frac{\partial r}{\partial y} &= \frac{y}{\sqrt{x^2+y^2+z^2}} & &= \sin\theta\sin\phi \\
\frac{\partial r}{\partial z} &= \frac{z}{\sqrt{x^2+y^2+z^2}} & &= \cos\theta
\end{align}$$
$$\begin{align}
\frac{\partial \theta}{\partial x} &= \frac{zx}{\sqrt{x^2+y^2}\left(x^2+y^2+z^2\right)} & &= \frac{\cos\theta\cos\phi}{r} \\
\frac{\partial \theta}{\partial y} &= \frac{zy}{\sqrt{x^2+y^2}\left(x^2+y^2+z^2\right)} & &= \frac{\cos\theta\sin\phi}{r} \\
\frac{\partial \theta}{\partial z} &= -\frac{\sqrt{x^2+y^2}}{x^2+y^2+z^2} & &= -\frac{\sin\phi}{r}
\end{align}$$
$$\begin{align}
\frac{\partial \phi}{\partial x} &= -\frac{y}{x^2+y^2} & &= -\frac{\sin\phi}{r\sin\theta} \\
\frac{\partial \phi}{\partial y} &= \frac{x}{x^2+y^2} & &= \frac{\cos\phi}{r\sin\theta} \\
\frac{\partial \phi}{\partial z} &= 0 & &= 0
\end{align}$$
Our Jacobian is then:
$$
J = \begin{bmatrix}
\sin\theta\cos\phi & \sin\theta\sin\phi & \cos\theta \\
\frac{\cos\theta\cos\phi}{r} & \frac{\cos\theta\sin\phi}{r} & -\frac{\sin\phi}{r} \\
-\frac{\sin\phi}{r\sin\theta} & \frac{\cos\phi}{r\sin\theta} & 0
\end{bmatrix}
$$
Alternatively, we could have calculated the inverse Jacobian (which is straightforward) and then inverted it (which is a nightmare). We can use Wolfram Alpha to confirm that it gives the same result:

Finally, we use the dot product to find the coefficients $r'$, $\theta'$, and $\phi'$:
$$
r' = \vec\nabla f\cdot\mathbf{\hat e}_r = \begin{bmatrix}\frac{\partial f}{\partial r} & \frac{\partial f}{\partial \theta} & \frac{\partial f}{\partial \phi}\end{bmatrix} J \begin{bmatrix}\sin\theta\cos\phi \\ \sin\theta\sin\phi \\ \cos\theta\end{bmatrix} = \begin{bmatrix}\frac{\partial f}{\partial r} & \frac{\partial f}{\partial \theta} & \frac{\partial f}{\partial \phi}\end{bmatrix} \begin{bmatrix}1 \\ 0 \\ 0\end{bmatrix} = \frac{\partial f}{\partial r}
$$
$$
\theta' = \vec\nabla f\cdot\mathbf{\hat e}_\theta = \begin{bmatrix}\frac{\partial f}{\partial r} & \frac{\partial f}{\partial \theta} & \frac{\partial f}{\partial \phi}\end{bmatrix} J \begin{bmatrix}\cos\theta\cos\phi \\ \cos\theta\sin\phi \\ -\sin\theta\end{bmatrix} = \begin{bmatrix}\frac{\partial f}{\partial r} & \frac{\partial f}{\partial \theta} & \frac{\partial f}{\partial \phi}\end{bmatrix} \begin{bmatrix}0 \\ 1/r \\ 0\end{bmatrix} = \frac{1}{r}\frac{\partial f}{\partial \theta}
$$
$$
\phi' = \vec\nabla f\cdot\mathbf{\hat e}_\phi = \begin{bmatrix}\frac{\partial f}{\partial r} & \frac{\partial f}{\partial \theta} & \frac{\partial f}{\partial \phi}\end{bmatrix} J \begin{bmatrix}-\sin\phi \\ \cos\phi \\ 0\end{bmatrix} = \begin{bmatrix}\frac{\partial f}{\partial r} & \frac{\partial f}{\partial \theta} & \frac{\partial f}{\partial \phi}\end{bmatrix} \begin{bmatrix}0 \\ 0 \\ 1/(r\sin\theta)\end{bmatrix} = \frac{1}{r\sin\theta}\frac{\partial f}{\partial \phi}
$$
Therefore:
$$
\vec\nabla f = \frac{\partial f}{\partial r}\mathbf{\hat e}_r + \frac{1}{r}\frac{\partial f}{\partial \theta}\mathbf{\hat e}_\theta + \frac{1}{r\sin\theta}\frac{\partial f}{\partial \phi}\mathbf{\hat e}_\phi
$$
A much better way
We're gonna need a new notation to avoid having to use different letters for $x$, $y$, and $z$, for example. Instead let's use indices from $1$ to $3$. For Cartesian coordinates we'll use the letter $x$, and for spherical coordinates we'll use the letter $r$. The following should be self explanatory:
$$
\vec{\mathbf p} = \sum_i x_i\mathbf{\hat x}^i = \sum_k r_k\mathbf{\hat r}^k
$$
From the definition of the basis vectors:
$$
\mathbf{\vec r}^k = \frac{\partial \vec{\mathbf p}}{\partial r_k}, \quad \mathbf{\hat r}^k = \frac{\mathbf{\vec r}^k}{||\mathbf{\vec r}^k||} = \frac{1}{h_k}\frac{\partial \vec{\mathbf p}}{\partial r_k}
$$
Where $h_k \triangleq ||\mathbf{\vec r}^k||$. Expanding in the $x$ basis:
$$
\mathbf{\hat r}^k = \sum_j\frac{1}{h_k}\frac{\partial x_j}{\partial r_k}\mathbf{\hat x}^j
$$
Now the gradient is just:
$$
\vec\nabla f = \mathbf{\vec F} = \sum_i F_i\mathbf{\hat x}^i = \sum_i \frac{\partial f}{\partial x_i}\mathbf{\hat x}^i
$$
To get the $k$'th component in spherical coordinates ($F'_k$), use the dot product:
$$\begin{align}
F'_k &= \mathbf{\vec F} \cdot \mathbf{\hat r}^k \\
&= \left(\sum_i \frac{\partial f}{\partial x_i}\mathbf{\hat x}^i\right) \cdot \left(\sum_j\frac{1}{h_k}\frac{\partial x_j}{\partial r_k}\mathbf{\hat x}^j\right) \\
&= \frac{1}{h_k}\sum_i \frac{\partial f}{\partial x_i}\frac{\partial x_i}{\partial r_k} \\
&= \frac{1}{h_k}\frac{\partial f}{\partial r_k}
\end{align}$$
and we're done.