开发者

inverse of 'predict' function

Using predict() one can obtain the predicted value of the dependent variable (y) for a certain value of the independent variable (x) for a given model. Is there any function that predicts x for a given y?

For example:

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 - kalyt开发者_JAVA百科hos$y)
model <- glm(Ymat ~ x, family = binomial, data = kalythos)

If we want to know the predicted value of the model for x=50:

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

I want to know which x makes y=30, for example.


Saw the previous answer is deleted. In your case, given n=50 and the model is binomial, you would calculate x given y using:

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

But when doing so, you better consult a statistician to show you how to calculate the inverse prediction interval. And please, take VitoshKa's considerations into account.


Came across this old thread but thought I would add some other info. Package MASS has function dose.p for logit/probit models. SE is via delta method.

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

Fitting the inverse model (x~y) would not makes sense here because, as @VitoshKa says, we assume x is fixed and y (the 0/1 response) is random. Besides, if the data weren’t grouped you’d have only 2 values of the explanatory variable: 0 and 1. But even though we assume x is fixed it still makes sense to calculate a confidence interval for the dose x for a given p, contrary to what @VitoshKa says. Just as we can reparameterize the model in terms of ED50, we can do so for ED60 or any other quantile. Parameters are fixed, but we still calculate CI's for them.


The chemcal package has an inverse.predict() function, which works for fits of the form y ~ x and y ~ x - 1


You just have to rearrange the regression equation, but as the comments above state this may prove tricky and not necessarily have a meaningful interpretation.

However, for the case you presented you can use:

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

Note I did the division by the x coefficient first to allow the name attribute to be correct.


For just a quick view (without intervals and considering additional issues) you could use the TkPredict function in the TeachingDemos package. It does not do this directly, but allows you to dynamically change the x value(s) and see what the predicted y-value is, so it would be fairly simple to move x until the desired Y is found (for given values of additional x's), this will also show possibly problems with multiple x's that would work for the same y.

0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜