217 lines
11 KiB
Python
217 lines
11 KiB
Python
# scipy-yli: Helpful SciPy utilities and recipes
|
|
# Copyright © 2022 Lee Yingtong Li (RunasSudo)
|
|
#
|
|
# This program is free software: you can redistribute it and/or modify
|
|
# it under the terms of the GNU Affero General Public License as published by
|
|
# the Free Software Foundation, either version 3 of the License, or
|
|
# (at your option) any later version.
|
|
#
|
|
# This program is distributed in the hope that it will be useful,
|
|
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
# GNU Affero General Public License for more details.
|
|
#
|
|
# You should have received a copy of the GNU Affero General Public License
|
|
# along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
|
|
from pytest import approx
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
import statsmodels.api as sm
|
|
|
|
import yli
|
|
from yli.regress import CategoricalTerm
|
|
|
|
def test_regress_ols_ol11_4():
|
|
"""Compare yli.regress for Ott & Longnecker (2016) example 11.4/11.7"""
|
|
|
|
df = pd.DataFrame({
|
|
'SoilPh': [3.3, 3.4, 3.4, 3.5, 3.6, 3.6, 3.7, 3.7, 3.8, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 5.0, 5.1, 5.2],
|
|
'GrowthRet': [17.78, 21.59, 23.84, 15.13, 23.45, 20.87, 17.78, 20.09, 17.78, 12.46, 14.95, 15.87, 17.45, 14.35, 14.64, 17.25, 12.57, 7.15, 7.50, 4.34]
|
|
})
|
|
|
|
result = yli.regress(sm.OLS, df, 'GrowthRet', 'SoilPh')
|
|
|
|
assert result.dof_model == 1
|
|
assert result.dof_resid == 18
|
|
assert result.f_statistic == approx(52.01, abs=0.01)
|
|
assert result.ftest().pvalue < 0.0005
|
|
assert result.rsquared == approx(0.7429, abs=0.0001)
|
|
|
|
assert result.terms['(Intercept)'].beta.point == approx(47.48, abs=0.01)
|
|
assert result.terms['(Intercept)'].pvalue < 0.0005
|
|
assert result.terms['SoilPh'].beta.point == approx(-7.86, abs=0.01)
|
|
assert result.terms['SoilPh'].pvalue < 0.0005
|
|
|
|
assert result.terms['SoilPh'].beta.ci_lower == approx(-10.15, abs=0.01)
|
|
assert result.terms['SoilPh'].beta.ci_upper == approx(-5.57, abs=0.01)
|
|
|
|
def test_regress_ols_ol13_5():
|
|
"""Compare yli.regress for Ott & Longnecker (2016) chapter 13.5"""
|
|
|
|
df = pd.DataFrame({
|
|
'C': [460.05, 452.99, 443.22, 652.32, 642.23, 345.39, 272.37, 317.21, 457.12, 690.19, 350.63, 402.59, 412.18, 495.58, 394.36, 423.32, 712.27, 289.66, 881.24, 490.88, 567.79, 665.99, 621.45, 608.8, 473.64, 697.14, 207.51, 288.48, 284.88, 280.36, 217.38, 270.71],
|
|
'D': [68.58, 67.33, 67.33, 68, 68, 67.92, 68.17, 68.42, 68.42, 68.33, 68.58, 68.75, 68.42, 68.92, 68.92, 68.42, 69.5, 68.42, 69.17, 68.92, 68.75, 70.92, 69.67, 70.08, 70.42, 71.08, 67.25, 67.17, 67.83, 67.83, 67.25, 67.83],
|
|
'T1': [14, 10, 10, 11, 11, 13, 12, 14, 15, 12, 12, 13, 15, 17, 13, 11, 18, 15, 15, 16, 11, 22, 16, 19, 19, 20, 13, 9, 12, 12, 13, 7],
|
|
'T2': [46, 73, 85, 67, 78, 51, 50, 59, 55, 71, 64, 47, 62, 52, 65, 67, 60, 76, 67, 59, 70, 57, 59, 58, 44, 57, 63, 48, 63, 71, 72, 80],
|
|
'S': [687, 1065, 1065, 1065, 1065, 514, 822, 457, 822, 792, 560, 790, 530, 1050, 850, 778, 845, 530, 1090, 1050, 913, 828, 786, 821, 538, 1130, 745, 821, 886, 886, 745, 886],
|
|
'PR': [0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1],
|
|
'NE': [1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
|
'CT': [0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0],
|
|
'BW': [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1],
|
|
'N': [14, 1, 1, 12, 12, 3, 5, 1, 5, 2, 3, 6, 2, 7, 16, 3, 17, 2, 1, 8, 15, 20, 18, 3, 19, 21, 8, 7, 11, 11, 8, 11],
|
|
'PT': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]
|
|
})
|
|
df['LNC'] = np.log(df['C'])
|
|
|
|
result = yli.regress(sm.OLS, df, 'LNC', 'D + T1 + T2 + S + PR + NE + CT + BW + N + PT')
|
|
|
|
assert result.dof_model == 10
|
|
assert result.dof_resid == 21
|
|
assert result.f_statistic == approx(13.28, abs=0.01)
|
|
assert result.ftest().pvalue < 0.00005
|
|
assert result.rsquared == approx(0.8635, abs=0.0001)
|
|
|
|
assert result.terms['(Intercept)'].beta.point == approx(-10.63398, abs=0.00001)
|
|
assert result.terms['(Intercept)'].pvalue == approx(0.0766, abs=0.0001)
|
|
assert result.terms['D'].beta.point == approx(0.22760, abs=0.00001)
|
|
assert result.terms['D'].pvalue == approx(0.0157, abs=0.0001)
|
|
assert result.terms['T1'].beta.point == approx(0.00525, abs=0.00001)
|
|
assert result.terms['T1'].pvalue == approx(0.8161, abs=0.0001)
|
|
assert result.terms['T2'].beta.point == approx(0.00561, abs=0.00001)
|
|
assert result.terms['T2'].pvalue == approx(0.2360, abs=0.0001)
|
|
assert result.terms['S'].beta.point == approx(0.00088369, abs=0.00000001)
|
|
assert result.terms['S'].pvalue < 0.0001
|
|
assert result.terms['PR'].beta.point == approx(-0.10813, abs=0.00001)
|
|
assert result.terms['PR'].pvalue == approx(0.2094, abs=0.0001)
|
|
assert result.terms['NE'].beta.point == approx(0.25949, abs=0.00001)
|
|
assert result.terms['NE'].pvalue == approx(0.0036, abs=0.0001)
|
|
assert result.terms['CT'].beta.point == approx(0.11554, abs=0.00001)
|
|
assert result.terms['CT'].pvalue == approx(0.1150, abs=0.0001)
|
|
assert result.terms['BW'].beta.point == approx(0.03680, abs=0.00001)
|
|
assert result.terms['BW'].pvalue == approx(0.7326, abs=0.0001)
|
|
assert result.terms['N'].beta.point == approx(-0.01203, abs=0.00001)
|
|
assert result.terms['N'].pvalue == approx(0.1394, abs=0.0001)
|
|
assert result.terms['PT'].beta.point == approx(-0.22197, abs=0.00001)
|
|
assert result.terms['PT'].pvalue == approx(0.1035, abs=0.0001)
|
|
|
|
def test_regress_logit_ol12_23():
|
|
"""Compare yli.regress for Ott & Longnecker (2016) chapter 12.23"""
|
|
|
|
df = pd.DataFrame({
|
|
'Unhealthy': [False, False, False, False, False, False, False, True, False, False, False, True, False, False, False, False, False, False, True, False, True, False, False, False, False, False, True, False, False, True, False, False],
|
|
'Fibrinogen': [2.52, 2.46, 2.29, 3.15, 2.88, 2.29, 2.99, 2.38, 2.56, 3.22, 2.35, 3.53, 2.65, 2.15, 3.32, 2.23, 2.19, 2.21, 5.06, 2.68, 2.09, 2.54, 2.18, 2.68, 3.41, 3.15, 3.35, 2.60, 2.28, 3.93, 2.60, 3.34],
|
|
'GammaGlobulin': [38, 36, 36, 36, 30, 31, 36, 37, 31, 38, 29, 46, 46, 31, 35, 37, 33, 37, 37, 34, 44, 28, 31, 39, 37, 39, 32, 38, 36, 32, 41, 30]
|
|
})
|
|
|
|
result = yli.regress(sm.Logit, df, 'Unhealthy', 'Fibrinogen + GammaGlobulin')
|
|
|
|
# Some numerical differences as intercept term is very negative
|
|
lrtest_result = result.lrtest_null()
|
|
assert lrtest_result.statistic == approx(7.9138, rel=0.01)
|
|
assert lrtest_result.dof == 2
|
|
assert lrtest_result.pvalue == approx(0.0191, rel=0.02)
|
|
|
|
# Can't validate intercept as very negative
|
|
|
|
expbeta_fib = np.exp(result.terms['Fibrinogen'].beta)
|
|
assert expbeta_fib.point == approx(6.756, rel=0.01)
|
|
assert expbeta_fib.ci_lower == approx(1.007, rel=0.01)
|
|
assert expbeta_fib.ci_upper == approx(45.308, rel=0.02)
|
|
assert result.terms['Fibrinogen'].pvalue == approx(0.0491, rel=0.01)
|
|
|
|
expbeta_gam = np.exp(result.terms['GammaGlobulin'].beta)
|
|
assert expbeta_gam.point == approx(1.169, abs=0.001)
|
|
assert expbeta_gam.ci_lower == approx(0.924, abs=0.001)
|
|
assert expbeta_gam.ci_upper == approx(1.477, abs=0.001)
|
|
assert result.terms['GammaGlobulin'].pvalue == approx(0.1925, rel=0.01)
|
|
|
|
expected_summary = ''' Logistic Regression Results
|
|
======================================================
|
|
Dep. Variable: Unhealthy | No. Observations: 32
|
|
Model: Logit | Df. Model: 2
|
|
Method: MLE | Df. Residuals: 29
|
|
Date: {0:%Y-%m-%d} | Pseudo R²: 0.26
|
|
Time: {0:%H:%M:%S} | LL-Model: -11.47
|
|
Std. Errors: Non-Robust | LL-Null: -15.44
|
|
| p (LR): 0.02*
|
|
======================================================
|
|
exp(β) (95% CI) p
|
|
-----------------------------------------------
|
|
(Intercept) 0.00 (0.00 - 0.24) 0.03*
|
|
Fibrinogen 6.80 (1.01 - 45.79) 0.049*
|
|
GammaGlobulin 1.17 (0.92 - 1.48) 0.19
|
|
-----------------------------------------------'''.format(result.fitted_dt)
|
|
|
|
assert result.summary() == expected_summary
|
|
|
|
def test_regress_logit_ol10_18():
|
|
"""Compare odds ratios via yli.regress for Ott & Longnecker (2016) example 10.18"""
|
|
|
|
data = [
|
|
(False, False, 250),
|
|
(True, False, 750),
|
|
(False, True, 400),
|
|
(True, True, 1600)
|
|
]
|
|
|
|
df = pd.DataFrame({
|
|
'Response': np.repeat([d[0] for d in data], [d[2] for d in data]),
|
|
'Stress': np.repeat([d[1] for d in data], [d[2] for d in data])
|
|
})
|
|
|
|
result = yli.regress(sm.Logit, df, 'Stress', 'Response', bool_baselevels=True)
|
|
|
|
assert isinstance(result.terms['Response'], CategoricalTerm)
|
|
assert result.terms['Response'].ref_category == False
|
|
|
|
expbeta = np.exp(result.terms['Response'].categories['True'].beta)
|
|
assert expbeta.point == approx(1.333, abs=0.001)
|
|
assert expbeta.ci_lower == approx(1.113, abs=0.001)
|
|
assert expbeta.ci_upper == approx(1.596, abs=0.001)
|
|
|
|
def test_regress_penalisedlogit_kleinman():
|
|
"""Compare yli.regress with yli.PenalisedLogit for http://sas-and-r.blogspot.com/2010/11/example-815-firth-logistic-regression.html"""
|
|
|
|
df = pd.DataFrame({
|
|
'Pred': [1] * 20 + [0] * 220,
|
|
'Outcome': [1] * 40 + [0] * 200
|
|
})
|
|
|
|
result = yli.regress(yli.PenalisedLogit, df, 'Outcome', 'Pred', exp=False)
|
|
|
|
assert result.dof_model == 1
|
|
assert result.terms['(Intercept)'].beta.point == approx(-2.280389)
|
|
assert result.terms['(Intercept)'].beta.ci_lower == approx(-2.765427)
|
|
assert result.terms['(Intercept)'].beta.ci_upper == approx(-1.851695)
|
|
assert result.terms['(Intercept)'].pvalue < 0.0001
|
|
assert result.terms['Pred'].beta.point == approx(5.993961)
|
|
assert result.terms['Pred'].beta.ci_lower == approx(3.947048)
|
|
assert result.terms['Pred'].beta.ci_upper == approx(10.852893)
|
|
assert result.terms['Pred'].pvalue < 0.0001
|
|
|
|
lrtest_result = result.lrtest_null()
|
|
assert lrtest_result.statistic == approx(78.95473)
|
|
assert lrtest_result.dof == 1
|
|
assert lrtest_result.pvalue < 0.0001
|
|
|
|
expected_summary = ''' Penalised Logistic Regression Results
|
|
=========================================================
|
|
Dep. Variable: Outcome | No. Observations: 240
|
|
Model: Logit | Df. Model: 1
|
|
Method: Penalised ML | Pseudo R²: 0.37
|
|
Date: {0:%Y-%m-%d} | LL-Model: -66.43
|
|
Time: {0:%H:%M:%S} | LL-Null: -105.91
|
|
Std. Errors: Non-Robust | p (LR): <0.001*
|
|
=========================================================
|
|
β (95% CI) p
|
|
---------------------------------------------
|
|
(Intercept) -2.28 (-2.77 - -1.85) <0.001*
|
|
Pred 5.99 (3.95 - 10.85) <0.001*
|
|
---------------------------------------------'''.format(result.fitted_dt)
|
|
|
|
assert result.summary() == expected_summary
|
|
|
|
# TODO: Test for logit_then_regress
|