You need a probability model.
The idea behind a ranking system is that a single number adequately characterizes a player's ability. We might call this number their "strength" (because "rank" already means something specific in statistics). We would predict that player A will beat player B when strength(A) exceeds strength(B). But this statement is too weak because (a) it is not quantitative and (b) it does not account for the possibility of a weaker player occasionally beating a stronger player. We can overcome both problems by supposing the probability that A beats B depends only on the difference in their strengths. If this is so, then we can re-express all the strengths is necessary so that the difference in strengths equals the log odds of a win.
Specifically, this model is
$$\mathrm{logit}(\Pr(A \text{ beats } B)) = \lambda_A - \lambda_B$$
where, by definition, $\mathrm{logit}(p) = \log(p) - \log(1-p)$ is the log odds and I have written $\lambda_A$ for player A's strength, etc.
This model has as many parameters as players (but there is one less degree of freedom, because it can only identify relative strengths, so we would fix one of the parameters at an arbitrary value). It is a kind of generalized linear model (in the Binomial family, with logit link).
The parameters can be estimated by Maximum Likelihood. The same theory provides a means to erect confidence intervals around the parameter estimates and to test hypotheses (such as whether the strongest player, according to the estimates, is significantly stronger than the estimated weakest player).
Specifically, the likelihood of a set of games is the product
$$\prod_{\text{all games}}{\frac{\exp(\lambda_{\text{winner}} - \lambda_{\text{loser}})}{1 + \exp(\lambda_{\text{winner}} - \lambda_{\text{loser}})}}.$$
After fixing the value of one of the $\lambda$, the estimates of the others are the values that maximize this likelihood. Thus, varying any of the estimates reduces the likelihood from its maximum. If it is reduced too much, it is not consistent with the data. In this fashion we can find confidence intervals for all the parameters: they are the limits in which varying the estimates does not overly decrease the log likelihood. General hypotheses can similarly be tested: a hypothesis constrains the strengths (such as by supposing they are all equal), this constraint limits how large the likelihood can get, and if this restricted maximum falls too far short of the actual maximum, the hypothesis is rejected.
In this particular problem there are 18 games and 7 free parameters. In general that is too many parameters: there is so much flexibility that the parameters can be quite freely varied without changing the maximum likelihood much. Thus, applying the ML machinery is likely to prove the obvious, which is that there likely are not enough data to have confidence in the strength estimates.
If you try graphing the equations on https://www.desmos.com/calculator they will make much more sense. For the NFL equation with different (ELOW-ELOL) values, say 50,100,200 for example you get get a graph probably 95% correlated to ln(x+1) which is the main part of the equation. Adding abs(x) makes sure the point differential is always expressed as a positive number so you can ignore the graph where x is negative.
The final bit of the equation, (2.2/((ELOW-ELOL)*.001+2.2)) which we'll refer to as auto_m is the most complicated and confusing term, but it has very little impact on the overall equation. This is the bit that accounts for autocorrelation. There are two cases
1)ELOW - ELOL will be positive when the favorite wins and auto_m will be a number slightly < 1
2)ELOW - ELOL will be negative when the underdog wins and auto_m will be a number slightly > 1.
The net affect is to decrease the margin of victory multiplier when the favorite wins and increase the margin of victory multiplier when the underdog wins.
For the NBA equation ((x+3)^.8)/7.5 it looks like they are using a more linear equation (x^1 would be exactly linear) so winning by 20 points vs 18 points is almost as significant as winning by 4 points vs 2 where as in the NFL they use a much quicker diminishing reward for blowouts. The general shape of the base equations ((x+3)^.8)/7.5 and ln(x+1) are very similar for the values of x>0. For a value of x=1 the multiplier is 0.40 for the NBA and 0.69 for the NFL. Using an exactly linear equation would have too big of a multiplier(as large as 3.0) for large blowouts. The multiplier should be kept between 0 and 2.0 as a rule of thumb.
A similar equation using the NFL type equation for the NBA would be something like ln(x+4)1/((ELOW-ELOL).006+3) but as I said they wanted a more linear equation which is why they use x^n instead of ln(x) though the general shape of the curve is the same r shape where x>0.
For games with less scoring like hockey or soccer where there are usually only going to be a few margin of victory outcome possibilities 0,1,2,3,4,5,6 etc you could hard code the multipliers you wanted to use for each one using a slight multiplier when the underdog achieves this margin of victory and slight discount when the favorite does.
Best Answer
Regarding "how to do it in R", the prefmod package http://cran.r-project.org/web/packages/prefmod/index.html is meant for preference analysis with paired comparisons, rankings and ratings. It fits Bradley-Terry models and pattern models with object and subject covariates. See my answer here How to fit Bradley–Terry–Luce model in R, without complicated formula? for a short intro, or this paper http://www.jstatsoft.org/v48/i10 for more info.