具有稳健线性模型、分位数回归和机器学习方法的逆回归过程答案

【问题标题】：Inverse regression procedures with robust linear models, quantile regression, and machine learning methods具有稳健线性模型、分位数回归和机器学习方法的逆回归过程
【发布时间】：2019-08-12 20:52:38
【问题描述】：

上下文

在剂量反应模型中，我们通常会根据反应变量回归某些剂量范围，但我们真正感兴趣的是确定引发特定反应所需的剂量。通常这是通过逆回归技术完成的（即after-fitting / reparameterisation）。编辑：澄清 - 当您需要估计杀死 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)

【问题讨论】：

这听起来更像是一个统计问题而不是编程问题。如果是这样，这将更适合Cross Validated
它确实位于边界上......我把它放在这里是因为我正在寻找一个假设它不是一个新概念的编程解决方案。如果需要，很乐意迁移它。
“有兴趣确定引发特定反应所需的剂量” 我不确定我是否理解您的要求。通常人们会将 parametric（剂量-响应）模型拟合到数据中，这将允许您直接获得例如的估计值。反应量减少 50% 所需的剂量（IC50 值）；所以我不确定你为什么要使用“逆回归技术”。剂量反应模型将允许您确定“引发”任何特定反应的浓度。
我试图澄清我上面的问题；希望我在问什么更清楚，但请随时询问是否仍有任何意义。 @MauritsEvers 我从附加到我的问题的论文中的理解是，总是需要逆回归技术 - 在检疫研究中，Fieller 的公式通常用于这些目的来计算估计的间隔。您如何“直接获得例如将反应量（IC50 值）降低 50% 所需的剂量的估计值”？
至少您需要阅读?formula，这样您才能理解为什么您的所有公式都无法达到您的预期。我还质疑您是否知道使用最后两种（也许是最后三种）方法在做什么。我认为公式在这些情况下没有任何意义。

标签： r machine-learning regression prediction

【解决方案1】：

除上述我的 cmets 之外，您的数据与典型的剂量反应测量数据并不真正相似

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

这里我假设Covariate 是剂量/浓度。

每个Covariate 的不同测量值是否与不同的实验/组相关？您是否计划对不同组的多剂量反应曲线进行拟合以比较它们？

一种可能的分析策略

这里有一些东西可能会给你一些想法。我在这里使用drc，因为它允许我将“合理的”剂量反应曲线拟合到您的数据中。合理的剂量反应模型具有dose → 0 和dose → ∞ 的水平渐近线。

在这个特定示例中，我们为您的数据拟合了一个四参数 Weibull 函数。
```
library(drc)
model <- drm(Response ~ Covariate, data = df, fct = W2.4())
```

让我们绘制原始数据和模型预测（包括置信区间）

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)

我们现在可以使用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
}

这基本上是模型的反演，所以也许这就是您在原始帖子中链接到的论文的作者所称的“逆回归技术”。

让我们计算一下LD50值（即减少50%反应所需的剂量）
```
getLDx(model, x = 0.5)
#[1] 9.465188
```
从图的检查中可以看出，该值确实对应于响应为最大响应值的 50% 的剂量。

【讨论】：

我知道这不是典型的——这是我公司唯一允许我在线发布的。协变量的每个值都是受试者接受的剂量。
@André.B 好的，从统计建模的角度来看，我想说你可以做的很少。考虑到您在特定剂量下的受试者变异性，我可以保证任何模型拟合都会很差。换句话说，您可以将任何数据拟合到数据中，结果估计会有很大的不确定性，因此很难确定任何一个模型是“最佳”模型。
我明白，但幸运的是，这不是我实际使用的数据——只是我被允许发布的唯一数据（可能是因为你提到的原因）。老实说，感谢您参与我的问题并努力理解它！我知道在基于文本的讨论中可能会丢失这种语气，但我真的很感谢您的帮助。
不客气@André.B;这是一个有趣的讨论。在我的工作（生物医学药物疗效研究）中，剂量反应分析通常按照drc 论文中描述的内容进行，因此这是我最熟悉的。无论如何，祝你在进一步的探索中好运。
@André.B 我添加了一个可能的分析策略，可以为您提供一些关于如何进行的想法。