# Установка и подключение необходимых библиотек R
if(!require(simstudy)) install.packages("simstudy") #для создания набора данных
if(!require(epiR)) install.packages("epiR") #для задачи №1
if(!require(DescTools)) install.packages("DescTools") #для задачи №1
if(!require(pROC)) install.packages("pROC") #для задач №2,3
if(!require(ggplot2)) install.packages("ggplot2") # для задачи №4
if(!require(scales)) install.packages("scales") # для задачи №4
library(simstudy) #для создания набора данных
library(epiR) #для задачи №1
library(DescTools) # для задачи №1
library(pROC) #для задач №2,3
library(ggplot2) #для задачи №4
library(scales) #для задачи №4

# R коды для создания набора исходных данных.
# Для создания набора данных используется пакет "simstudy". Для воспроизводимости результатов этой задачи задается параметр генерации (333). Генерируются 1000 записей с идентификатором "id" каждой строки, затем в набор добавляется две коррелирующих (r=0.9) друг с другом непрерывных переменных "TFLW" и "TRES" с заданными значениями их средних (120,-0.5) и стандартных отклонений (30,0.6). На основе переменной "TFLW" с помощью готовой логистической модели генерируется дихотомическая переменая "VR", а на основе ее значений - дихотомическая переменая "OUT". В конце кода изменяется порядок переменных для нового набора "df" и для контроля выводятся первые строки сформированного набора данных.
set.seed(333)  
num <- 1000
df<-genData(num, dtDefs = NULL, id = "id")
dtAdd<-addCorData(df, idname="id", mu=c(120,-0.5), sigma=c(30,0.6), corMatrix = NULL, rho=0.9,  corstr = "ind", cnames = c("TFLW","TRES"))
dtAdd$TFLW<-round(dtAdd$TFLW,0)
dtAdd$TRES<-round(dtAdd$TRES,2)
nFrml <- genFormula(c(-6.36295, 0.05477), c("TFLW"))
defVR <- defDataAdd(varname="VR", formula = nFrml, dist = "binary", link = "logit")
dtAdd <- addColumns(defVR, dtAdd)
defOUT <- defCondition(condition = "VR == 1",
                       formula = ".1;.9", dist = "categorical")
defOUT <- defCondition(defOUT, condition = "VR == 0",
                       formula = ".9;.1", dist = "categorical")
dtAdd <- addCondition(defOUT, dtAdd, "OUT")
dtAdd$OUT <- ifelse(dtAdd$OUT==1,0,1)
df<-as.data.frame(dtAdd)
View(df)
df<-df[,c("id","OUT","VR","TFLW","TRES")]
View(df)
head(df)

#Задача 1.  Сравнение исследуемого диагностического теста (переменная  - "VR") с тестом «золотого стандарта» (переменная  - "OUT") по показателям клинической информативности с помощью рассчитанных Se, Sp, NPV, PPV, Acc, LR+ и LR-.
# Для создания таблицы 2x2 из исходного набора данных «df»,  требуется подсчет количества строк в наборе данных «df» по четырем шаблонам: если «11», то это – TP, «00» – TN, «10» – FP, и «01» – FN. Например, если у одного пациента исследуемый тест "VR" дал положительный результат (1), а тест «золотого стандарта» "OUT" отрицательный (0), то эта запись соответствует шаблону «10» и относится к  ложноположительному результату FP. Подсчет частот каждого из вариантов результата (с помощью функции табуляции «table()») позволяет получить матрицу результатов (или решений) теста, транспонировать ее в привычный вид, когда значения TP находятся в верхней левой ячейке (с помощью функции «Rev()» из пакета «DescTools»), дать названия строкам и столбцам матрицы с помощью функций colnames() и rownames(), применить функцию «epi.tests()» из пакета «epiR» для расчета показателей клинической информативности ДМТ (их точечных оценок и 95% доверительных интервалов) и вывести результат с помощью «summary()» в виде машинограммы матрицы исходов. 
dat<-table(df$VR, df$OUT)
dat<-Rev(dat, margin=c(1,2))
colnames(dat) <- c("Dis+","Dis-")
rownames(dat) <- c("Test+","Test-")
rval <- epi.tests(dat, conf.level = 0.95)
print(rval); summary(rval)

#Задача 2. Построение ROC кривых для двух диагностических медицинских технологий – использование показателей воздушного потока "TFLW"  и показателей давления "TRES"; расчет и представление AUC (на одном рисунке), нахождение диагностических пороговых точек для TFLW и TRES (на каждом свою). 
# Для построения ROC диаграммы была использована функция plot.roc() из пакета "pROC". На квадратной диаграмме (ее форма определяется параметром pty = "s") отображена пороговая точка со значениями специфичности и чувствительности теста при ее использовании. Вертикальные планки погрешности отображают 95% доверительный интервал для чувствительности теста в этой точке, горизонтальные – 95% доверительный интервал для специфичности теста в этой точке. R код для показателя сопротивления "TRES" аналогичен. 
#ROC кривая для df$TFLW - переменная "TFLW"
par(pty = "s")
rocTFLW<-plot.roc(df$OUT, df$TFLW,
percent=TRUE,
ci=TRUE, of="thresholds",
thresholds="best",
print.thres="best",
auc.polygon=TRUE, grid=TRUE, auc.polygon.col=rgb(.3,0,.8,0.2))
coords(rocTFLW, "best", ret=c("threshold", "sens", "spec", "ppv", "npv"))
#ROC кривая для df$TRES - переменная "TRES"
par(pty = "s")
rocTRES<-plot.roc(df$OUT, as.numeric(df$TRES),
percent=TRUE,
ci=TRUE, of="thresholds",
thresholds="best",
print.thres="best",
auc.polygon=TRUE, grid=TRUE, auc.polygon.col=rgb(.3,0,.8,0.2))
coords(rocTRES, "best", ret=c("threshold", "sens", "spec", "ppv", "npv"))

#Задача 3. Проверка статистической гипотезы о равенстве AUC двух указанных тестов.
# Для построения кривых ROC на одной диаграмме была также использована функция plot.roc() с дополнительными аргументами, отвечающими за цвет линии, цвет надписи со значением AUC и координатами этой надписи для каждой из кривой. Для проверки статистической гипотезы о равенстве AUC двух указанных тестов использовалась функция roc.test() из того же пакета "pROC" с выводимым в виде текстовой надписи значением уровня значимости p – "p-value". 
roctTFLW <- plot.roc(df$OUT, df$TFLW,
                     print.auc = TRUE, col = "green", print.auc.col = "green", 
                     print.auc.y = 0.7, print.auc.x = 0.7)
roctTRES <- plot.roc(df$OUT, as.numeric(df$TRES), 
                     print.auc = TRUE, col = "blue", print.auc.col = "blue", add = TRUE, 
                     print.auc.y = 0.5, print.auc.x = 0.5)
testroc <- roc.test(roctTFLW, roctTRES)
text(1, 1, labels=paste("p-value =", format.pval(testroc$p.value)), adj=c(0, .5))
legend("bottomright", legend=c("Total Flow", "Total Resist"), col=c("green", "blue"), lwd=2)

#Задача 4. Таблица №4 с R кодом для создания номограммы Фагана.
# Готовая функция "nomogrammer.r" загружена из внешнего источника. В качестве аргумента для построения линий, связывающих предтестовую и последестовую вероятность болезни, указаны значения чувствительности (Sens = 0.92) и специфичности (Spec = 0.89) теста,полученные при решения задачи №1. Произвольно выбранное для данного примера значение предтестовой вероятности болезни составляет 5% (Prevalence=0.05).
source("https://raw.githubusercontent.com/achekroud/nomogrammer/master/nomogrammer.r")
nomogrammer(Prevalence = .05, Detail = TRUE, NullLine=TRUE, LabelSize=4,
            Sens = 0.92,
            Spec = 0.89)

#Ниже приведен R код, не описанный в статье для создания объединенной ROC диаграммы двух тестов с cut-off point, AUC и p-level
roctTFLW <- plot.roc(df$OUT, df$TFLW,
                     percent=TRUE,transpose = TRUE,
                     ci=TRUE, of="thresholds",
                     thresholds="best",
                     print.thres="best",
                     print.auc = TRUE, col = "green", print.auc.col = "green", 
                     print.auc.y = 40, print.auc.x = 50)
roctTRES <- plot.roc(df$OUT, as.numeric(df$TRES), 
                     percent=TRUE, transpose = TRUE,
                     ci=TRUE, of="thresholds",
                     thresholds="best",
                     print.thres="best",
                     print.auc = TRUE, col = "blue", print.auc.col = "blue", add = TRUE, 
                     print.auc.y = 35, print.auc.x = 50)
testroc <- roc.test(roctTFLW, roctTRES)
text(50, 25, labels=paste("p-value =", format.pval(testroc$p.value)), adj=c(0, .5))
legend("bottomright", legend=c("Total Flow", "Total Resist"), col=c("green", "blue"), lwd=2)