Question to advanced users, math pros/programists. Can you help me understand this explanation how to use GAM outside R?

As I understand very vaguely, a matrix of parameters is being exported and then new predictions can be extrapolated. But, what is the content of lpmatrix? How it is accomplished in this code? How can it be generalized to any gam model to predict new data https://stackoverflow.com/questions/15584541/how-to-extract-fitted-splines-from-a-gam-mgcvgam/15587786#15587786

I see this question repeats often in different forms like a question how a formula/function can be extracted to use outside R, in database etc: extract specific regression equation from a GAM model? https://stackoverflow.com/questions/15584541/how-to-extract-fitted-splines-from-a-gam-mgcvgam/15587786#15587786

  • $\begingroup$ Some specific database system? Some other language You van use within the base? Then IT should be possible to generator in R code for the estimated modell in that language. Maybe a nettet for for stavk overflow $\endgroup$ Commented May 4, 2019 at 16:26
  • $\begingroup$ @kjetilbhalvorsen Say the output is grid of parameters (lpmatrix) and it can be used for prediction outside stat program, any program, can be Excel. I wonder what lpmatrix stands for and how it can be used to predict. $\endgroup$ Commented May 4, 2019 at 16:35

1 Answer 1


The way this works is that you create a matrix $\mathbf{X}_p$ such that

$$\hat{\boldsymbol{\eta}}_p = \mathbf{X}_p\hat{\boldsymbol{\beta}}$$

where $\hat{\boldsymbol{\eta}}_p$ are the predicted values on the linear predictor (i.e. on the link scale in GLM terminology). So $\mathbf{X}_p$ is basically a design matrix with the each of the smooths in the model evaluated at the covariate values wanted.

Say your model has a single smooth of covariate x, then you can create $\mathbf{X}_p$ for that model by creating say 100 or 200 values equally-spaced over the range of x, and pass that to predict() via argument newdata, with type = "lpmatrix":

Xp <- predict(mod, newdata = my_new_data, type = 'lpmatrix')

Then, given that matrix you can do a matrix multiplication with it and the vector of model coefficients (extracted by

betas <- coef(mod)

) to yield predicted values at the covariate values of x that you requested:

fit_vals <- Xp %*% betas

and back-transform if using a non-identity link via the inverse link function.

For more complex models with multiple predictors, you need to provide all the combinations of values for predictors that you want to evaluate the model at.

Then you export the matrix Xp and the vector of coefficients to whatever system you want. Or precompute the fitted values directly and export that.

How this allows you to do predictions elsewhere is that you need to map the covariate values for the new observations you want predictions for to the values that you evaluated the smooths at. Returning to the simple x example of a single smooth, you find which value of x you generated in my_new_data that is closest to the value of x for the new observation and then either do the multiplication of that row of $\mathbf{X}_p$ with the vector of coefficients to get a prediction for that new observation.

For more complex models you map the values for the new observation to the combination of values in Xp that is closest to the observed set of values for the new observation.

If the grid of values over which you evaluated the smooth(s) at and which you embedded in $\mathbf{X}_p$ is fine enough, then it won't matter too much that you are predicting at not quite the exact values for the new observation.

You can get the standard errors also if you export the variance-covriance matrix of the coeficients $\mathbf{V}_p$ or just create the standard errors directly from the square root of the diagonals of


which in R code is:

Vp <- vcov(mod)
se <- diag(Xp %*% Vp %*% t(Xp))
## or more efficiently via
se <- rowSums(Xp %*% (Vp %*% Xp))

You can use those to add an error to the prediction, by finding the index of the row of $\mathbf{X}_p$ that is closest to the new observation and then extracting the value of se indexed in the same way (say row 10 is the closest, you extract the 10th element of se, in whatever software you are doing the external predictions).

If you precalculate the fitted values for the gird of covariate values in $\mathbf{X}_p$ and their associated uncertainties, or 95% credible interval, you could simply export those values are use a lookup to identify which value in the precomputed data is closest to the new observation.

Otherwise you need to be able to do the mathematical operations implied by multiplication $\mathbf{X}_p\hat{\boldsymbol{\beta}}$ (and the other bits for the standard errors) in whatever system you specify.


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.