I've seen the correct way of finding $r(\theta,\phi)$ for the purposes of integrating for area, but that left me wondering why we can't just use:
$$ x = a \cdot \cos{\theta} \cdot \sin{\phi} \\ y = b \cdot \sin{\theta} \cdot \sin{\phi} \\ z = c \cdot \cos{\phi} $$
and then say that:
$$ r(\theta,\phi) = \sqrt{ x^{2} + y^{2} + z^{2}} \\ r(\theta,\phi) = \sqrt{ (a \cdot \cos{\theta} \cdot \sin{\phi})^{2} + (b \cdot \sin{\theta} \cdot \sin{\phi})^{2} + (c \cdot \cos{\phi})^{2}} $$
Numerically computing A = $\int_{0}^{2\pi}\int_{0}^{\pi}r(\theta,\phi)^{2}\sin(\phi)d\phi d\theta$ for area returns an incorrect value. Obviously something is wrong, but what?
$\newcommand{xphi}{\underline\phi}\newcommand{xtheta}{\underline\theta}$ The formula in the question is a formula that generates points on the ellipsoid from arbitrary parameters $\xphi$ and $\xtheta$ as follows:
First, consider a point $P$ on the unit sphere with spherical coordinates $(r,\phi,\theta) = (1, \xphi,\xtheta).$ Then stretch (or compress) the sphere along the $x,$ $y,$ and $z$ axes by factors $a,$ $b,$ and $c$ respectively. The image of $P$ under this transformation is a point $P'$ on the ellipsoid with Cartesian coordinates $$(x,y,z) = (a \cos{\xtheta} \sin{\xphi}, b \sin{\xtheta} \sin{\xphi}, c \cos{\xphi}).$$
The result of the transformation is typically that the point $P'$ is not only at a different distance from the origin than $P$ was, but is also not on the same radial from the origin as $P$ was. That is, the distortion of the sphere into the ellipsoid tends to change both the distance and direction to each of the points on the surface. If we now take the spherical coordinates $(r,\phi,\theta)$ of the transformed point $P',$ we will (typically) find that $r \neq 1,$ $\phi \neq \xphi,$ and $\theta \neq \xtheta.$ In short, although the formulas in the question use symbols $\phi$ and $\theta$ that often represent spherical coordinates, those symbols are not spherical coordinates in those formulas.
Because the parametric formulas do not use spherical coordinates, you will get incorrect results if you pretend that the parameters are spherical coordinates and integrate over them as if they were spherical coordinates.
The reason I rewrote the formulas from the question using $\xphi$ rather than $\phi$ and $\xtheta$ rather than $\theta$ was in order to be able to clearly distinguish them from the usual symbols for spherical coordinates while keeping some indication of how my formulas over these parameters correspond to the formulas in the question. Note that the article in Wolfram Mathworld goes even further in relabeling the parameters: that article uses $u$ where you used $\theta$ and $v$ where you used $\phi.$