32 广义可加模型

相比于广义线性模型，广义可加模型可以看作是一种非线性模型，模型中含有非线性的成分。

注释

多元适应性（自适应）回归样条 multivariate adaptive regression splines
Friedman, Jerome H. 1991. Multivariate Adaptive Regression Splines. The Annals of Statistics. 19(1):1–67. https://doi.org/10.1214/aos/1176347963
earth: Multivariate Adaptive Regression Splines http://www.milbo.users.sonic.net/earth
Friedman, Jerome H. 2001. Greedy function approximation: A gradient boosting machine. The Annals of Statistics. 29(5):1189–1232. https://doi.org/10.1214/aos/1013203451
Friedman, Jerome H., Trevor Hastie and Robert Tibshirani. Additive Logistic Regression: A Statistical View of Boosting. The Annals of Statistics. 28(2): 337–374. http://www.jstor.org/stable/2674028
Flexible Modeling of Alzheimer’s Disease Progression with I-Splines PDF 文档
Implementation of B-Splines in Stan 网页文档

32.1 频率派

MASS 包的 mcycle 数据集

data(mcycle, package = "MASS")
library(ggplot2)
ggplot(data = mcycle, aes(x = times, y = accel)) +
  geom_point() +
  theme_bw()

样条回归

library(mgcv)
fgam <- gam(accel ~ s(times), data = mcycle, method = "REML")
summary(fgam)

#> 
#> Family: gaussian 
#> Link function: identity 
#> 
#> Formula:
#> accel ~ s(times)
#> 
#> Parametric coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  -25.546      1.951  -13.09   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>            edf Ref.df    F p-value    
#> s(times) 8.625  8.958 53.4  <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.783   Deviance explained = 79.7%
#> -REML = 616.14  Scale est. = 506.35    n = 133

方差成分

gam.vcomp(fgam, rescale = FALSE)

#> 
#> Standard deviations and 0.95 confidence intervals:
#> 
#>            std.dev     lower      upper
#> s(times) 807.88726 480.66162 1357.88215
#> scale     22.50229  19.85734   25.49954
#> 
#> Rank: 2/2

plot(fgam)

32.2 贝叶斯派

library(cmdstanr)

bgam <- brms::brm(accel ~ s(times),
  data = mcycle, family = gaussian(), cores = 2, seed = 20232023,
  iter = 4000, warmup = 1000, thin = 10, refresh = 0,
  control = list(adapt_delta = 0.99)
)
summary(bgam)

brms::fixef(bgam) # 固定效应
# brms::ranef(bgam) # 随机效应

ms_bgam <- brms::conditional_smooths(bgam) # 平滑效应
plot(ms_bgam)

brms::bayes_R2(bgam)
brms::waic(bgam)
# brms::log_lik(bgam)
brms::loo(bgam)

brms::pp_check(bgam, ndraws = 50)

brms::make_stancode(accel ~ s(times), data = mcycle, family = gaussian())