The mean density of the star is really only defined by the formula $\bar\rho=M/V=3M/4\pi R^3$. The radius of a star is a generally a very complicated function of a star's other properties. When we determine the radius in stellar models, it's only because we've solved equations that describe the structure of the whole star, and read off the value at what we define as the surface. So no simple formula in general.
That said, one can derive the approximate functional dependence for stars of various evolutionary states through the principle of homology. i.e. assuming that stars of a certain type are just rescaled versions of each other. Glancing at my old course notes, on the upper main sequence, where stars burn hydrogen principally through the CNO cycle and have radiative envelopes dominated by electron-scattering opacity, we derived $R\propto M^{15/19}$. The same principle (but with different assumptions about the star) is used to determine the location of the Hayashi track for pre-main-sequence stars, along which $R\propto M^{-7}T^{49}$. Particular formulae can be found for different types of star but the relationships between $M$ and $R$ vary wildly.
Neither the two stars you mentioned are typical main-sequence stars. R136a1 is a Wolf-Rayet star, which is basically a star that has blasted away most of its hydrogen envelope. Mass-radius relations are usually strongly dependent on mean molecular weight, which is higher without hydrogen, so the relations break down (or, rather, would have to be derived separately). But usually higher mean molecular weight gives a more compact star. UY Scuti has probably finished burning hydrogen in its core and has moved off the main sequence. So again, it'll follow a different relation.
