I was going through the Stan documentation which can be downloaded from here. I was particularly interested in their implementation of the Gelman-Rubin diagnostic. The original paper Gelman & Rubin (1992) define the the potential scale reduction factor (PSRF) as follows:
Let $X_{i,1}, \dots , X_{i,N}$ be the $i$th Markov chain sampled, and let there be overall $M$ independent chains sampled. Let $\bar{X}_{i\cdot}$ be the mean from the $i$th chain, and $\bar{X}_{\cdot \cdot}$ be the overall mean. Define, $$W = \dfrac{1}{M} \sum_{m=1}^{M} {s^2_m}, $$ where $$s^2_m = \dfrac{1}{N-1} \sum_{t=1}^{N} (X_{m t} - \bar{X}_{m \cdot})^2\,. $$ And define $B$ $$B = \dfrac{N}{M-1} \sum_{m=1}^{M} (\bar{X}_{m \cdot} - \bar{X}_{\cdot \cdot})^2 \,.$$
Define $$\hat{V} = \left(\dfrac{N-1}{N} \right)W + \left( \dfrac{M+1}{MN} \right)B\,.$$ The PSRF is estimate with $\sqrt{\hat{R}}$ where $$ \hat{R} = \dfrac{\hat{V}}{W} \cdot \dfrac{df+3}{df+1}\,,$$ where $df = 2\hat{V}/Var(\hat{V})$.
The Stan documentation on page 349 ignores the term with $df$ and also removes the $(M+1)/M$ multiplicative term. This is their formula,
The variance estimator is $$\widehat{\text{var}}^{+}(\theta \, | \, y) = \frac{N-1}{N} W + \frac{1}{N} B\,.$$ Finally, the potential scale reduction statistic is defined by $$ \hat{R} = \sqrt{\frac{\widehat{\text{var}}^{+}(\theta \, | \, y) }{W}}\,. $$
From what I could see, they don't provide a reference for this change of formula, and neither do they discuss it. Usually $M$ is not too large, and can often be as low so as $2$, so $(M+1)/M$ should not be ignored, even if the $df$ term can be approximated with 1.
So where does this formula come from?
EDIT: I have found a partial answer to the question "where does this formula come from?", in that the Bayesian Data Analysis book by Gelman, Carlin, Stern, and Rubin (Second edition) has exactly the same formula. However, the book does not explain how/why it is justifiable to ignore those terms?