在R语言mgcv包gam模型中同时使用s()与I()的方法

问题重现 #

我在使用mgcv包的gam()函数拟合数据时，调用s()定义平滑项，编写了如下代码：

y =  data.frame(n = c(1:5),b = runif(5),c = runif(5,min = 2,max = 3),d = c(5:9),e = scale(c(11:15)),f = c(12:16))

model1 = gam(n ~ b + c + I(s(sqrt(d),bs = 'tp',df=2)) +
                I(s(sqrt(e),bs = 'tp',df=2)) +
                I(s(sqrt(f),bs = 'tp',df=2)),
              family=quasipoisson(),
              data = y,
              method = 'REML',
              select = TRUE)

运行后触发错误：

Error in terms.formula(reformulate(term[i])) : 
  invalid model formula in ExtractVars

需要知道正确的公式写法，以及s()与I()在公式中的使用规则。

正确写法 #

无需用I()包裹s()函数，直接在s()内部处理变量变换即可，修正后的代码：

model1 = gam(n ~ b + c + 
                s(sqrt(d), bs = 'tp', df = 2) +
                s(sqrt(e), bs = 'tp', df = 2) +
                s(sqrt(f), bs = 'tp', df = 2),
              family=quasipoisson(),
              data = y,
              method = 'REML',
              select = TRUE)

s()与I()的使用规则 #

• s()是mgcv专门用于定义平滑项的特殊函数，公式解析器会直接识别并处理它内部的变量变换逻辑，不需要额外用I()包裹
• I()的作用是强制stats::formula将括号内的内容当作纯算术表达式解析，避免公式语法（如+、~）干扰计算。但用它包裹s()会让解析器错误地将s(...)视为普通算术表达式，引发语法错误
• 只有当你需要将变换后的变量作为普通线性项加入模型时，才需要用I()，例如：n ~ b + c + I(sqrt(d))

会话信息 #

> sessionInfo()
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] Hmisc_4.8-0       Formula_1.2-5     survival_3.4-0    lattice_0.20-45   mgcv_1.8-41       nlme_3.1-160      dlnm_2.4.7        rvest_1.0.3       plotly_4.10.1     data.table_1.14.8 arrow_11.0.0.2   
[12] magrittr_2.0.3    lubridate_1.9.2   forcats_1.0.0     stringr_1.5.0     dplyr_1.1.0       purrr_1.0.1       readr_2.1.4       tidyr_1.3.0       tibble_3.1.8      ggplot2_3.4.1     tidyverse_2.0.0  
[23] readxl_1.4.2     

loaded via a namespace (and not attached):
  [1] minqa_1.2.5             colorspace_2.1-0        deldir_1.0-6            sjlabelled_1.2.0        htmlTable_2.4.1         estimability_1.4.1      snakecase_0.11.0        base64enc_0.1-3        
  [9] rstudioapi_0.14.0-9000  bit64_4.0.5             fansi_1.0.4             mvtnorm_1.1-3           xml2_1.3.3              splines_4.2.2           cachem_1.0.6            knitr_1.42             
 [17] sjmisc_2.8.9            jsonlite_1.8.4          splitstackshape_1.4.8   nloptr_2.0.3            ggeffects_1.2.1         broom_1.0.4             cluster_2.1.4           png_0.1-8               
 [25] clipr_0.8.0             compiler_4.2.2          httr_1.4.5              sjstats_0.18.2          emmeans_1.8.5           backports_1.4.1         assertthat_0.2.1        Matrix_1.5-1           
 [33] fastmap_1.1.0           lazyeval_0.2.2          cli_3.6.0               htmltools_0.5.4         tools_4.2.2             coda_0.19-4             gtable_0.3.3            glue_1.6.2              
 [41] Rcpp_1.0.10             santoku_0.9.0           cellranger_1.1.0        vctrs_0.5.2             sjPlot_2.8.12           insight_0.19.1          lemon_0.4.6             xfun_0.37               
 [49] metR_0.14.0             openxlsx_4.2.5.2        lme4_1.1-31             timechange_0.2.0        lifecycle_1.0.3         MASS_7.3-58.1           scales_1.2.1            hms_1.1.3               
 [57] RColorBrewer_1.1-3      curl_5.0.0              memoise_2.0.1           gridExtra_2.3           padr_0.6.2              rpart_4.1.19            latticeExtra_0.6-30     stringi_1.7.12         
 [65] bayestestR_0.13.0       checkmate_2.1.0         boot_1.3-28.1           zip_2.2.2               repr_1.1.6              rlang_1.1.0             pkgconfig_2.0.3         tsModel_0.6-1           
 [73] htmlwidgets_1.6.2       bit_4.0.5               tidyselect_1.2.0        plyr_1.8.8              R6_2.5.1                generics_0.1.3          foreign_0.8-83          pillar_1.9.0           
 [81] withr_2.5.0             nnet_7.3-18             datawizard_0.7.1        performance_0.10.2      janitor_2.2.0           modelr_0.1.11           crayon_1.5.2            interp_1.1-3           
 [89] utf8_1.2.3              tzdb_0.3.0              extraInserts_0.0.0.9003 viridis_0.6.2           jpeg_0.1-10             grid_4.2.2              digest_0.6.31           xtable_1.8-4           
 [97] ggeasy_0.1.4            munsell_0.5.0           viridisLite_0.4.1       skimr_2.1.5 

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