# Загрузка и подключение библиотеки BayesFactor
if(!require(BayesFactor)) install.packages("BayesFactor")
if(!require(pBrackets)) install.packages("pBrackets") #использется для вывода специального вида скобок на диаграммах статьи
if(!require(yarrr)) install.packages("yarrr") ## библиотека функций, поддерживающих электронную книгу «YaRrr! The Pirate's Guide to R» (www.thepiratesguidetor.com)
if(!require(simstudy)) install.packages("simstudy") # библиотека для генерации синтетических данных

# Т.к. библиотека BayesianFirstAid загружается из GitHub, следует использовать процесс инсталляции, описанный в разделе "Installation" ресурса https://www.rdocumentation.org/packages/BayesianFirstAid. 
library(BayesianFirstAid)
library(BayesFactor)
library(pBrackets) 
library(yarrr) 
library(simstudy) 

###############
#Вывод рисунка 1
x=seq(0,100,length=200)
muCt<-30 # среднее показателя по гипотезе группы "С"
muTr<-50 # среднее по гипотезе группы "T"
sdCt<-8 # стандартное отклонение показателя по гипотезе группы "С"
sdTr<-6 # стандартное отклонение показателя по гипотезе группы "T"
Case<-48 # среднее показателя по проведенному клиническому исследованию
yCt=dnorm(x,mean=muCt,sd=sdCt)
plot(x,yCt,type="l",lwd=2,col="red",xlim = c(0, 100), ylim = c(0, 0.08), 
     ylab = "Density", xlab = "Observed Effect Size, ES", 
     axes=FALSE)
axis(side=2, labels=TRUE)
axis(1, at = c(0,20,40,60,80,100), labels = c("0", "0.2","0.4","0.6","0.8","1"))
yTr=dnorm(x,mean=muTr,sd=sdTr)
lines(x,yTr,type="l",lwd=2,col="blue")
polygon(x,pmin(yCt,yTr), col="grey")
my.funCt <- function(x) {return(dnorm(x, mean = muCt, sd = sdCt))}
my.funTr <- function(x) {return(dnorm(x, mean = muTr, sd = sdTr))}
segments(x0=Case, y0=0, x1 = Case, y1 = my.funCt(Case),
         col = "red", lty = par("lty"), lwd = 4)
segments(x0=Case, y0=0, x1 = Case, y1 = my.funTr(Case),
         col="blue", lty = par("lty"), lwd = 2)
abline(h = my.funCt(Case), col="grey", lty=3)
abline(h = my.funTr(Case), col="grey", lty=3)
abline(v=Case, col="grey", lty=3)
fTrCs<-as.character(round(my.funTr(Case),3))
fCtCs<-as.character(round(my.funCt(Case),3))
text(70, my.funTr(Case)+0.001, bquote(bold("P")[t]~"="~(.(fTrCs))), font=3, pos = 4, col="black", cex=0.8)
text(70, my.funCt(Case)+0.001, bquote(bold("P")[c]~"="~(.(fCtCs))), font=3, pos = 4, col="black", cex=0.8)
BF<-as.character(round(my.funTr(Case)/my.funCt(Case),2))
text(65, 0.02, bquote(bold("BF")[tc]~"="~paste(frac(P[t], P[c])~ " = "~frac((.(fTrCs)), (.(fCtCs))),sep="")~" = "~bold(.(BF))), font=1, pos = 4, col="black", cex=1.1)
text(Case,-0.002,label="0.48", font=3, pos = 2, col="black", cex=0.8)

# Вывод рисунка 2 (используются значения переменных для рисунка 1)
my.funBF <- function(x) {return(my.funTr(x)/my.funCt(x))}
optBF<-optimize(my.funBF, interval=c(0, 150), maximum=TRUE)
maxoptBF<-round(optBF$maximum,2)
plot(x,my.funBF(x),type="l",lwd=3,ylab = "BF(ES)", xlab = "Observed Effect Size, ES", axes=FALSE, log="y", ylim=c(0.001,10000)) 
at.y <- outer(1:9, 10^(log10(0.0001):log10(1000)))
lab.y <- ifelse(log10(at.y) %% 1 == 0, sapply(at.y, function(i) as.expression(bquote(10^.(log10(i))))), NA)
axis(2, at = at.y, labels = lab.y, las = 1)
axis(1, at = c(0,20,40,60,80,100), labels = c("0", "0.2","0.4","0.6","0.8","1"))# ось x переназначаем метки)
abline(h = 1,lwd=1,lty=1, col="red")
abline(v = Case,lwd=1,lty=3)
abline(v = optBF$maximum,lwd=1,lty=3, col="grey")
arrows(optBF$maximum,my.funBF(optBF$maximum),optBF$maximum+5,my.funBF(optBF$maximum)+3*10^3,length = .1, lwd = 1.5,code=1)
arrows(Case,my.funBF(Case),Case-10,my.funBF(Case)+20,length = .1, lwd = 1.5,code=1)
text(optBF$maximum+15, my.funBF(optBF$maximum)+5*10^3, labels = paste("MAX BF=",round(my.funBF(optBF$maximum),0),"\n"," ES: ",round(maxoptBF/100,2), sep=""), font=1)
text(Case-20,my.funBF(Case)+20, labels = paste("BF=",round(my.funBF(Case),2),"\n"," ES: ",Case/100, sep=""), font=1)

#Задача 1. Исследование эффекта от посещения сауны у больных аллергическим ринитом. 
PTrmu<-161.9 # средняя PNIF  в группе лечения, L/min
PCtmu<-103.0 # средняя PNIF  в группе контроля, L/min
PTrSD<-46.7 # SD PNIF  в группе лечения
PCtSD<-31.3 # SD PNIF  в группе контроля
n<-13 # число испытуемых в каждой группе
set.seed(1234)
df <- data.frame(GROUP=factor(rep(c("SAUN", "CTRL"), each=n)), pnif=round(c(rnorm(n, mean=PTrmu, sd=PTrSD), rnorm(n, mean=PCtmu, sd=PCtSD))))
t.test(pnif~GROUP, data = df,var.eq=TRUE) # традиционный тест проверки 
bf = ttestBF(formula = pnif~GROUP, data = df) с
bf

#Задача 2. Исследование ассоциации частоты нарушения слуха у новорожденных и временем года
deafS <- matrix(c(53,214,1736,5349), 
				nrow = 2, 
				dimnames = list(c("Summer","Not Summer"), #row
				c("Deaf","Not Deaf"))) # создание матрицы данных из двух строк, определение названия строк и столбцов
deafS
# вычисление фактора Байеса (BF) для таблиц сопряжености функцией contingencyTableBF() пакета BayesFactor
bf = contingencyTableBF(deafS, sampleType = "indepMulti", fixedMargin = "rows")
bf
# выполнение традиционного Хи-квадрат теста и альтернативного теста на основе Байесовского подхода на сравнение пропорций новорожденных с глухотой, выявленных в летний период и в другие времена года с помощью пакета BayesianFirstAid.
prop.test(deafS) # традиционный тест Хи-квадрат (аналогичный результат получается при выполнении встроенной функции chisq.test(deafS))
bayes.prop.test(deafS) # байесовская альтернатива традиционному тесту

###############
#Рисунок 6. Генерация диаграмм с значениями апостериорной вероятности разности относительных частот глухоты у родившихся в разное время года
# а) с помощью функций библиотеки BayesianFirstAid 
fit <- bayes.prop.test(deafS)
plot(fit)
# б) с помощью функций библиотеки BayesFactor
chains = posterior(bf, iterations = 10000)
SumDeafs = chains[,"pi[1,1]"] / chains[,"pi[1,*]"]
NoSumDeafs = chains[,"pi[2,1]"] / chains[,"pi[2,*]"]
plot(mcmc(SumDeafs - NoSumDeafs), main = "Change in probability of deaf birth in summer compared to other seasons")