Вкратце о маленьком пакете cir.

Минутка теории

Изотоническая регрессия - это крайне специфический вид регрессии, который применяется при жестком требовании неубывания значения зависимой переменной при возрастании значения независимой переменной

Есть два основных применения данного вида регрессии в реальных задачах:

1. В фармакологических задачах, где нужно найти взаимосвязь типа "отклик - доза" (например, "концентрация препарата - доля умерших")

2. При моделировании распределений, когда в качестве зависимой переменной выступает квантиль функции распределения

Классический алгоритм расчетов дает в результате кусочно-постоянную неубывающую функцию (Pool-Adjacent-Violators Algorithm, его схема представлена ниже):

Недавно алгоритм был модифицирован. Две работы позволили немного модифицировать алгоритм:

1. Oron A. P., Flournoy N. Centered Isotonic Regression: Point and Interval Estimation for Dose-Response Studies / Statistics in Biopharmaceutical Research. 2017. № 3 (9) DOI:10.1080/19466315.2017.1286256 - был представлен новый, сглаживающий алгоритм

2. Oron A. P., Flournoy N. Bias induced by adaptive dose-finding designs // Journal of Applied Statistics. 2019. DOI: 10.1080/02664763.2019.1649375 - была представлена методика расчета доверительных интервалов при предсказании значений зависимых переменных

Новый алгоритм основан на сглаживании кусочно-постоянной регрессии:

Собственно, на всем этом и основан пакет CIR, с которым мы будем знакомиться

Методология

Основа для данной статьи - база данных из исследования (что-то связано с инфекциями в человеческой популяции), которая отдельно доступна по ссылке.

В качестве независимой переменной выступает логарифм дозы носителей инфекции, в качестве зависимой переменной - вероятность заболеть.

База данных и все сопутствующие материалы есть на GitHub.

Расчеты

Data <- read_csv(".../Dataforfigure8.csv", 
                           col_names = FALSE)
nam<-t(Data[,1])
Base<-as.data.frame(t(Data[,2:24]))
colnames(Base)<-nam
ggplot(data=Base, aes(x=Base$log_dose, y=Base$Probability_of_infection_1)) +
  geom_point() + xlab("Логарифм дозы") + ylab("Вероятность инфицирования")

Если посмотреть на точечный график, то это будет что-то типа

В пакете представлены три алгоритма регрессии - их реализуют функции cirPAVA() - это новый алгоритм, oldPAVA - старый алгоритм и iterCIR() - итеративная версия нового алгоритма. Но мы воспользуемся другой функцией для построения регрессии и доверительных интервалов

x1<-Base[,2]
names(x1)<-c()
x2<-Base[,1]
names(x2)<-c()
dat<-doseResponse(y=x1,x=x2) 
quick1<-quickIsotone(dat)  # Быстрая регрессия - получение доверительных интервалов для значений у
ggplot(data=quick1, aes(x=x, y=y)) +
  geom_point() + xlab("Логарифм дозы") + ylab("Вероятность инфицирования") + 
  geom_line(data=quick1,aes(x=x, y=lower90conf))+
  geom_line(data=quick1,aes(x=x, y=upper90conf))

При этом расчет доверительных интервалов может осуществляться тремя способами, выбрать которые можно при использовании функции isotInterval. Все вместе будет выглядеть так:

slow1<-cirPAVA(dat,full=TRUE) # Построение регрессии по алгоритму CIR
slow1$output # Предсказанные значения
slow1$input # Исходные данные
slow1$shrinkage # Сокращенные данные для построения графика зависимости
int1_0<-isotInterval(slow1,narrower=FALSE) # Расчет доверительных интервалов по готовой модели
int1_0
ggplot(data=quick1, aes(x=x, y=y)) +
  geom_point() + xlab("Логарифм дозы") + ylab("Вероятность инфицирования") +
  geom_line(data=slow1$shrinkage,aes(x=x, y=y), color = "green") +
  geom_line(data=int1_0,aes(x=quick1$x, y=ciLow))+
  geom_line(data=int1_0,aes(x=quick1$x, y=ciHigh))

Функция quickIsotone() может также осуществлять и прогнозирование. Для этого есть параметр outx

quickIsotone(dat, outx = c(2.15,7.75)) # Предсказание с доверительными интервалами для значения