R语言 “预测”函数的逆

c8ib6hqw  于 2023-07-31  发布在  其他
关注(0)|答案(6)|浏览(125)

使用predict(),可以获得给定模型的自变量(x)的某个值的因变量(y)的预测值。有没有一个函数可以预测给定yx
举例来说:

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

hiz5n14c

hiz5n14c1#

看到前面的答案就删除了。在你的例子中,给定n=50,模型是二项分布的,你可以用下面的公式计算x,给定y:

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

字符串
但是在这样做的时候,你最好咨询统计学家,让他告诉你如何计算逆预测区间。请你考虑一下Vitoshka的意见

vhmi4jdf

vhmi4jdf2#

偶然发现了这个旧线程,但我想我会添加一些其他信息。包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响应)是随机的。此外,如果数据没有分组,那么解释变量只有2个值:但是,即使我们假设x是固定的,对于给定的p,计算剂量x的置信区间仍然是有意义的,与@VitoshKa所说的相反。就像我们可以根据ED 50重新参数化模型一样,我们可以为ED 60或任何其他分位数这样做。参数是固定的,但我们仍然为它们计算CI。

jm2pwxwz

jm2pwxwz3#

chemcal包具有inverse.predict()函数,该函数适用于y ~ xy ~ x - 1形式的拟合

bvjxkvbb

bvjxkvbb4#

你只需要重新排列回归方程,但正如上面的评论所说,这可能很棘手,不一定有意义的解释。
但是,对于您提供的案例,您可以用途:

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

字符串
注意,我首先除以x系数,以确保name属性正确。

cxfofazt

cxfofazt5#

为了快速查看(没有间隔和考虑其他问题),您可以使用TeachingDemos包中的TkPredict函数。它并不直接这样做,但允许您动态更改x值并查看预测的y值,因此移动x直到找到所需的Y(对于给定的附加x值)是相当简单的,这也将显示多个x可能存在的问题,这些x可能适用于同一个y。

6vl6ewon

6vl6ewon6#

假设您使用的是二项式模型,则模型为:


的数据
因此,逆模型为:

您可以在Estimate列下的glm的摘要报告中找到α和β值,或使用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))
}

相关问题