Как извлечь подогнанные сплайны из GAM (`mgcv::gam`)

Question

Я использую GAM для моделирования временных трендов в логистической регрессии. Тем не менее, я хотел бы извлечь из него встроенный сплайн, чтобы добавить его к другой модели, которая не может быть встроена в GAM или GAMM.

Таким образом, у меня есть 2 вопроса:

1. Как я могу установить сглаживатель со временем, чтобы один узел находился в определенном месте, позволяя модели находить другие узлы?

2. Как я могу извлечь матрицу из встроенного GAM, чтобы я мог использовать ее в качестве вменения для другой модели?

Типы моделей, которые я запускаю, имеют следующую форму:

gam <- gam(mortality.under.2~ maternal_age_c+ I(maternal_age_c^2)+
           s(birth_year,by=wealth2) + wealth2 + sex +
           residence + maternal_educ + birth_order,
           data=colombia2, family="binomial")

Я прочитал обширную документацию для GAM, но я все еще не уверен. Любое предложение действительно ценится.

Answer 1

В mgcv::gam есть способ сделать это (ваш Q2) через predict.gam метод и type = "lpmatrix",

?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)

Это проходит весь процесс, даже шаг прогнозирования, который будет сделан за пределами R или модели GAM. Вам нужно будет немного изменить пример, чтобы сделать то, что вы хотите, так как пример оценивает все термины в модели, и у вас есть два других термина, кроме сплайна - по сути, вы делаете то же самое, но только для сплайн-терминов, которые включает в себя поиск соответствующих столбцов и строк Xp матрица для сплайна. Кроме того, вы также должны заметить, что сплайн отцентрирован, так что вы можете или не хотите отменить это тоже.

Для вашего Q1 выберите соответствующие значения для xn вектор / матрица в примере. Они соответствуют значениям для nчлен в модели. Поэтому установите те, которые вы хотите зафиксировать, на некоторое среднее значение, а затем измените значение, связанное со сплайном.

Если вы делаете все это в R, было бы проще просто оценить сплайн по значениям ковариации сплайна, для которой у вас есть данные, которые входят в другую модель. Вы делаете это, создавая фрейм данных значений, на которых можно прогнозировать, а затем использовать

predict(mod, newdata = newdat, type = "terms")

где mod это установленная модель GAM (через mgcv::gam), newdat является фреймом данных, содержащим столбец для каждой переменной в модели (включая параметрические термины; установите термины, которые вы не хотите изменять, на некоторое постоянное среднее значение [скажем, среднее значение переменной в наборе данных] или определенный уровень, если фактор). type = "terms" часть вернет матрицу для каждой строки в newdat с "вкладом" в подогнанное значение для каждого члена в модели, включая сплайн-член. Просто возьмите столбец этой матрицы, который соответствует сплайну - снова он центрирован.

Возможно, я неправильно понял ваш Q1. Если вы хотите контролировать узлы, см. knots аргумент mgcv::gam, По умолчанию, mgcv::gam помещает узел в крайности данных, а затем оставшиеся "узлы" распределяются равномерно по интервалу. mgcv::gam не находит узлы - он помещает их для вас, и вы можете контролировать, где он размещает их через knots аргумент.