我有这样的数据集:
x <- c(0, 40, 80, 120, 160, 200)
y <- c(6.52, 5.10, 4.43, 3.99, 3.75, 3.60)
我使用lm()
计算了一个线性模型:
model <- lm(y ~ x)
我想知道如果我有新的y
值,例如ynew < - c(5.5, 4.5, 3.5)
,那么x
的预测值是多少,但如果我使用predict()
函数,它只计算新的y
值。如何预测新的x
值如果我有新的y
值?我认为你只需要使用代数将y=a+b*x
反转为x=(y-a)/b
:
cc <- coef(model)
(xnew <- (ynew-cc[1])/cc[2])
# [1] 31.43007 104.76689 178.10372
plot(x,y
abline(model)
points(xnew,ynew,col=2)
看着你这里的“数据”,我认为非线性回归可能更好...
由于这是化学中的典型问题(从校准中预测值),因此chemCal
软件包提供了inverse.predict
功能。但是,此功能仅适用于“具有模型公式y〜x或y〜x-1的单变量模型对象类lm或rlm。”
x <- c(0, 40, 80, 120, 160, 200)
y <- c(6.52, 5.10, 4.43, 3.99, 3.75, 3.60)
plot(x,y)
model <- lm(y ~ x)
abline(model)
require(chemCal)
ynew <- c(5.5, 4.5, 3.5)
xpred<-t(sapply(ynew,function(y) inverse.predict(model,y)[1:2]))
# Prediction Standard Error
#[1,] 31.43007 -38.97289
#[2,] 104.7669 -36.45131
#[3,] 178.1037 -39.69539
points(xpred[,1],ynew,col="red")
警告:如果您需要逆向预测大量值,则此函数速度较慢且不适合使用。
如果我没记错的话,负SE的出现是因为该函数期望斜率始终为正。 SE的绝对值仍然是正确的。
library(sos); findFn("{inverse prediction}")
可以找到这个函数,以及 quantchem
包中的一个类似函数(该函数似乎也可以进行非线性反演...) - Ben Bolker如果你的关系是非单调的,或者你有多个预测值,那么对于给定的y值可能会有多个x值,你需要决定如何处理。
一种慢速的方法(也许是其他提到的包中使用的方法)是使用uniroot函数:
x <- runif(100, min=-1,max=2)
y <- exp(x) + rnorm(100,0,0.2)
fit <- lm( y ~ poly(x,3), x=TRUE )
(tmp <- uniroot( function(x) predict(fit, data.frame(x=x)) - 4, c(-1, 2) )$root)
library(TeachingDemos)
plot(x,y)
Predict.Plot(fit, 'x', data=data.frame(x=x), add=TRUE, ref.val=tmp)
TeachingDemos
软件包中的TkPredict
函数来估算一个解决方案。approxfun
或splinfun
函数来快速得到近似值。tmpx <- seq(min(x), max(x), length.out=250)
tmpy <- predict(fit, data.frame(x=tmpx) )
tmpfun <- splinefun( tmpy, tmpx )
tmpfun(4)