scipy-yli/tests/test_bayes_factors.py

#   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'
	
	expected_summary = '''BF10 = 1.23, Anecdotal evidence in favour of H1
H1: Fibrinogen ≠ 0
H0: Fibrinogen = 0'''
	assert bf.summary() == expected_summary
	assert bf._repr_html_() == 'BF<sub>10</sub> = 1.23, Anecdotal evidence in favour of H<sub>1</sub><br>H<sub>1</sub>: Fibrinogen ≠ 0<br>H<sub>0</sub>: 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'
	
	expected_summary = '''BF01 = 2.42, Anecdotal evidence in favour of H0
H0: GammaGlobulin = 0
H1: GammaGlobulin ≠ 0'''
	assert bf.summary() == expected_summary
	assert bf._repr_html_() == 'BF<sub>01</sub> = 2.42, Anecdotal evidence in favour of H<sub>0</sub><br>H<sub>0</sub>: GammaGlobulin = 0<br>H<sub>1</sub>: GammaGlobulin ≠ 0'