'predict'函数的反函数

Question

使用predict()函数可以得到给定模型中独立变量(x)某个值时，因变量(y)的预测值。那么是否有一种函数可以根据给定的y值来预测x呢？

例如：

kalythos <- data.frame(x = c(20,35,45,55,70), 
    n = rep(50,5), y = c(6,17,26,37,44))
kalythos$Ymat <- cbind(kalythos$y, kalythos$n - kalythos$y)
model <- glm(Ymat ~ x, family = binomial, data = kalythos)

如果我们想知道模型对于 x=50 的预测值：

predict(model, data.frame(x=50), type = "response")

我想知道哪个x使得y=30，举个例子。

- danilinares

---

预测始终是在某种统计模型的背景下进行的。在变量可以被“预测”之前，需要有分布和结构假设。对于像 lm 和 glm 这样的函数，独立变量被假定为固定的（即确定性的），因此对它们的预测是没有意义的。如果您想对 X 进行推断，则必须使用某种层次方法使 X 成为随机的。最有可能的是，您将最终进入贝叶斯框架，这将为您提供 X 的后验，进而可以用于预测。 - VitoshKa

你最好明确你想要什么。有一个x，那是可行的。有两个x，你就有无限多的可能答案。所以我真的很想知道你为什么需要反向预测。是为了校准目的吗？- 编辑：请参阅VitoshKa的评论。 - Joris Meys

你可以构建一个反向模型，类似于 invM1 <- lm(x ~ y, data)，然后在你的新预测器 y 上使用 predict。现在，在你跳进去这么做之前，我建议考虑一下 @vitoshKa 在上面评论中提到的内容。 - PavoDive

你也可以使用approx函数来进行这种校准/反向预测，https://stackoverflow.com/questions/23957486/calibration-inverse-prediction-from-loess-object-in-r - Tom Wenseleers

@PavoDive 这仅适用于简单的线性回归（一个 x 变量）。 - NoName

Answer 1

之前的回答已被删除。在您的情况下，假设n=50且模型是二项式分布，则可以使用以下公式计算y给定x：

f <- function (y,m) {
  (logit(y/50) - coef(m)[["(Intercept)"]]) / coef(m)[["x"]]
}
> f(30,model)
[1] 48.59833

但是在这样做时，最好咨询一位统计学家，让他向您展示如何计算逆预测间隔。请注意VitoshKa的考虑。

- Joris Meys

---

y是什么，m又是什么？ - NoName

Answer 2

我看到了这个旧帖子，但是想补充一些其他的信息。包MASS有用于logit/probit模型的函数dose.p。标准误差(SE)通过delta方法计算。

> dose.p(model,p=.6)
             Dose       SE
p = 0.6: 48.59833 1.944772

在这里，拟合反函数模型（x~y）是没有意义的，因为正如@VitoshKa所说，我们假设x是固定的，而y（0/1响应）是随机的。此外，如果数据未分组，则解释变量只有两个值：0和1。但是，即使我们假设x是固定的，计算给定p剂量x的置信区间仍然是有意义的，与@VitoshKa所说的相反。正如我们可以重新参数化模型以ED50为单位，我们还可以将其用于ED60或任何其他分位数。参数是固定的，但我们仍然为它们计算置信区间。

- davep

Answer 3

“chemcal”包中有一个inverse.predict()函数，适用于形如y~x和y~x-1的拟合。其中，chemcal是一个链接。

- Mark Dayel

Answer 4

你只需要重新排列回归方程，但正如上面的评论所述，这可能会证明棘手，并且不一定具有有意义的解释。

然而，对于你提出的情况，你可以使用：

(1/coef(model)[2])*(model$family$linkfun(30/50)-coef(model)[1])

请注意，我首先通过 x 系数进行了除法运算，以便使名称属性正确。

- James

Answer 5

如果只是快速查看（不考虑间隔和其他问题），您可以使用TeachingDemos包中的TkPredict函数。它并不直接执行此操作，但允许您动态更改x值，并查看预测的y值，因此将x移动到找到所需的Y值（对于给定的其他x值）将相当简单，这也将显示可能存在的多个x值的问题，这些x值可用于相同的y。

- Greg Snow

---

这看起来根本不可扩展。想象一下用100个x做这件事。哎呀！ - NoName

Answer 6

鉴于您正在使用二项式模型，该模型如下：

因此，逆模型为：

你可以在你的glm的摘要报告中找到α和β值，它们位于“Estimate”列下，或者使用“model$coefficients”。α是截距系数，β是你的变量（在你的情况下是x）的系数。

因此，你需要以下代码，其中y是你想要的概率值，而model是你的回归模型。

binomial.f.inv <- function(y, model) {
  (log(y / (1 - y)) - model$coefficients[[1]]) / model$coefficients[[2]]
}

模型方程的代码如下：

binomial.f <- function(x, model) {
  1 / (1 + exp(-model$coefficients[[1]] - model$coefficients[[2]] * x))
}

- Tom