Best way to plot a heavy tailed distribution? Log-log seems more conventional to plot a probability distribution to look for evidence of a heavy tail. Why is this the case? For data with a heavier tail than an exponential distribution, wouldn't log-linear be best, i.e., $log(f(x))$ vs. $x$? On that plot, an exponential function would be linear.

The classical concept of "heavy tails" that we are usually interested in is when the tails are sufficiently heavy to give infinite variance.  This occurs under power-law tail behaviour when one or both tails decay at a rate that is no faster than cubic decay.  (Conversely, the variance will be finite if both tails decay faster than cubic decay.)  You can find a detailed explanation of the log-log properties of a distribution with power-law tails in this related question.
If a distribution has a power-law tail of the form $f(x) \rightarrow c x^{-\omega-1}$ then it can be shown that the log-tail-probability and log-deviation from a fixed point $x_0$ are related in the tail by:
$$\ln \mathbb{P}(X>x) \rightarrow \text{const} - \omega \ln(x-x_0) \quad \quad \text{as } x \rightarrow \infty.$$
Now, if we let $x_{(1)} \geqslant \cdots \geqslant x_{(n)}$ be the ordered sample values in our data, we can estimate the log-tail-probability by $\ln \hat{\mathbb{P}}(X>x_{(i)}) = \ln(2i-1)- \ln(2n)$.  For values in the tails of the distribution we should therefore expect the values to follow the relationships:
$$\begin{matrix}
\text{Right tail } & & & & \quad \quad \quad \ln(2i-1) \approx \text{const} - \omega \ln|x_{(i)}-x_0|, \\
\text{Left tail} \quad & & & & \text{ } \ln(2(n-i)-1) \approx \text{const} - \omega \ln|x_{(i)}-x_0|. \\
\end{matrix}$$
As you can see, this gives us a log-log comparison that allows us to see whether the data exhibits power-law decay in the tails, and lets us estimate the rate of decay.  If $\omega \leqslant 2$ then there is decay no faster than cubic decay, which means the distribution has infinite variance --- i.e., it is heavy tailed.

Generating the tail-plot: You can generate the tail-plot for a set of data using the tailplot function in the utilities package in R.  This function takes in a set of data and produces the tail-plot for one or both tails.  (The function can also produce the Hill-plot and/or the DeSousa-Michailidis-plot, both of which are also useful for estimating rates of tail decay.  Below we generate some data from an exponential distribution and show the tail-plot for both tails.  We can see that the left-tail decays extremely rapidly (which is unsurprising, since it is bounded) and the right-tail also decays faster than cubic decay.  This plot provides strong empirical evidence that the distribution is not heavy tailed ---i.e., it has finite variance.
#Generate data from an exponential distribution
set.seed(1)
DATA <- rexp(1000, rate = 1)

#Generate the tail-plot for the data (both tails)
library(utilities)
tailplot(DATA)


