1
$\begingroup$

I would like to get the inverse function of my GAM. Here is my modeling code:

y <- c(0.0000943615,0.0074918919,0.0157332851,0.0783308615,
       0.1546375803,0.5558444681,0.8583806898,0.9617216854,
       0.9848004112,0.9964662546)
x <- log(c(0.05, 0.1, 0.15, 0.2, 0.4, 0.8, 1.6, 3.2, 4.5, 6.4))
fit.gam <- gam(y~s(x,k=-1, bs="cr"),family=betar(link="logit"))

I used the natural cubic regression spline with bs="cr".

I would like to know the x and its confidence interval that leads to y=0.1.

The most straightforward solution is to find the corresponding knots (x_j and x_j+1) and the formula when x_j < x < x_j+1, then x=f(beta|y,x_j,x_j+1). However, when I check the formula in the book written by Dr. Wood, it is not easy to write the formula down.

On Pg.145, the book indicates the formula for cubic regression splines (after formula (4.3)) to include a matrix of F, which is consist of other two matrices B and D. I can write the other parts down, but I am not sure about B and D matrix. For B, I think it is a 2x2 matrix, but for D, the illustration looks like it is a 1x3 matrix?

Could you give me an insight on whether it is a good way to find the inverse function? If yes, how should I construct the D matrix?

Thanks!

Edits based on comments:

To BalaGiezh:

Below are the screenshots I got from the book Generalized Additive Models: an introduction with R by Simon N. Wood. enter image description here

enter image description here

B and D matrices are used to construct the F matrix. My goal is to write the f(x) formula (the one under equation (4.3)) out so that I can manually do the inverse function. I hope this helps.

To Roland:

Yes. I would like to do inverse prediction with the coefficients fitted by GAM.

To Firebug:

Here, I assume my regression model is monotone. It comes with some outliers, so I use GAM to help me penalize the parameter estimates.

$\endgroup$
3
  • $\begingroup$ Could you maybe elaborate a little bit more on what you mean by matrix B, D etc.? $\endgroup$
    – BalaGizeh
    Oct 13 at 5:47
  • $\begingroup$ "I would like to know the x and its confidence interval that leads to y=0.1." Do you want to do inverse prediction? $\endgroup$
    – Roland
    Oct 14 at 6:47
  • $\begingroup$ It's not always possible to do that. The reason is that the function can be surjective, i.e. the same y can have multiple x. $\endgroup$
    – Firebug
    Oct 14 at 14:55
0
$\begingroup$

This is a very crude, practical solution but it may be good enough while we wait for a proper answer... Predict the values of y at a grid of points on x and then find the value on the grid closest to your target.

y <- c(0.0000943615,0.0074918919,0.0157332851,0.0783308615,
       0.1546375803,0.5558444681,0.8583806898,0.9617216854,
       0.9848004112,0.9964662546)
x <- log(c(0.05, 0.1, 0.15, 0.2, 0.4, 0.8, 1.6, 3.2, 4.5, 6.4))

fit.gam <- mgcv::gam(y ~ s(x,k=-1, bs="cr"))

Predict at a grid of x values:

x_grid <- seq(min(x), max(x), by= 0.01)
pd <- as.data.frame(predict(fit.gam, data.frame(x= x_grid), se.fit= TRUE))
pd$x_grid <- x_grid

Which is the value of x giving y = 0.1?

target_y <- 0.1
idx <- which.min(abs(pd$fit - target_y))
pd[idx,]
       fit se.fit x_grid
178 0.0997 0.0176  -1.23

answer: at x = -1.23 you get y = 0.0997 (close enough?) with s.e. of the prediction 0.0176.

It could be scripted to get closer and closer to your target y.

$\endgroup$
2
  • $\begingroup$ Thanks, dariober! I am sure this is effective, but my concern here is I want to know the se of x, so I have to get the inverse function x=f(beta|y,x_j,x_j+1). $\endgroup$
    – XFkuma
    Oct 14 at 16:52
  • $\begingroup$ @XFkuma I don't think the following is correct but I'll throw it in... At y = 0.1 you have s.e. = 0.0176. So an approximate 95% CI could be the value of x giving y = 0.1 - 1.96 * 0.0176 = 0.07, that is x_low_ci = -1.57; and similarly for the upper bound you want the value of x giving y = 0.1 + 1.96 * 0.0176 = 0.13, that is x_up_ci = -1.05. That is: the interval (-1.57 -1.05) is compatible with y = 0.1. (Again, while we wait for the proper answer...) $\endgroup$
    – dariober
    Oct 14 at 16:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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