Извлечение сплайнов из модели класса GAM (`mgcv::gam`)

#r

                    
Прошу подсказать по следующему вопросу, который является продолжением данного вопроса.
Строим аддитивную модель, как показано на примере из справки ?predict.gam:

 library(mgcv)
 n <- 200
 sig <- 2
 dat <- gamSim(1,n=n,scale=sig)

 b <- gam(y ~ s(x0) + s(I(x1^2)) + s(x2) + offset(x3), data = dat)

 newd <- data.frame(x0=(0:30)/30, x1=(0:30)/30, x2=(0:30)/30, x3=(0:30)/30)

 Xp <- predict(b, newd, type="lpmatrix")

 ##################################################################
 ## The following shows how to use use an "lpmatrix" as a lookup 
 ## table for approximate prediction. The idea is to create 
 ## approximate prediction matrix rows by appropriate linear 
 ## interpolation of an existing prediction matrix. The additivity 
 ## of a GAM makes this possible. 
 ## There is no reason to ever do this in R, but the following 
 ## code provides a useful template for predicting from a fitted 
 ## gam *outside* R: all that is needed is the coefficient vector 
 ## and the prediction matrix. Use larger `Xp'/ smaller `dx' and/or 
 ## higher order interpolation for higher accuracy.  
 ###################################################################

 xn <- c(.341,.122,.476,.981) ## want prediction at these values
 x0 <- 1         ## intercept column
 dx <- 1/30      ## covariate spacing in `newd'
 for (j in 0:2) { ## loop through smooth terms
   cols <- 1+j*9 +1:9      ## relevant cols of Xp
   i <- floor(xn[j+1]*30)  ## find relevant rows of Xp
   w1 <- (xn[j+1]-i*dx)/dx ## interpolation weights
   ## find approx. predict matrix row portion, by interpolation
   x0 <- c(x0,Xp[i+2,cols]*w1 + Xp[i+1,cols]*(1-w1))
 }
 dim(x0)<-c(1,28) 
 fv <- x0%*%coef(b) + xn[4];fv    ## evaluate and add offset
 se <- sqrt(x0%*%b$Vp%*%t(x0));se ## get standard error
 ## compare to normal prediction
 predict(b,newdata=data.frame(x0=xn[1],x1=xn[2],
         x2=xn[3],x3=xn[4]),se=TRUE)


Возможно ли каким-то образом извлечь из модели в явном виде саму формулу сплайна
s(x), которая бы соответствовала форме записи:

yi = β0 + β1 b1 (xi) + β2 b2 (xi) + · · · + βK+3 bK+3(xi) 

(James G. et al. - An Introduction to Statistical Learning with Applications in R)

Заметил, что путем умножения коэффициентов coef(mod) на predict(mod, type="lpmatrix")
(матрица модели) можно получить предсказанные значения, возвращаемые функцией predict(mod,
type="response"). Собственно, проблема в том, что я не могу поставить эти коэффициенты
в соответствие степеням переменной x и (x-knot), как это предполагается при записи
сплайна в виде комбинации базисных функций. 
Приведенный пример предсказания значений на основе интерполяции значений из "lpmatrix"
служит, согласно справке, для использования модели вне среды R. Означает ли это, что
данная реализация GAM не предполагает получения записи модели в явном виде? Есть ли
отличия в этом плане в функции gam() из пакета gam, написанного Хасти и Тибширани -
создателями методам обобщенных аддитивных моделей?
Спасибо.
    

Ответы

Ответ 1

В пакете rms, который является приложением к известной книге Ф. Харрелла, есть функция
Function(), которая выдает уравнение поданной на нее модели. К сожалению, объекты класса
gam эта функция не принимает. Но в состав rms входят другие функции, которые позволяют
подгонять сплайн-модели - возможно, они подойдут и для Ваших целей. Примеры можно посмотреть
здесь и здесь.