Is R-squared truly an invalid metric for non-linear models? I have read that R-squared is invalid for non-linear models, because the relationship that SSR + SSE = SSTotal no longer holds.  Can somebody explain why this is true?
SSR and SSE are just the squared norms of the regression  and residual vectors, whose $i^{th}$ components are $(\hat{Y_i}-\bar{Y})$ and $(Y_i-\hat{Y_i})$, respectively.  As long as these vectors are orthogonal to one another, then shouldn't the above relationship always hold, no matter what kind of function that is used to map predictor values to fitted ones? 
Furthermore, shouldn't the regression and residual vectors associated with any least-squares model be orthogonal, by definition of least-squares?
The residual vector is the difference between vector $(Y_i-\bar{Y_i})$ and the regression vector.  If the regression vector is such that the residual/difference vector is not orthogonal to it, then the regression vector can be multiplied by a constant so that it is now orthogonal to the residual/difference vector.  This should also reduce the norm of the residual/difference vector.
If I have explained this poorly, then please tell me and I will try to clarify. 
 A: Why should SSE + SSR be equal to SST? It just happened to be the case for the linear model. There are many ways to show that it should hold for $y=X\beta+\varepsilon$ under Gauss-Markov conditions. However, it doesn't need to hold in general case. The burden is to prove that it holds, not that it doesn't
A: The sums-of-squares in linear regression are special cases of the more general deviance values in the generalised linear model.  In the more general model there is a response distribution with mean linked to a linear function of the explanatory variables (with an intercept term).  The three deviance statistics in a GLM are defined as:
$$\begin{matrix}
\text{Null Deviance} \quad \quad \text{ } \text{ } & & \text{ } D_{TOT} = 2(\hat{\ell}_{S} - \hat{\ell}_0), \\[6pt]
\text{Explained Deviance} & & D_{REG} = 2(\hat{\ell}_{p} - \hat{\ell}_0), \\[6pt]
\text{Residual Deviance}^\dagger \text{ } & & \text{ } D_{RES} = 2(\hat{\ell}_{S} - \hat{\ell}_{p}). \\[6pt]
\end{matrix}$$
In these expressions the value $\hat{\ell}_S$ is the maximised log-likelihood under a saturated model (one parameter per data point), $\hat{\ell}_0$ is the maximised log-likelihood under a null model (intercept only), and $\hat{\ell}_{p}$ is the maximised log-likelihood under the model (intercept term and $p$ coefficients).
These deviance statistics play a role analogous to scaled versions of the sums-of-squares in linear regression.  It is easy to see that they satisfy the decomposition $D_{TOT} = D_{REG} + D_{RES}$, which is analogous to the decomposition of the sums-of-squares in linear regression.  In fact, in the case where you have a normal response distribution with a linear link function you get a linear regression model, and the deviance statistics reduce to the following:
$$\begin{equation} \begin{aligned}
D_{TOT} = \frac{1}{\sigma^2} \sum_{i=1}^n (y_i - \bar{y})^2 = \frac{1}{\sigma^2} \cdot SS_{TOT}, \\[6pt]
D_{REG} = \frac{1}{\sigma^2} \sum_{i=1}^n (\hat{y}_i - \bar{y})^2 = \frac{1}{\sigma^2} \cdot SS_{REG}, \\[6pt]
D_{RES} = \frac{1}{\sigma^2} \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \frac{1}{\sigma^2} \cdot SS_{RES}. \\[6pt]
\end{aligned} \end{equation}$$
Now, the coefficient of variation in a linear regression model is a goodness-of-fit statistic that measures the proportion of the total variation in the response that is attributable to the explanatory variables.  A natural extension in the case of a GLM is to form the statistic:
$$R_{GLM}^2 = 1-\frac{D_{RES}}{D_{TOT}} = \frac{D_{REG}}{D_{TOT}}.$$
It is easily seen that this statistic reduces to the coefficient of variation in the special case of linear regression, since the scaling values cancel out.  In the broader context of a GLM the statistic has a natural interpretation that is analogous to its interpretation in linear regression: it gives the proportion of the null deviance that is explained by the explanatory variables in the model.
Now that we have seen how the sums-of-squares in linear regression extend to the deviances in a GLM, we can see that the regular coefficient of variation is inappropriate in the non-linear model, since it is specific to the case of a linear model with a normally distributed error term.  Nevertheless, we can see that although the standard coefficient of variation is inappropriate, it is possible to form an appropriate analogy using the deviance values, with an analogous interpretation.

$^\dagger$ The residual deviance is sometimes just called the deviance.
A: While R-squared may still be a flawed measurement in non-linear models for other reasons, I believe I have sufficiently shown that the SSR + SSE = SSTotal relationship still holds in a least-squares model for certain non-linear functions, especially those that allow for a constant term, such as polynomial models.  I believe that this conclusion is compatible with what has been posted in this discussion, including what I read from the ncbi link provided, although I was unable to access the full report.
If one has a series of fitted values $\hat y_i$ with respect to a series of observations $y_i$, where $\hat y_i$ $ = A + f(X) = $ $\bar Y$ $ + (A-\bar Y)$ $+ f(X) $, with $A$ being a constant term and $f(X)$ a function of predictor variables, in which the vector of $(\hat{Y_i} - \bar{Y})$ is not orthogonal to $(Y_i - \hat{Y_i})$, one can create a new set of fitted values $Z_i$ such that $Z_i = c*(\hat{Y_i} - \bar{Y}) + \bar{Y}$, where c = $\sum{(\hat{Y_i}-\bar{Y})*(Y_i-\hat{Y_i})} / \sum{(\hat{Y_i} - \bar{Y})^2}$.  With new fitted values $Z_i$, the vector $(Z_i - \bar{Y})$ will be orthogonal to the error vector and this new error vector $(Y_i - Z_i)$ will have a smaller sum of squares than the original $(Y_i-\hat{Y_i})$.  The $Z_i$ were simply obtained by multiplying the original estimated model by a constant $"c"$ and adding a multiple of the observations' mean, which is compatible with the model having a constant term.  Therefore a least-squares model should always have orthogonal regression and error vectors in these circumstances, which means that $SSE + SSR = SSTotal$.
I have created polynomial models on a handful of datasets at work and this relationship has held with all of them.  I am just saying.
A: $R^2$ is of limited use in nonlinear regression. We make it available in GraphPad Prism, but suggest it be used in only one way:
Look at $R^2$ when you run a series of experiments, and you want to make sure that today's experiment is consistent with other runs of the experiment. For example, if you always get $R^2$ between 0.90 and 0.95 but today you got $R^2$=0.75, then you should be suspicious and look carefully to see if something went wrong with the methods or reagents used in that particular experiment. And if a new employee brings you results showing $R^2$ of 0.99 using that same system, you should look carefully at how many "outliers" were removed, and whether some data were made up. 
More. 
