1

In fitting GAMs to some data, my interest is in using inverse estimation to find the value of the predictor $X$ that corresponds to a given value of the response $Y$.

I have written some R code to compute $X$ from $Y$ with corresponding large-sample confidence interval.

My issue now is in using simulation of GAMs to find estimates of bias and standard error for the inverse estimates. However, when I apply bootstrapping, sometimes R throws errors that can't seem to be debugged, such as "replacement has length zero".

Is this problem due to the inherent structure of GAM (in mgcv)?

I read this previous post

Can I use bootstrapping to estimate the uncertainty in a maximum value of a GAM?

and it seems similar to what I would like to be able to do (predict the $X$ given $Y$). This approach uses lpmatrix to simulate from the posterior of the GAM covariates. Can the same be done using jagam (though this approach is in itself unstable)?

compbiostats
  • 1,557
  • What do you mean "unstable" in relation to fully Bayesian posterior sampling via jagam()? Also, the specific error sounds like a bug in your wrapper; I've seen & coded several bootstrap-based approaches for the sorts of things you mention and never had a problem (beyond the know issues of bootstrapping GAMs). – Gavin Simpson Sep 20 '17 at 15:27
  • @GavinSimpson If my memory serves me correctly, Simon Wood mentions in his book (2nd ed.) that jagam can be unpredictable (I can't remember which words he used). In the post I mention above, how can I modify your solution to output the x value corresponding to any y value on would like (not necessarily the maximum) within the last for loop (i.e., in opt[i])? – compbiostats Sep 20 '17 at 15:39
  • I haven't seen anything as regards the approach being unstable, in Simon's new book nor in the note paper in which he introduces and discusses jagam(). The only drawback mentioned is that the general gibbs sampling done by JAGS will not be as efficient as using GAM-specific fully-Bayesian implementations like BayesX. – Gavin Simpson Sep 20 '17 at 15:58
  • To find x given y, you will need some criterion by which to choose x; I have seen examples where we search for the value of x that minimises the log-likelihood of observing the stated y. Basically, if you have a way to find x for a given y, rewrite it so it uses predictions via the $Xp$ matrix and then all you need to do is replace the usage of coefficients from the fitted model with draws from a multivariate normal with mean vector given by the estimated model coefficients and covariance matrix given by the covariance matrix of the estimated parameters. – Gavin Simpson Sep 20 '17 at 16:02
  • @GavinSimpson You can see a recent question I asked on Cross Validated (though it was immediately migrated to Stack Overflow): https://stackoverflow.com/questions/46352798/comparison-of-boot-outputs-in-r-numeric-vs-integer-data-types regarding the bootstrap GAM simulation. The routine works as expected for large sample sizes (and other supplied input parameters), but throws errors for smaller sample sizes. Perhaps posterior simulation will ameliorate this problem. – compbiostats Sep 22 '17 at 04:12

1 Answers1

2

The GAM, as represented in mgcv via penalized likelihood with quadratic penalties, is an empirical Bayes model. The penalties on splines can be thought of as improper Gaussian priors on (IIRC) variance parameters.

If you use the $Xp$ matrix and draws from the posterior distribution of the parameter estimates you can generate draws from the posterior distribution.

jagam() allows exactly the same thing to be done but using the full Bayesian machinery (via JAGs) and with proper priors on all aspects of the model.

The Gibbs sampling approach employed in jagam() is not specifically designed for smooth model terms, so it is not the most efficient approach, but neither is bootstrapping GAMs. Software with samplers specially tailored to estimating GAMs would be more efficient such as BayesX and INLA, which have R interfaces.

Gavin Simpson
  • 47,626