scipy-yli/tests/test_bayes_factors.py
2023-04-17 13:24:01 +10:00

49 lines
2.2 KiB
Python

# scipy-yli: Helpful SciPy utilities and recipes
# Copyright © 2022–2023 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 pandas as pd
import yli
def test_afbf_logit_beta_zero():
"""Compare RegressionModel.bayesfactor_beta_zero for Ott & Longnecker (2016) chapter 12.23 with R BFpack"""
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(yli.Logit, df, 'Unhealthy', 'Fibrinogen + GammaGlobulin')
# model <- glm(Unhealthy ~ Fibrinogen + GammaGlobulin, data=df, family=binomial())
# bf_fit <- BF(model, hypothesis="Fibrinogen = 0")
# summary(bf_fit)
# bf_fit <- BF(model, hypothesis="GammaGlobulin = 0")
# summary(bf_fit)
bf = result.bayesfactor_beta_zero('Fibrinogen')
assert bf.factor == approx(1.229, abs=0.001)
assert bf.num_desc == 'Fibrinogen ≠ 0'
assert bf.denom_desc == 'Fibrinogen = 0'
bf = result.bayesfactor_beta_zero('GammaGlobulin')
assert bf.factor == approx(2.417, abs=0.001)
assert bf.num_desc == 'GammaGlobulin = 0'
assert bf.denom_desc == 'GammaGlobulin ≠ 0'