Jak wyodrębnić dopasowane splajmy z GAM (`mgcv :: gam`)

Question

Używam GAM do modelowania trendów czasowych w regresji logistycznej. Chciałbym jednak wyciągnąć z niego dopasowany splajn, aby dodać go do innego modelu, którego nie można zamontować w GAM lub GAMM.

Zatem mam 2 pytania:

1. Jak mogę dopasować gładsza upływem czasu, tak aby wymusić jeden węzeł będzie w określonym miejscu, a jednocześnie przepuszczają model, aby znaleźć inne węzły?

2. W jaki sposób mogę wyodrębnić matrycę z dopasowanego GAM, aby można go było wykorzystać jako przypięcie dla innego modelu?

Rodzaje modeli używam mają następującą postać:

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

czytałem obszerną dokumentację dla GAM, ale nie jestem pewien jeszcze. Każda sugestia jest naprawdę doceniana.

Answer 1

W mgcv::gam jest sposób, aby to zrobić (twój Q2), za pomocą metody predict.gam i type = "lpmatrix".

?predict.gam ma nawet przykład, który mi reprodukowania poniżej:

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) 

który przechodzi przez cały proces nawet kroku przewidywania, które byłyby wykonane poza R lub modelu GAM. Będziesz musiał nieco zmodyfikować przykład, aby zrobić to, co chcesz, ponieważ przykład ocenia wszystkie terminy w modelu, a oprócz splajnu masz dwa inne terminy - w zasadzie robisz to samo, ale tylko dla terminów splajnu, które polega na znalezieniu odpowiednich kolumn i rzędów matrycy Xp dla splajnu. Następnie należy pamiętać, że splajn jest wyśrodkowany, więc możesz lub nie chcesz tego cofnąć.

Dla swojego Q1, wybierz odpowiednie wartości dla wektora/macierzy xn w przykładzie. Odpowiadają one wartościom z modelu n th. Ustaw te, które chcesz ustawić, na pewną średnią wartość, a następnie zmień tę, która jest powiązana z splajnem.

Jeśli robisz to wszystko w R, łatwiej będzie po prostu ocenić splajn na wartości współzmiennej splajnu, dla której masz dane, które przechodzą do innego modelu. To zrobić tworząc ramkę danych wartości, w których do przewidzenia, a następnie użyć

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

gdzie mod jest wyposażona modelu GAM (przez mgcv::gam) newdat jest ramka danych zawiera kolumny dla każdej zmiennej w model (w tym terminy parametryczne, ustaw terminy, których nie chcesz zmieniać, do pewnej stałej wartości średniej [podaj średnią zmiennej w zestawie danych] lub pewien poziom, jeśli czynnik). Część type = "terms" zwróci matrycę dla każdego wiersza w newdat z "wkładem" do dopasowanej wartości dla każdego terminu w modelu, włączając w to splajn. Po prostu weź kolumnę tej matrycy, która odpowiada splajnowi - znowu jest wyśrodkowana.

Być może źle zrozumiałem twój Q1. Jeśli chcesz kontrolować węzły, zobacz argument knots dla mgcv::gam.Domyślnie mgcv::gam umieszcza węzeł na krańcach danych, a pozostałe "węzły" rozkładają się równomiernie w przedziale. mgcv::gam nie szuka znaleźć węzłów - umieszcza je dla Ciebie i można kontrolować, gdzie umieszcza je za pośrednictwem argumentu knots.