15

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 - kalythos$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.

5
  • 5
    Prediction is always in the context of some statistical model. One needs a distributional and structural assumptions before the variable can be "predicted". In the case of functions as lm and glm the independent variables are assumed to be fixed (i.e. deterministic) so the prediction of those is meaningless. If you want to draw inference on X then you will have to use some kind of hierarchical approach to make X stochastic. Most likely you will end up in a Bayesian framework which will give you the posterior for your X, which in turn you can use for predictions.
    – VitoshKa
    Commented Nov 16, 2010 at 9:26
  • 2
    You better specify what exactly you want. With 1 x, that's doable. With 2 x's, you have an infinite amount of possible answers. So I really wonder why exactly you need the inverse prediction. Is it for calibration purposes or so? - edit: see also VitoshKa's comment.
    – Joris Meys
    Commented Nov 16, 2010 at 9:27
  • you could build an inverse model, something like invM1 <- lm(x ~ y, data) and then use predict on your new predictor y. Now, before you jump in and do so, I recommend taking into account what @vitoshKa commented above.
    – PavoDive
    Commented Jun 6, 2016 at 16:43
  • You could also use the approx function to do this kind of calibration / inverse prediction, stackoverflow.com/questions/23957486/… Commented Mar 6, 2018 at 9:28
  • 1
    @PavoDive That only works for simple linear regression (one x variable).
    – NoName
    Commented Sep 12, 2021 at 18:56

6 Answers 6

10

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.

1
  • what is y and what is m?
    – NoName
    Commented Sep 12, 2021 at 18:57
6

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.

0
3

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

2

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.

0

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.

1
  • This does not seem scalable at all. Imagine doing this with 100 x's. OOF!
    – NoName
    Commented Sep 12, 2021 at 19:02
0

Given you are using the binomial model, the model is:

f(x)=\frac{1}{1+e^{-(\alpha + \beta x)}

Therefore, the inverse model is:

f^{-1}(y)=\frac{log(\frac{y}{1-y}) - \alpha}{\beta}

You can find the α and β values in the summary report of your glm under the Estimate column, or with model$coefficients. α is the Intercept coefficient and β the coefficient of your variable (x in your case).

Hence, the code you need is the following, with y the probability for which you would like the value, and model your regression model.

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

The code for the model equation is:

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.