Correct R² formulas and add Efron's R²

src/statmodels.jl

@@ -197,7 +197,13 @@ Coefficient of determination (R-squared).
 For a linear model, the R² is defined as ``ESS/TSS``, with ``ESS`` the explained sum of squares
 and ``TSS`` the total sum of squares.
 """
-r2(obj::StatisticalModel) = mss(obj) / deviance(obj)
+function r2(obj::StatisticalModel)
+    Base.depwarn("The default r² method for linear models is deprecated. " *
+                 "Packages should define their own methods.", :r2)
+
+    μ = meanresponse(obj)
+    mss(obj) / sum(x -> abs2(x - μ), response(obj))
+end
 
 """
     r2(obj::StatisticalModel, variant::Symbol)
@@ -207,27 +213,37 @@ Pseudo-coefficient of determination (pseudo R-squared).
 
 For nonlinear models, one of several pseudo R² definitions must be chosen via `variant`.
 Supported variants are:
-- `:MacFadden` (a.k.a. likelihood ratio index), defined as ``1 - \\log L/\\log L0``.
-- `:CoxSnell`, defined as ``1 - (L0/L)^{2/n}``
-- `:Nagelkerke`, defined as ``(1 - (L0/L)^{2/n})/(1 - L0^{2/n})``, with ``n`` the number
-of observations (as returned by [`nobs`](@ref)).
+- `:MacFadden` (a.k.a. likelihood ratio index), defined as ``1 - \\log (L)/\\log (L_0)``;
+- `:CoxSnell`, defined as ``1 - (L_0/L)^{2/n}``;
+- `:Nagelkerke`, defined as ``(1 - (L_0/L)^{2/n})/(1 - L_0^{2/n})``;
+- `:Efron`, defined as ``1 - \\sum_i (y_i - \\hat{y}_i)^2 / \\sum_i (y_i - \\bar{y})^2``.
 
-In the above formulas, ``L`` is the likelihood of the model, ``L0`` that of the null model
-(the model including only the intercept). These two quantities are taken to be minus half
-`deviance` of the corresponding models.
+In the above formulas, ``L`` is the likelihood of the model,
+``L_0`` is the likelihood of the null model (the model with only an intercept),
+``n`` is the number of observations, ``y_i`` are the responses, 
+``\\hat{y}_i`` are fitted values and ``\\bar{y}`` is the average response.
+
+Cox and Snell's and Efron's R² should match the classical R² for linear models.
 """
-function r2(obj::StatisticalModel, variant::Symbol)
-    ll = -deviance(obj)/2
-    ll0 = -nulldeviance(obj)/2
 
-    if variant == :McFadden
-        1 - ll/ll0
-    elseif variant == :CoxSnell
-        1 - exp(2/nobs(obj) * (ll0 - ll))
-    elseif variant == :Nagelkerke
-        (1 - exp(2/nobs(obj) * (ll0 - ll)))/(1 - exp(2/nobs(obj) * ll0))
+function r2(obj::StatisticalModel, variant::Symbol)
+    if variant == :Efron
+        y = response(obj)
+        ŷ = fitted(obj)
+        μ = meanresponse(obj)
+        1 - sum(x -> abs2(x[1] - x[2]), zip(y, ŷ)) / sum(x -> abs2(x - μ), response(obj))
     else
-        error("variant must be one of :McFadden, :CoxSnell or :Nagelkerke")
+        ll = loglikelihood(obj)
+        ll0 = nullloglikelihood(obj)
+        if variant == :McFadden
+            1 - ll/ll0
+        elseif variant == :CoxSnell
+            1 - exp(2 * (ll0 - ll) / nobs(obj))
+        elseif variant == :Nagelkerke
+            (1 - exp(2 * (ll0 - ll) / obs(obj))) / (1 - exp(2 * ll0 / nobs(obj)))
+        else
+            error("variant must be one of :McFadden, :CoxSnell, :Nagelkerke or :Efron")
+        end
     end
 end
 
@@ -252,17 +268,16 @@ adjr2(obj::StatisticalModel) = error("adjr2 is not defined for $(typeof(obj)).")
 Adjusted pseudo-coefficient of determination (adjusted pseudo R-squared).
 
 For nonlinear models, one of the several pseudo R² definitions must be chosen via `variant`.
-The only currently supported variant is `:MacFadden`, defined as ``1 - (\\log L - k)/\\log L0``.
+The only currently supported variant is `:MacFadden`, defined as ``1 - (\\log (L) - k)/\\log (L0)``.
 In this formula, ``L`` is the likelihood of the model, ``L0`` that of the null model
 (the model including only the intercept). These two quantities are taken to be minus half
 `deviance` of the corresponding models. ``k`` is the number of consumed degrees of freedom
 of the model (as returned by [`dof`](@ref)).
 """
 function adjr2(obj::StatisticalModel, variant::Symbol)
-    ll = -deviance(obj)/2
-    ll0 = -nulldeviance(obj)/2
+    ll = loglikelihood(obj)
+    ll0 = nullloglikelihood(obj)
     k = dof(obj)
-
     if variant == :McFadden
         1 - (ll - k)/ll0
     else