Would the real adjusted R-squared formula please step forward? The following question, What is the adjusted R-squared formula in lm in R and how should it be interpreted?, presents different formulas for adjusted R$^2$, which were, Quote:

*

*Wherry’s formula: $1-(1-R^2)\frac{(n-1)}{(n-v)}$

*McNemar’s formula: $1-(1-R^2)\frac{(n-1)}{(n-v-1)}$

*Lord’s formula: $1-(1-R^2)\frac{(n+v-1)}{(n-v-1)}$

*Stein's formula: $1-\big[\frac{(n-1)}{(n-k-1)}\frac{(n-2)}{(n-k-2)}\frac{(n+1)}{n}\big](1-R^2)$
End quote.
Here is a real life example for which only what the OP of that prior question refers to as Wherry's formula makes any sense, i.e., $1-(1-R^2)\frac{(n-1)}{(n-v)}$
I have obtained a fit a model with 9 parameters and 10 $(x,y)$ data elements. The R$^2$ value was 0.99971297 and immediately we see the problem. If one were to divide by $n-v-1$ we would be dividing by zero, which is not an adjustment, it is a destruction.
So using the only formula that makes sense in this case, I wish to compare with a 6 parameter model and the same 10 data elements having an R$^2$-value of 0.99944407. Thus, running the numbers we find that in the first case we have an adjusted R$^2$-value of 0.9974168, and for the second case 0.9987491, which is better, but only following adjustment.
So what is the argument in favor of the Wherry formula? In point of fact if I have 10 elements and 9 parameters, I still have one degree of freedom for any model that is not exactly determined. From this I infer that somewhere along the line of the derivation of those formulas with $n-v-1$ in the denominator that the assumption of exact determinism was made for $n-v-1$, but I see evidence that this is not the case for regression. Comments welcome, as I would really like to get to the bottom of this.
In followup, @whuber asks what the formulas are and whether or not there is an intercept. The answer to the first part is that they are hitherto undocumented statistical functions being used during non-linear regression with a proportional norm error structure, and proportional norm regression. The data are censored, and are not a candidate for maximum likelihood regression. For the second part, the intercepts are, coincidentally $(0,0)$, because these are initially CDF's of semi-infinite support starting at $t=0$. However, these are piecewise statistical functions as mandated by the physics, and the concept of an intercept does not refer to a parameter in that context. If it helps, these are not $y=m x+b+ error$ but variations of $G[F(v_1,v_2,v_3,\dots,v_5,v_6)]+error*y$, if that is what you are asking.
The above question about "intercepts" may refer to a different context, we have from the answer to Adjusted 2 regression without an intercept that
$$R_{\text{adj.}}^{2}= 1 - \frac{n}{n-p} (1 - R^2).$$
This is not the same as Wherry's formula. Moreover, I do not see how the $n-1$ in the numerator changed to be $n$, nor would I ignore in $y=m x+b$ that $b$ is a parameter. So, what gives? I am clearly missing something. Why? Consider polynomial interpolation, which states that a $P_n$, which has $n+1$ coefficients, w.l.o.g., $a_0$ through $a_n$, can exactly pass through $n+1$ 2D points in the form of a Lagrange basis function. We can call this property of exact solution to occur at polynomial efficiency. Note that at polynomial efficiency the df and number of coefficients, including the constant term are equal, and that our R-value is exactly 1. Certain other functions have exact polynomial efficiency, for example, power function products, because when a power function product's logarithm is taken, it is a polynomial in $\ln(y) = P_n[\ln (x)]$. However, many other functions do not have the ability to map to exact points at polynomial efficiency. A inability to map $m$ points with $m$ parameters we can call less than polynomial mapping efficiency (in $\mathbb{R}$). Nevertheless, when we compare fit functions for adjusted R-squared, it would behoove us to use polynomial efficiency of fitting as a basis for counting df because in the polynomial case we have a deterministic formula of $m$ coefficients for $m$ 2D points. So, where does the denominator $n-v-1$ come from and which type of function with greater than polynomial efficiency is this supposed to imply?
 A: Standard implementation in R
One potential confusion comes from the role of the intercept in the correction factor used by lm. If present, it shows up in both the numerator and the denominator of the correction factor.
Quoting from an answer in the linked question:

Extracting the most relevant line you get:
ans$adj.r.squared <- 1 - (1 - ans$r.squared) * ((n - df.int)/rdf)

which corresponds in mathematical notation to:
$$R^2_{adj} = 1 - (1 - R^2) \frac{n-1}{n-p-1}$$
assuming that there is an intercept (i.e., df.int=1), $n$ is your
sample size, and $p$ is your number of predictors. Thus, your error
degrees of freedom (i.e., rdf) equals n-p-1.

Note the "assuming that there is an intercept" phrase in the first quote above. The above description holds for a situation where the number of predictors p does not include the intercept. That can be checked in the code for the lm.fit() function. The residual degrees of freedom (called df.residual in the lm.fit code but rdf in the summary.lm code) are calculated as
df.residual = n - z$rank

where n is the number of observations and z$rank is the rank of the fitted model z, the number of linearly independent parameter estimates returned. If the model has sufficient n and there is an intercept, that rank is the number of predictors other than the intercept plus 1 for the intercept. Without an intercept, that rank it is just the number of predictors.
So if p is interpreted as the number of predictors not including the intercept, the formula in the above quote holds if the model includes an intercept. If there is no intercept, then df.int = 0 in the code quoted above and the rdf in the denominator is simply n-p. Thus without an intercept you get instead:
$$R^2_{adj} = 1 - (1 - R^2) \frac{n}{n-p}$$
So one could say that are two "real adjusted $R^2$ formulas" used in lm depending on whether the intercept is included in the model. Recognizing that the denominator in either case is df.residual = n - z$rank, where z$rank is the rank of your model, should help reduce any confusion.
Where the correction factors come from
The Wikipedia page provides an explanation for where the $n-p-1$ in the denominator and the $n-1$ in the numerator come from, when the model includes an intercept:

The principle behind the adjusted ''R''2 statistic can be seen by rewriting the ordinary ''R''2 as
$$R^2 = {1-{\textit{VAR}_\text{res} \over \textit{VAR}_\text{tot}}}$$
where $\text{VAR}_\text{res} = SS_\text{res}/n$ and $\text{VAR}_\text{tot} = SS_\text{tot}/n$ are the sample variances of the estimated residuals and the dependent variable respectively, which can be seen as biased estimates of the population variances of the errors and of the dependent variable. These estimates are replaced by statistically unbiased versions: $\text{VAR}_\text{res} = SS_\text{res}/(n-p-1)$ and $\text{VAR}_\text{tot} = SS_\text{tot}/(n-1)$.

The extension to the above formula for a no-intercept model seems direct.
Other adjustments
Although the above standard $R_{adj}^2$ effectively incorporates unbiased estimates of total and residual variance, it is not itself unbiased. In particular, that $R_{adj}^2$ estimate can take on negative values even though the underlying population multiple squared correlation coefficient (between observed and linearly-predicted values) is necessarily non-negative. That point is raised in the OP and in the answer from @Carl. But as Olkin and Pratt (1958) put it in their work on the unique uniformly minimum variance unbiased estimator of $\rho^2$ (page 211):

We cannot hope for a non-negative unbiased estimator, since there is no region in the sample space having zero probability for $\rho^2=0$ and positive probability for $\rho^2 > 0$.

There is an extensive literature on attempts to obtain improved adjusted estimates. Karch (2020) recently compared 20 different versions of adjusted $R^2$, and provided an exact implementation of the Olkin-Pratt estimator. The comparison included both the original estimators and the corresponding "positive-part shrinkage estimators" proposed by Shieh (2008), which keep non-negative values as is and set negative values to 0. The Olkin-Pratt estimator was the clear choice for an unbiased minimum mean-square error (MSE), but unbiased estimators could improve on MSE. Karch notes (page 7):

No estimator had uniformly lowest MSE across all conditions. Even stronger, no estimator was uniformly best according to the maximum or average MSE perspectives. However, across all conditions, the positive-part version of the most widely used Ezekiel (adjusted $R^2$) estimator performed best both according to the maximum MSE as well as average MSE perspective.

What Karch calls the Ezekiel estimator is the formula used in the standard R implementation described above, called McNemar's formula in the OP.
"The real adjusted $R^2$ formula"
There is no single "real" formula that will work the best in all situations, as Karch demonstrates. The choice is based on what tradeoffs one wishes to make with respect to bias and MSE, and how you wish to deal with the possibility of negative sample estimates.
Nevertheless, Karch has also shown that the standard adjusted $R^2$ formula works remarkably well in practical situations with reasonable ratios of observations to predictors and non-negative population values of $\rho^2$ as low as 0.01. For more extreme situations, with even lower population $\rho^2$ values or barely more observations than predictors, the proposal in the answer from @Carl (inverse Fisher transformation of an adjusted Fisher z-transformed correlation coefficient) might have some advantages.
