如何在R的brms包中获取样条项基函数（类mgcv包实现）

1. 指定样条类型 #

brms的公式接口完全兼容mgcv的样条语法，直接在s()中通过bs参数指定样条类型，用法和mgcv一致：

• 三次回归样条：bs = "cr"
• 薄板样条：bs = "tp"
• 自然样条：bs = "ns"

示例代码：

# 指定三次回归样条的模型公式
formula_cr <- brmsformula(wage ~ year + s(age, bs = "cr"))
# 拟合模型
fit_cr_brms <- brm(formula = formula_cr, family = gaussian(), 
                   data = Wage, seed = 1, chains = 1, iter = 2000)

2. 提取样条基函数矩阵 #

brms拟合时会将样条设计矩阵存储在模型的standata中，也可以通过专用函数提取：

• standata(fit_cr_brms)$Zs_1_1：这是中心化后的样条基矩阵，每行对应一个观测，每列对应一个样条基函数。
• 若需要和mgcv中predict(..., type = "lpmatrix")结构一致的原始基矩阵，可使用brms::get_z()：

# 方法1：从standata提取中心化基矩阵
z_matrix <- standata(fit_cr_brms)$Zs_1_1

# 方法2：用get_z获取对应样条的设计矩阵
z_list <- get_z(fit_cr_brms)
age_spline_basis <- z_list$sage  # 提取age对应的样条基矩阵

3. 绘制样条基函数 #

生成覆盖age范围的序列，构造对应的基矩阵后即可可视化每一个基函数：

library(ggplot2)
library(mgcv)
library(tidyr)

# 生成age序列
new_age <- seq(min(Wage$age), max(Wage$age), length.out = 100)
new_data <- data.frame(age = new_age, year = mean(Wage$year))  # 固定year消除干扰

# 构造与brms一致的样条基矩阵
smooth_obj <- smoothCon(s(age, bs = "cr"), data = new_data, knots = NULL)[[1]]
basis_matrix <- smooth_obj$X

# 转换为绘图格式
basis_df <- as.data.frame(basis_matrix)
colnames(basis_df) <- paste0("Basis_", 1:ncol(basis_matrix))
basis_df$age <- new_age
basis_long <- pivot_longer(basis_df, cols = starts_with("Basis_"), 
                           names_to = "Basis", values_to = "Value")

# 绘制基函数
ggplot(basis_long, aes(x = age, y = Value, color = Basis)) +
  geom_line() +
  theme_minimal() +
  labs(title = "三次回归样条基函数", x = "Age", y = "基函数值")

4. 获取特定age值的基函数值 #

构造包含目标age的新数据，生成对应基矩阵即可提取值：

# 目标age值
target_ages <- c(5, 120)
target_data <- data.frame(age = target_ages, year = mean(Wage$year))

# 构造基矩阵
smooth_obj_target <- smoothCon(s(age, bs = "cr"), data = target_data, knots = NULL)[[1]]
target_basis <- smooth_obj_target$X

# 查看结果
rownames(target_basis) <- target_ages
print(target_basis)

内容的提问来源于stack exchange，提问作者Kana