如何访问和重用 `R` 中的 `mgcv` 包中的 smooths？
How to access and reuse the smooths in the `mgcv` package in `R`?
我正在查看 R 中的 mgcv 包，我想知道
如何根据新数据更新模型。例如，假设我有
以下数据，我有兴趣拟合三次回归样条。
# Load library.
library(mgcv)
# Set seed.
set.seed(2022)
# Data
x <- seq(10, 100, by = 10)
y <- sort(runif(10))
我可以使用 mgcv::s() 函数拟合模型来转换我的预测变量，其中 bs = "cr" 代表 三次回归样条 ，如
文档（即 ?mgcv::s）。
# Fit.
model <- mgcv::gam(y ~ s(x, bs = "cr"))
# Print model.
model
# Family: gaussian
# Link function: identity
#
# Formula:
# y ~ s(x, bs = "cr")
#
# Estimated degrees of freedom:
# 7.51 total = 8.51
#
# GCV score: 0.001123237
我假设 mgcv::s() 习惯于
确定样条基函数的内结？如果我想
插值 x 的整个范围，看来我可以使用 predict
功能。
# Prepare range of `x` for interpolation.
x_new <- 10:100
# Interpolate.
mgcv_interpolation <- predict(model, type = "link", newdata = data.frame(x = x_new))
# Plot.
plot(x, y, pch = 19)
lines(x_new, mgcv_interpolation, lwd = 2, col = "red")
我不清楚的是当有新数据时如何更新模型（即y）
进来。例如，假设我的新数据看起来像这样。
# Set seed.
set.seed(2022)
y_new <- sort(sample(y, size = length(y), replace = TRUE))
我的理解是我可以简单地使用之前创建的基础矩阵，
但我不确定如何用 mgcv 做到这一点。例如，这就是我能做的
它手动使用 B 样条基础。
# ...
# Suppose that based on some cross-validation procedure `df = 6` is selected.
df <- 6
# Create B-Spline basis functions.
basis <- splines::bs(x, df = df, degree = 3, intercept = TRUE)
# Estimate spline coefficients.
coefficients <- lm.fit(basis, y)$coef
# Compute fitted values.
fitted <- basis %*% coefficients
# Create extended basis for `x_new`.
basis_x_new <- splines::bs(x_new, df = df, degree = 3, intercept = TRUE)
# Interpolate.
bs_interpolation <- basis_x_new %*% coefficients
# Add to previous plot.
lines(x_new, bs_interpolation, lwd = 2, col = "blue")
# Update model based on `y_new`.
coefficients_y_new <- lm.fit(basis, y_new)$coef
# Add points and lines to the previous plot.
points(x, y_new, pch = 19, col = "orange")
lines(x_new, basis_x_new %*% coefficients_y_new, lwd = 2, col = "orange")
我想我的问题是如何找到 mgcv::s() 创建和
在 mgcv::gam 的后续调用中重用它？或者，还有更多
mgcv-这样做的惯用方式？
编辑 1.
仔细研究后，我发现我可以使用带有参数 type = "lpmatrix" 的 mgcv::predict.gam() 来提取基矩阵。但是，我仍然无法复制 mgcv::gam() 提供的确切系数。差异不大，但我想知道它们来自哪里。例如：
# Extract the basis matrix from the `gam` object.
basis_gam <- mgcv::predict.gam(model, type = "lpmatrix")
# Fit the model using the basis matrix.
model_basis_gam <- mgcv::gam(y ~ basis_gam - 1)
# Compare the coefficients.
round(data.frame(
difference = coef(model) - coef(model_basis_gam)
), 4)
# difference
# (Intercept) 0.0000
# s(x).1 0.0004
# s(x).2 -0.0033
# s(x).3 0.0071
# s(x).4 -0.0020
# s(x).5 -0.0016
# s(x).6 -0.0054
# s(x).7 0.0103
# s(x).8 -0.0064
# s(x).9 0.0017
编辑 2.
似乎有一个函数 mgcv::bam.update() 可以为新日期更新 GAM 模型，但对于通过 mgcv::bam() 而不是 mgcv::gam() 拟合的模型。尽管如此，S3 方法 update 似乎与 mgcv::gam() 对象一起工作，可能是因为 class(model) 包含 "gam" "glm" "lm"，但是，在文档。例如：
# Update the model for `y_new`.
model_y_new_via_update <- update(model, data = data.frame(y = y_new))
# Extract the basis matrices for `model` and `model_y_new_via_update`.
basis_model <- mgcv::predict.gam(model, type = "lpmatrix")
basis_model_y_new_via_update <- mgcv::predict.gam(model_y_new_via_update, type = "lpmatrix")
# Check that both models used the same basis matrix.
all(basis_model == basis_model_y_new_via_update)
# TRUE
此外，还有一些我无法解释的系数差异。
# Fit the model using the extracted basis matrix.
model_y_new_via_basis <- mgcv::gam(y_new ~ basis_model - 1)
# Eyeballing the coefficients.
round(data.frame(
via_update = coef(model_y_new_via_update),
via_basis = coef(model_y_new_via_basis),
difference = coef(model_y_new_via_update) - coef(model_y_new_via_basis),
row.names = names(coef(model))
), 4)
# via_update via_basis difference
# (Intercept) 0.4420 0.4420 0.0000
# s(x).1 -0.2385 -0.2333 -0.0052
# s(x).2 -0.1901 -0.1689 -0.0212
# s(x).3 -0.0854 -0.1689 0.0835
# s(x).4 0.1315 0.1902 -0.0586
# s(x).5 0.2666 0.2821 -0.0155
# s(x).6 0.2907 0.2821 0.0085
# s(x).7 0.2855 0.2821 0.0033
# s(x).8 0.3119 0.2936 0.0183
# s(x).9 0.3917 0.4036 -0.0120
这是一个简单的例子
- 创建您的
smoothCon 对象，使用 x
sm = smoothCon(s(x, bs="cr"), data=data.frame(x))[[1]]
- 创建简单函数以获取给定
y 和您的 smoothCon 对象的 beta 系数
get_beta <- function(y,sm) {
as.numeric(coef(lm(y~sm$X-1)))
}
- 给定
x、y 和 smoothCon 对象 ，创建简单的函数来获取预测
get_pred <- function(x,y,sm) {
PredictMat(sm, data.frame(x=x)) %*% get_beta(y, sm)
}
- 用红色绘制原始 x,y 点，用蓝色绘制新的 x,y 点
plot(x,y, col="red")
points(x,y_new, col="blue")
- 添加线条，仅使用新的 x 范围 (
x_new)、旧的 (y) 和新的 (y_new) y 值以及 smoothCon 对象
lines(x_new, get_pred(x_new,y, sm), col="red")
lines(x_new, get_pred(x_new,y_new, sm), col="blue")
我正在查看 R 中的 mgcv 包，我想知道
如何根据新数据更新模型。例如，假设我有
以下数据，我有兴趣拟合三次回归样条。
# Load library.
library(mgcv)
# Set seed.
set.seed(2022)
# Data
x <- seq(10, 100, by = 10)
y <- sort(runif(10))
我可以使用 mgcv::s() 函数拟合模型来转换我的预测变量，其中 bs = "cr" 代表 三次回归样条 ，如
文档（即 ?mgcv::s）。
# Fit.
model <- mgcv::gam(y ~ s(x, bs = "cr"))
# Print model.
model
# Family: gaussian
# Link function: identity
#
# Formula:
# y ~ s(x, bs = "cr")
#
# Estimated degrees of freedom:
# 7.51 total = 8.51
#
# GCV score: 0.001123237
我假设 mgcv::s() 习惯于
确定样条基函数的内结？如果我想
插值 x 的整个范围，看来我可以使用 predict
功能。
# Prepare range of `x` for interpolation.
x_new <- 10:100
# Interpolate.
mgcv_interpolation <- predict(model, type = "link", newdata = data.frame(x = x_new))
# Plot.
plot(x, y, pch = 19)
lines(x_new, mgcv_interpolation, lwd = 2, col = "red")
我不清楚的是当有新数据时如何更新模型（即y）
进来。例如，假设我的新数据看起来像这样。
# Set seed.
set.seed(2022)
y_new <- sort(sample(y, size = length(y), replace = TRUE))
我的理解是我可以简单地使用之前创建的基础矩阵，
但我不确定如何用 mgcv 做到这一点。例如，这就是我能做的
它手动使用 B 样条基础。
# ...
# Suppose that based on some cross-validation procedure `df = 6` is selected.
df <- 6
# Create B-Spline basis functions.
basis <- splines::bs(x, df = df, degree = 3, intercept = TRUE)
# Estimate spline coefficients.
coefficients <- lm.fit(basis, y)$coef
# Compute fitted values.
fitted <- basis %*% coefficients
# Create extended basis for `x_new`.
basis_x_new <- splines::bs(x_new, df = df, degree = 3, intercept = TRUE)
# Interpolate.
bs_interpolation <- basis_x_new %*% coefficients
# Add to previous plot.
lines(x_new, bs_interpolation, lwd = 2, col = "blue")
# Update model based on `y_new`.
coefficients_y_new <- lm.fit(basis, y_new)$coef
# Add points and lines to the previous plot.
points(x, y_new, pch = 19, col = "orange")
lines(x_new, basis_x_new %*% coefficients_y_new, lwd = 2, col = "orange")
我想我的问题是如何找到 mgcv::s() 创建和
在 mgcv::gam 的后续调用中重用它？或者，还有更多
mgcv-这样做的惯用方式？
编辑 1.
仔细研究后，我发现我可以使用带有参数 type = "lpmatrix" 的 mgcv::predict.gam() 来提取基矩阵。但是，我仍然无法复制 mgcv::gam() 提供的确切系数。差异不大，但我想知道它们来自哪里。例如：
# Extract the basis matrix from the `gam` object.
basis_gam <- mgcv::predict.gam(model, type = "lpmatrix")
# Fit the model using the basis matrix.
model_basis_gam <- mgcv::gam(y ~ basis_gam - 1)
# Compare the coefficients.
round(data.frame(
difference = coef(model) - coef(model_basis_gam)
), 4)
# difference
# (Intercept) 0.0000
# s(x).1 0.0004
# s(x).2 -0.0033
# s(x).3 0.0071
# s(x).4 -0.0020
# s(x).5 -0.0016
# s(x).6 -0.0054
# s(x).7 0.0103
# s(x).8 -0.0064
# s(x).9 0.0017
编辑 2.
似乎有一个函数 mgcv::bam.update() 可以为新日期更新 GAM 模型，但对于通过 mgcv::bam() 而不是 mgcv::gam() 拟合的模型。尽管如此，S3 方法 update 似乎与 mgcv::gam() 对象一起工作，可能是因为 class(model) 包含 "gam" "glm" "lm"，但是，在文档。例如：
# Update the model for `y_new`.
model_y_new_via_update <- update(model, data = data.frame(y = y_new))
# Extract the basis matrices for `model` and `model_y_new_via_update`.
basis_model <- mgcv::predict.gam(model, type = "lpmatrix")
basis_model_y_new_via_update <- mgcv::predict.gam(model_y_new_via_update, type = "lpmatrix")
# Check that both models used the same basis matrix.
all(basis_model == basis_model_y_new_via_update)
# TRUE
此外，还有一些我无法解释的系数差异。
# Fit the model using the extracted basis matrix.
model_y_new_via_basis <- mgcv::gam(y_new ~ basis_model - 1)
# Eyeballing the coefficients.
round(data.frame(
via_update = coef(model_y_new_via_update),
via_basis = coef(model_y_new_via_basis),
difference = coef(model_y_new_via_update) - coef(model_y_new_via_basis),
row.names = names(coef(model))
), 4)
# via_update via_basis difference
# (Intercept) 0.4420 0.4420 0.0000
# s(x).1 -0.2385 -0.2333 -0.0052
# s(x).2 -0.1901 -0.1689 -0.0212
# s(x).3 -0.0854 -0.1689 0.0835
# s(x).4 0.1315 0.1902 -0.0586
# s(x).5 0.2666 0.2821 -0.0155
# s(x).6 0.2907 0.2821 0.0085
# s(x).7 0.2855 0.2821 0.0033
# s(x).8 0.3119 0.2936 0.0183
# s(x).9 0.3917 0.4036 -0.0120
这是一个简单的例子
- 创建您的
smoothCon对象，使用x
sm = smoothCon(s(x, bs="cr"), data=data.frame(x))[[1]]
- 创建简单函数以获取给定
y和您的smoothCon对象的 beta 系数
get_beta <- function(y,sm) {
as.numeric(coef(lm(y~sm$X-1)))
}
- 给定
x、y和smoothCon对象 ，创建简单的函数来获取预测
get_pred <- function(x,y,sm) {
PredictMat(sm, data.frame(x=x)) %*% get_beta(y, sm)
}
- 用红色绘制原始 x,y 点，用蓝色绘制新的 x,y 点
plot(x,y, col="red")
points(x,y_new, col="blue")
- 添加线条，仅使用新的 x 范围 (
x_new)、旧的 (y) 和新的 (y_new) y 值以及smoothCon对象
lines(x_new, get_pred(x_new,y, sm), col="red")
lines(x_new, get_pred(x_new,y_new, sm), col="blue")