Процедуры обратной регрессии с надежными линейными моделями, квантильная регрессия и методы машинного обучения

Question

Context

Часто в моделях доза-ответ мы регрессируем некоторый диапазон доз от переменной ответа, но мы действительно заинтересованы в определении дозы, необходимой для получения конкретного ответа. Обычно это делается с помощью методов обратной регрессии (т. Е. после подгонки / репараметризации ). РЕДАКТИРОВАТЬ: Чтобы уточнить - это обычно делается, когда вам нужно оценить дозу, необходимую для уничтожения, скажем, 50%, или 99,99% для карантинных протоколов. Чтобы получить эти оценки, люди используют методы обратной регрессии - вышеупомянутая ссылка проходит через это более тщательно (см. Страницу 9).

Вопрос

Как я могу выполнить эти процедуры обратной регрессии, используя такие методы, как устойчивые линейные модели, квантильная регрессия или модели машинного обучения (то есть нейронные сети или машины опорных векторов)? РЕДАКТИРОВАТЬ: Чтобы уточнить, я хочу программное решение для того, как я могу оценить дозу, необходимую для получения ответа 99,99, когда модель, которую я установил, является одной из вышеупомянутых. Ниже приведены примерные модели для этих целей.

Мои данные выглядят так:

  df <- structure(list(Response = c(100, 91.1242603550296, 86.9822485207101, 
100, 0, 0, 90.5325443786982, 95.8579881656805, 88.7573964497041, 
96.4497041420118, 82.2485207100592, 99.4082840236686, 99.4082840236686, 
98.8165680473373, 91.7159763313609, 59.1715976331361, 44.9704142011834, 
0, 100, 95.2662721893491, 100, 82.8402366863905, 7.69230769230769, 
81.6568047337278, 62.7218934911243, 97.6331360946746, 73.9644970414201, 
8.87573964497041, 0, 98.8165680473373, 78.1065088757396, 98.2248520710059, 
52.6627218934911, 96.4497041420118, 52.0710059171598, 0, 62.043795620438, 
84.6715328467153, 97.8102189781022, 4.37956204379562, 89.051094890511, 
99.2700729927007, 99.2700729927007, 97.0802919708029, 81.7518248175183, 
80.2919708029197, 90.5109489051095, 99.2700729927007, 96.3503649635037, 
0, 0, 94.8905109489051, 79.5620437956204, 67.8832116788321, 73.7226277372263, 
100, 97.0802919708029, 93.4306569343066, 86.8613138686131, 33.5766423357664, 
32.1167883211679, 46.7153284671533, 98.5401459854015, 95.6204379562044, 
86.1313868613139, 14.5985401459854, 92.7007299270073, 86.1313868613139, 
0, 77.3722627737226, 89.051094890511, 80.2919708029197, 98.1818181818182, 
96.3636363636364, 30.9090909090909, 0, 60.9090909090909, 100, 
0, 83.6363636363636, 88.1818181818182, 97.2727272727273, 0, 0, 
99.0909090909091, 100, 100, 91.8181818181818, 88.1818181818182, 
46.3636363636364, 50.9090909090909, 99.0909090909091, 97.2727272727273, 
100, 0, 92.7272727272727, 60.9090909090909, 90.9090909090909, 
57.2727272727273, 76.3636363636364, 94.5454545454545, 50, 98.1818181818182, 
16.3636363636364, 87.2727272727273, 92.7272727272727, 87.2727272727273, 
88.1818181818182, 10.7438016528926, 91.7355371900827, 98.3471074380165, 
60.3305785123967, 95.8677685950413, 0, 63.6363636363636, 71.900826446281, 
0, 74.3801652892562, 76.8595041322314, 0, 61.9834710743802, 0, 
0, 0, 84.297520661157, 47.1074380165289, 69.4214876033058, 97.5206611570248, 
100, 61.1570247933884, 90.0826446280992, 78.5123966942149, 10.7438016528926, 
100, 98.3471074380165, 100, 98.3471074380165, 93.3884297520661, 
90.9090909090909, 57.8512396694215, 57.8512396694215, 92.5619834710744, 
77.6859504132231, 69.4214876033058), Covariate = c(20, 14, 14, 
20, 0, 0, 14, 14, 14, 16, 10, 20, 20, 20, 16, 10, 10, 0, 16, 
16, 16, 10, 0, 12, 10, 12, 12, 0, 0, 20, 12, 16, 10, 12, 12, 
0, 14, 14, 16, 0, 14, 20, 16, 20, 14, 12, 12, 20, 20, 0, 0, 14, 
12, 10, 10, 20, 16, 16, 14, 10, 10, 10, 20, 16, 10, 0, 12, 12, 
0, 12, 16, 14, 16, 14, 0, 0, 12, 20, 0, 12, 14, 14, 0, 0, 20, 
20, 20, 14, 14, 10, 10, 20, 16, 16, 0, 12, 10, 10, 10, 16, 16, 
12, 20, 10, 12, 12, 16, 14, 0, 16, 20, 12, 14, 10, 10, 0, 0, 
12, 12, 10, 10, 0, 0, 0, 14, 12, 12, 20, 20, 14, 14, 14, 12, 
20, 20, 20, 16, 16, 14, 10, 10, 16, 16, 16)), row.names = 433:576, class = "data.frame")

моя формула обычно выглядит примерно так:

Ответ ~ Ковариат + I (Ковариат ^ 2)

Вот пример моделей, которые я установил:

#Robust linear model
MASS::rlm(Response ~ Covariate + I(Covariate^2), data = df)

#Quantile regression
quantreg::rq(Response ~ Covariate + I(Covariate^2), data = df, tau = c(0.5, 0.95)) # In this case I want to predict the specified quantiles for the dose required to elicit a given response, although I realised this code doesn't do that...

#Machine learning algorithms were trained with caret
TRControl <- trainControl(method = "cv")

#Neural Network
caret::train(Response  ~ Covariate, data = df, method = "neuralnet", trControl = TRControl)

#Support Vector Machine
caret::train(Response  ~ Covariate, data = df, method = "polySVM", trControl = TRControl)

Answer 1

В дополнение к моим комментариям выше, ваши данные на самом деле не похожи на данные типичного измерения доза-ответ

library(ggplot2)
ggplot(df, aes(Covariate, log10(Response))) +
    geom_point()

Здесь япредположим, что Covariate это доза / концентрация.

Относятся ли разные измерения для каждого Covariate к разным экспериментам / группам?Планируете ли вы подгонять кривые зависимости от нескольких доз для разных групп, чтобы сравнить их?

---

Возможная стратегия анализа

Вот кое-что, что может дать вам некоторые идеи.Я использую drc здесь, потому что это позволяет мне подогнать «разумную» кривую доза-эффект к вашим данным.Разумная модель доза-ответ имеет горизонтальные асимптоты для dose → 0 и dose → ∞.

1. . В этом конкретном примере мы подгоняем ваши данные по четырем параметрам Вейбулла.

   library(drc)
   model <- drm(Response ~ Covariate, data = df, fct = W2.4())
   

2. Построим исходные данные и прогнозы модели (включая доверительный интервал)

   library(tidyverse)
   df.pred <- data.frame(Covariate = 1.1 * seq(min(df$Covariate), max(df$Covariate), length.out = 20)) %>%
       bind_cols(as.data.frame(predict(model, data.frame(Covariate = Covariate), interval = "confidence"))) %>%
       rename(Response = Prediction)
   
   ggplot(df, aes(Covariate, Response)) +
       geom_point() +
       geom_line(data = df.pred, aes(Covariate, Response), color = "blue") +
       geom_ribbon(data = df.pred, aes(x = Covariate, ymin = Lower, ymax = Upper), fill = "blue", alpha = 0.2)
   

   

3. Теперь мы можем использовать uniroot для определения конкретных значений LDx, которые определяются как доза, необходимая для уменьшения максимальной реакции на x / 100.

   getLDx <- function(model, x = 0.5) {
       maxResponse <- max(predict(model, data.frame(x = c(0, Inf))))
       uniroot(
           function(Covariate) predict(model, newdata = data.frame(Covariate = Covariate)) - x * maxResponse,
           interval = range(Covariate))$root
   }
   

   Это в основном инверсиямодель, поэтому, возможно, именно это авторы статей, на которые вы ссылаетесь в своем первоначальном посте, называют «методами обратной регрессии».

4. Давайте вычислим значение LD50 (то есть требуемую дозучтобы уменьшить ответ на 50%)

   getLDx(model, x = 0.5)
   #[1] 9.465188
   

   Из проверки графика видно, что это значение действительно соответствует дозе, где ответ составляет 50% от максимального значения ответа.