У меня нет ваших данных, поэтому я беру следующий пример из ?adaptive.smooth
, чтобы показать вам, где вы можете найти нужную информацию. Обратите внимание, что хотя этот пример предназначен для гауссовских данных, а не для данных Пуассона, отличается только функция связи; все остальное просто стандартное.
x <- 1:1000/1000 # data between [0, 1]
mu <- exp(-400*(x-.6)^2)+5*exp(-500*(x-.75)^2)/3+2*exp(-500*(x-.9)^2)
y <- mu+0.5*rnorm(1000)
b <- gam(y~s(x,bs="ad",k=40,m=5))
Теперь вся информация о плавном строительстве хранится в b$smooth
, достаем:
smooth <- b$smooth[[1]] ## extract smooth object for first smooth term
узлы:
smooth$knots
показывает расположение узлов.
> smooth$knots
[1] -0.081161 -0.054107 -0.027053 0.000001 0.027055 0.054109 0.081163
[8] 0.108217 0.135271 0.162325 0.189379 0.216433 0.243487 0.270541
[15] 0.297595 0.324649 0.351703 0.378757 0.405811 0.432865 0.459919
[22] 0.486973 0.514027 0.541081 0.568135 0.595189 0.622243 0.649297
[29] 0.676351 0.703405 0.730459 0.757513 0.784567 0.811621 0.838675
[36] 0.865729 0.892783 0.919837 0.946891 0.973945 1.000999 1.028053
[43] 1.055107 1.082161
Обратите внимание, три внешних узла размещены с каждой стороны [0, 1]
для построения основы сплайна.
базовый класс
attr(smooth, "class")
сообщает вам тип сплайна. Как вы можете прочитать из ?adaptive.smooth
, для bs = ad
, mgcv
используйте P-шлицы, поэтому вы получите "pspline.smooth".
mgcv
используйте pspline 2-го порядка, вы можете проверить это, проверив матрицу различий smooth$D
. Ниже приведен снимок:
> smooth$D[1:6,1:6]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 -2 1 0 0 0
[2,] 0 1 -2 1 0 0
[3,] 0 0 1 -2 1 0
[4,] 0 0 0 1 -2 1
[5,] 0 0 0 0 1 -2
[6,] 0 0 0 0 0 1
коэффициенты
Вы уже знаете, что b$coefficients
содержат коэффициенты модели:
beta <- b$coefficients
Обратите внимание, что это именованный вектор:
> beta
(Intercept) s(x).1 s(x).2 s(x).3 s(x).4 s(x).5
0.37792619 -0.33500685 -0.30943814 -0.30908847 -0.31141148 -0.31373448
s(x).6 s(x).7 s(x).8 s(x).9 s(x).10 s(x).11
-0.31605749 -0.31838050 -0.32070350 -0.32302651 -0.32534952 -0.32767252
s(x).12 s(x).13 s(x).14 s(x).15 s(x).16 s(x).17
-0.32999553 -0.33231853 -0.33464154 -0.33696455 -0.33928755 -0.34161055
s(x).18 s(x).19 s(x).20 s(x).21 s(x).22 s(x).23
-0.34393354 -0.34625650 -0.34857906 -0.05057041 0.48319491 0.77251118
s(x).24 s(x).25 s(x).26 s(x).27 s(x).28 s(x).29
0.49825345 0.09540020 -0.18950763 0.16117012 1.10141701 1.31089436
s(x).30 s(x).31 s(x).32 s(x).33 s(x).34 s(x).35
0.62742937 -0.23435309 -0.19127140 0.79615752 1.85600016 1.55794576
s(x).36 s(x).37 s(x).38 s(x).39
0.40890236 -0.20731309 -0.47246357 -0.44855437
базовая матрица / матрица модели / матрица линейных предикторов (lpmatrix)
Вы можете получить матрицу модели из:
mat <- predict.gam(b, type = "lpmatrix")
Это матрица n-by-p
, где n
- количество наблюдений, а p
- количество коэффициентов. Эта матрица имеет имя столбца:
> head(mat[,1:5])
(Intercept) s(x).1 s(x).2 s(x).3 s(x).4
1 1 0.6465774 0.1490613 -0.03843899 -0.03844738
2 1 0.6437580 0.1715691 -0.03612433 -0.03619157
3 1 0.6384074 0.1949416 -0.03391686 -0.03414389
4 1 0.6306815 0.2190356 -0.03175713 -0.03229541
5 1 0.6207361 0.2437083 -0.02958570 -0.03063719
6 1 0.6087272 0.2688168 -0.02734314 -0.02916029
В первом столбце все 1, что дает перехват. В то время как s(x).1
предлагает первую базовую функцию для s(x)
. Если вы хотите увидеть, как выглядит отдельная базовая функция, вы можете построить столбец mat
для вашей переменной. Например:
plot(x, mat[, "s(x).20"], type = "l", main = "20th basis")
линейный предсказатель
Если вы хотите создать подгонку вручную, вы можете:
pred.linear <- mat %*% beta
Обратите внимание, что это именно то, что вы можете получить от b$linear.predictors
или
predict.gam(b, type = "link")
ответ / подогнанные значения
Для негауссовских данных, если вы хотите получить переменную ответа, вы можете применить функцию обратной связи к линейному предиктору, чтобы вернуться к исходному масштабу.
Семейная информация хранится в gamObject$family
, а gamObject$family$linkinv
- это функция обратной связи. Приведенный выше пример обязательно даст вам идентификационную ссылку, но для вашего подогнанного объекта x.gam
вы можете сделать:
x.gam$family$linkinv(x.gam$linear.predictors)
Обратите внимание: это то же самое, что и x.gam$fitted
, или
predict.gam(x.gam, type = "response").
Другие ссылки
Я только что понял, что подобных вопросов было довольно много.
- Этот ответ Гэвина Симпсона отлично подходит для
predict.gam( , type = 'lpmatrix')
.
- Этот ответ касается
predict.gam(, type = 'terms')
.
Но в любом случае лучший справочник - это ?predict.gam
, который включает обширные примеры.
person
Zheyuan Li
schedule
22.05.2016