Use bootstrap deviance for goodness of fit

notebooks/logit_mono-cubic4paper.ipynb

@@ -6160,7 +6160,7 @@
    ],
    "source": [
     "# goodness of fit measures\n",
-    "df_fit[['scale', 'dataset', 'AIC', 'BIC', 'A', 'n', 'k', '|cost_jac|', 'dof', 'cost','chi2', \"p-value(chi2)\", 'success']]"
+    "df_fit[['scale', 'dataset', 'AIC', 'BIC', 'A', 'n', 'k', '|cost_jac|', 'dof', 'cost', 'deviance', \"p-value(deviance_bootstrap)\", 'success']]"
    ]
   },
   {

tests/test_logistic.py

@@ -6,8 +6,8 @@ from irae_risk.logistic import LogisticPolyRegression
 
 @pytest.fixture
 def regression_data():
-    x = np.array([-1.0, 0.0, 1.0, 2.0])
-    y = np.array([0, 0, 1, 1])
+    x = np.linspace(-2.0, 2.0, 20)
+    y = np.array([0, 1] * 10)
     theta = np.array([-0.2, 0.8])
     return x, y, theta
 
@@ -16,7 +16,7 @@ def test_goodness_of_fit_uses_unpenalized_likelihood_by_default(regression_data)
     x, y, theta = regression_data
     model = LogisticPolyRegression(degree=1, lam=(0.0, 0.5))
 
-    result = model.goodness_of_fit(x, y, theta)
+    result = model.goodness_of_fit(x, y, theta, bootstrap_samples=20)
     expected_llf = -model.get_nllf(x, y, theta)
     k = len(theta)
     n = len(x)
@@ -30,12 +30,15 @@ def test_goodness_of_fit_can_include_regularization(regression_data):
     x, y, theta = regression_data
     model = LogisticPolyRegression(degree=1, lam=(0.0, 0.5))
 
-    unregularized = model.goodness_of_fit(x, y, theta)
+    unregularized = model.goodness_of_fit(
+        x, y, theta, bootstrap_samples=20
+    )
     regularized = model.goodness_of_fit(
         x,
         y,
         theta,
         regularization=True,
+        bootstrap_samples=20,
     )
     penalty = model.penalty(theta)
 
@@ -55,12 +58,42 @@ def test_regularization_switch_has_no_effect_without_penalty(regression_data):
     x, y, theta = regression_data
     model = LogisticPolyRegression(degree=1)
 
-    default = model.goodness_of_fit(x, y, theta)
+    default = model.goodness_of_fit(x, y, theta, bootstrap_samples=20)
     regularized = model.goodness_of_fit(
         x,
         y,
         theta,
         regularization=True,
+        bootstrap_samples=20,
     )
 
     assert regularized == pytest.approx(default)
+
+
+def test_goodness_of_fit_reports_reproducible_bootstrap_deviance(regression_data):
+    x, y, theta = regression_data
+    model = LogisticPolyRegression(degree=1)
+
+    first = model.goodness_of_fit(
+        x,
+        y,
+        theta,
+        bootstrap_samples=25,
+        bootstrap_seed=123,
+    )
+    second = model.goodness_of_fit(
+        x,
+        y,
+        theta,
+        bootstrap_samples=25,
+        bootstrap_seed=123,
+    )
+
+    assert "chi2" not in first
+    assert "p-value(chi2)" not in first
+    assert first["deviance"] == pytest.approx(2 * model.get_nllf(x, y, theta))
+    assert first["p-value(deviance_bootstrap)"] == second[
+        "p-value(deviance_bootstrap)"
+    ]
+    assert 0 < first["p-value(deviance_bootstrap)"] <= 1
+    assert first["deviance_bootstrap_samples"] == 25