# Загрузка и подключение библиотеки 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")