如何从ggplot的geom_smooth(method='gam')提取模型参数？

问题描述 #

我在ggplot绘图中使用了geom_smooth(method='gam')函数，想了解是否可以提取系数、偏差解释等模型参数，要求输出结果与mgcv::gam()函数summary()的输出格式一致。

测试数据 #

test <- structure(list(numDate = c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 
5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 
11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 
16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 
22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 27, 
27, 27, 28, 28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 31, 32, 32, 
32, 33, 33, 33, 34), mAODscale = c(0.0388813764134284, 0.0461877433707656, 
0.096886995092774, 0.0371438764134382, 0.0394419100374392, 0.0893994950927635, 
0.0413855430800965, 0.0418085767040992, 0.101328661759439, 0.0428313764134316, 
0.0403835767041016, 0.0978328284261067, 0.0444813764134295, 0.0566460767040979, 
0.136811995092771, 0.0404647097467716, 0.0461127433707702, 0.104611995092768, 
0.0391105430800991, 0.0388294100374367, 0.0883828284261057, 0.0377355430801032, 
0.0334335767040983, 0.0744786617594428, 0.0346647097467638, 0.0316752433707705, 
0.0691911617594343, 0.0365730430800966, 0.0329127433707725, 0.0721203284261094, 
0.0337897097467703, 0.0271960767041008, 0.0568703284261005, 0.0321188764134348, 
0.0226835767040967, 0.0467286617594311, 0.0389522097467676, 0.0317169100374315, 
0.0724911617594444, 0.0374147097467699, 0.0301335767041024, 0.066378661759444, 
0.0359688764134347, 0.0245710767041061, 0.0492828284261009, 0.0340355430801083, 
0.0208835767040938, 0.0401828284261114, 0.033193876413435, 0.0170627433707722, 
0.0363119950927739, 0.0320272097467722, 0.0147294100374324, 0.0298161617594417, 
0.0305563764134433, 0.0123752433707693, 0.0200744950927714, 0.0294522097467649, 
0.00966691003743847, 0.0144453284261061, 0.029439709746768, 0.00845441003743019, 
0.0127494950927769, 0.029218876413438, 0.00767941003742578, 0.0118203284260971, 
0.0283438764134303, 0.00608774337077023, 0.00807032842610056, 
0.027393235387791, 0.00498582029383954, 0.00612513611841337, 
0.0261313764134314, 0.004612743370771, 0.00541616175944171, 0.0260813764134298, 
0.00472524337077118, 0.00372449509276862, 0.0256522097467666, 
0.00535024337077061, 0.00356616175943714, 0.0250313764134376, 
0.0161044100374284, 0.00060366175944182, 0.0241147097467689, 
0.025246076704093, -0.00540050490722876, 0.0233022097467739, 
0.022454410037426, -0.00630883824055672, 0.0208772097467715, 
0.0152252433707645, -0.00819217157389573, 0.0194897097467646, 
-0.00307475662923196, -0.00896300490722979, 0.0178938764134386, 
0.00323774337076088, -0.00836717157389444, 0.0153022097467641
)), row.names = c(NA, -100L), class = c("data.table", "data.frame"
), .internal.selfref = <pointer: 0x000002cacef91f60>)

绘图代码 #

p <- ggplot(test, aes(x=numDate, y=mAODscale)) +
  geom_point() +
  stat_smooth(method='gam',
              formula=y~s(x,bs="cs",fx=TRUE,k=10))

解决方案 #

当然可以提取目标参数。ggplot运行stat_smooth(method='gam')时，会在内部完成GAM模型的拟合，你只需按以下步骤提取模型并生成标准输出：

1. 从ggplot对象中提取拟合好的GAM模型 #

使用ggplot_build()解析绘图对象，然后取出存储在其中的GAM模型：

# 解析绘图对象
p_build <- ggplot_build(p)
# 提取模型
gam_model <- p_build$data[[1]]$model

2. 生成与mgcv::gam()一致的summary输出 #

直接对提取到的模型调用summary()函数，输出格式和直接用mgcv::gam()拟合后调用summary的结果完全一致：

summary(gam_model)

验证一致性 #

你可以对比直接拟合的模型和从ggplot中提取的模型的summary，结果完全相同：

# 直接用mgcv拟合相同模型
direct_gam <- mgcv::gam(mAODscale ~ s(numDate, bs="cs", fx=TRUE, k=10), data=test)
summary(direct_gam)

# 对比从ggplot提取的模型summary
summary(gam_model)

通过这种方式，你就能获取到系数、偏差解释、显著性检验等所有标准GAM模型参数。

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