How is Goodness-of-fit different from a loss function? Specifically, can we use formulas like $R^2$ to train a model? If not, what stops us from Goodness-of-fit formulas as loss functions?

Depending on how you define $R^2$, it can be seen as a monotonic function of the square loss. A common definition that I like is: $$
R^2=1-\left(\dfrac{
\overset{n}{\underset{i=1}{\sum}}\left(
y_i-\hat y_i
\right)^2
}{
\overset{n}{\underset{i=1}{\sum}}\left(
y_i-\bar y
\right)^2
}\right)
$$
Square loss is in the numerator of the fraction in parentheses.
From calculus, monotonic transformations do not change the critical points (where maxima and minima occur), and using a minus sign like that $R^2$ equation does only interchanges the maxima and minima. Consequently, minimizing square loss is equivalent to maximizing $R^2$ (or minimizing $-R^2$).
