关于线性模型扩展：非线性回归方法

超越线性 (Moving Beyond Linearity)

本章内容继续放宽线性假设。
第六章主要是改进，这章主要是扩展问题，有下面的部分：

• 多项式回归 (Polynomial regression) 通过添加额外的预测变量来扩展线性模型，这些预测变量是通过将每个原始预测变量进行幂运算得到的。例如，一个三次回归使用三个变量 $X$, $X^2$ 和 $X^3$ 作为预测变量。这种方法提供了一种简单的方式来对数据进行非线性拟合。
• 阶梯函数 (Step functions) 将一个变量的范围切割成 K 个不同的区域，以产生一个定性变量。这相当于拟合一个分段常数函数。
• 回归样条 (Regression splines) 比多项式和阶梯函数更灵活，实际上是这两者的扩展。它们涉及将 $X$ 的范围划分为 K 个不同的区域。在每个区域内，对数据拟合一个多项式函数。然而，这些多项式受到约束，以便它们在区域边界（或称节点 (knots)）处平滑连接。只要区间被划分为足够多的区域，这就可以产生一个极其灵活的拟合。
• 平滑样条 (Smoothing splines) 类似于回归样条，但在略有不同的情况下出现。平滑样条是通过最小化一个带有平滑度惩罚项的残差平方和准则而得到的。
• 局部回归 (Local regression) 类似于样条，但在一个重要方面有所不同。其区域被允许重叠，并且它们确实以一种非常平滑的方式重叠。
• 广义可加模型 (Generalized additive models) 允许我们扩展上述方法以处理多个预测变量。

多项式回归 (Polynomial Regression)

当预测变量和响应变量之间的关系为非线性时，扩展线性回归的标准方法是用一个多项式函数来替换标准线性模型： $$y_i={\beta }_0+{\beta }_1{x}_i+{ϵ}_i$$ 替换为：
$$y_i={\beta }_0+{\beta }_1{x}_i+{\beta }_2{x}_i^2+{\beta }_3{x}_i^3+\cdots +{\beta }_d{x}_i^d+{ϵ}_i$$
其中 ${ϵ}_i$ 是误差项。
这个方法在前面有提到过，这种方法被称为多项式回归。虽然次数提高了，但他还是一个线性模型，预测变量为$x_i,x_i^2,x_i^3,\dots ,x_i^d$。
通常不常使用大于 3 或 4 的 $d$ 值，这会导致曲线过于灵活。

上图是复刻书里 Wage 数据集中工资与年龄的关系图，使用最小二乘法拟合一个 4 次多项式的结果（蓝色实线）。虚线表示估计的95%置信区间。
虚线的生成如下，假设我们已经计算了在特定年龄值 $x_0$ 处的拟合：
$$\hat{f}(x_0)={\hat{\beta }}_0+{\hat{\beta }}_1{x}_0+{\hat{\beta }}_2{x}_0^2+{\hat{\beta }}_3{x}_0^3+{\hat{\beta }}_4{x}_0^4)$$

接下来计算这个拟合值的方差$Var(\hat{f}(x_0))$，对每个参考点都如此计算就得到了。

图生成代码如下：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from ISLP
import load_data

Wage = load_data('Wage')
Wage['wage_k'] = Wage['wage'] / 1000
formula = 'wage_k ~ age + I(age**2) + I(age**3) + I(age**4)'
fit = smf.ols(formula, data = Wage).fit()
age_grid = np.linspace(Wage['age'].min(), Wage['age'].max(), 63)
X_test = pd.DataFrame(
{
    'age': age_grid
})
pred = fit.get_prediction(X_test)
pred_df = pred.summary_frame(alpha = 0.05)
plt.figure(figsize = (10, 6))
plt.scatter(Wage['age'], Wage['wage_k'], facecolors = 'none',
    edgecolors = 'gray', alpha = 0.5, label = 'Original Data')
plt.plot(age_grid, pred_df['mean'], color = 'blue', linewidth = 3,
    label = '4th Degree Polynomial Fit')
plt.plot(age_grid, pred_df['mean_ci_lower'], color = 'blue',
    linestyle = '--', linewidth = 1.5, label =
    '95% Confidence Interval')
plt.plot(age_grid, pred_df['mean_ci_upper'], color = 'blue',
    linestyle = '--', linewidth = 1.5)
plt.fill_between(age_grid, pred_df['mean_ci_lower'], pred_df[
    'mean_ci_upper'], color = 'lightblue', alpha = 0.2)
plt.title(
    'Wage (in thousands) vs. Age with 4th Degree Polynomial Fit')
plt.xlabel('Age')
plt.ylabel('Wage (in thousands of $)')
plt.legend()
plt.grid(True, linestyle = '--', alpha = 0.6)
plt.xlim(17, 83)
plt.ylim(0, 200 / 1000 * 1.05)
plt.show()

Wage 内的工资来自两个不同的群体：似乎有一个年收入超过 25 万美元的高收入群体，以及一个低收入群体。我们可以通过将工资分成这两个群体，将其视为一个二元变量，再使用逻辑回归来拟合，即拟合：
$$\mathrm{Pr}(y_i>250|x_i)=\frac{\exp ({\beta }_0+{\beta }_1{x}_i+{\beta }_2{x}_i^2+\cdots +{\beta }_d{x}_i^d)}{1+\exp ({\beta }_0+{\beta }_1{x}_i+{\beta }_2{x}_i^2+\cdots +{\beta }_d{x}_i^d)}$$

上图为拟合结果。图顶部和底部的灰色标记表示高收入者和低收入者的年龄。蓝色实线表示作为年龄函数的成为高收入者的拟合概率。图中也显示了估计的 95% 置信区间。
置信区间很宽，因为高收入者很少，导致系数的方差很大。

阶梯函数 (Step Functions)

本质上来说，就是用分段函数来拟合。将 $X$ 的范围分成多个区间 (bins)，并在每个区间内拟合一个不同的常数。这相当于将一个连续变量转换为一个有序分类变量。

比如我们在 $X$ 的范围内创建切点 $c_1,c_2,\dots ,c_K$，然后构造 $K+1$ 个新变量：
$$\begin{array}{rl}{C}_0(X)& =I(X<c_1),\\ C_1(X)& =I(c_1\le X<c_2),\\ C_2(X)& =I(c_2\le X<c_3),\\ & ⋮\\ C_{K-1}(X)& =I(c_{K-1}\le X<c_K),\\ C_K(X)& =I(c_K\le X),(7.4)\end{array}$$
其中 $I(\cdot )$ 是一个指示函数，当条件为真时返回 1，否则返回 0。例如，$I(c_K\le X)$ 在 $c_K\le X$ 时等于 1，否则等于 0。这些有时被称为虚拟变量 (dummy variables)。
对于任何 $X$ 的值，$C_0(X)+C_1(X)+\cdots +C_K(X)=1$，因为 $X$ 必须恰好位于这 $K+1$ 个区间中的一个。
最后使用最小二乘法，以 $C_1(X),C_2(X),\dots ,C_K(X)$ 作为预测变量来拟合一个线性模型：
$$y_i={\beta }_0+{\beta }_1{C}_1(x_i)+{\beta }_2{C}_2(x_i)+\cdots +{\beta }_K{C}_K(x_i)+{ϵ}_i$$
对于给定的 $X$ 值，在 $C_1,C_2,\dots ,C_K$ 中最多只有一个可以是非零的。
当 $X<c_1$ 时，可以知道${\beta }_0$ 解释为 $X<c_1$ 时 $Y$ 的平均值。
当 $c_j\le X<c_{j+1}$，响应变量为 ${\beta }_0+{\beta }_j$，因此 ${\beta }_j$ 代表了 $X$ 在 $c_j\le X<c_{j+1}$ 区间内相对于 $X<c_1$ 时响应的平均增量。

上图用了阶梯函数来复刻，当然我本人不知道值，我大概的取了四个数来分段。
代码如下：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from ISLP
import load_data
Wage = load_data('Wage')
cutpoints = [35, 50, 60, 65]
Wage['Age_Cut'] = pd.cut(Wage['age'],
    bins = [Wage['age'].min()] + cutpoints + [Wage['age'].max()],
    include_lowest = True,
    right = True)
X = Wage[['Age_Cut']]
y = Wage['wage']
X_dummies = pd.get_dummies(Wage['Age_Cut'], drop_first = True, dtype =
    int)
X_ols = sm.add_constant(X_dummies)
model_ols = sm.OLS(y, X_ols).fit()
age_grid_values = np.linspace(Wage['age'].min(), Wage['age'].max(),
    500)
bins = [Wage['age'].min()] + cutpoints + [Wage['age'].max()]
grid_cut = pd.cut(age_grid_values, bins = bins, include_lowest = True,
    right = True)
grid_dummies = pd.get_dummies(grid_cut, drop_first = True, dtype =
    int)
X_grid_ols = sm.add_constant(grid_dummies)
predictions_ols = model_ols.get_prediction(X_grid_ols)
summary_ols = predictions_ols.summary_frame(alpha = 0.05)
try:
pred_mean = summary_ols['mean']
lower_ci = summary_ols['mean_ci_lower']
upper_ci = summary_ols['mean_ci_upper']
except KeyError:
    pred_mean = summary_ols['predicted']
lower_ci = summary_ols['ci_lower']
upper_ci = summary_ols['ci_upper']
fig, ax = plt.subplots(figsize = (6, 8))
ax.scatter(Wage['age'], Wage['wage'], facecolors = 'none', edgecolors =
    'gray', alpha = 0.6)
ax.plot(age_grid_values, pred_mean, color = 'green', linewidth = 2,
    linestyle = '-', label = 'OLS Step Function Fit')
ax.plot(age_grid_values, lower_ci, color = 'green', linewidth = 2,
    linestyle = '--', label = '95% CI')
ax.plot(age_grid_values, upper_ci, color = 'green', linewidth = 2,
    linestyle = '--')
ax.set_xlabel('Age', fontsize = 12)
ax.set_ylabel('Wage', fontsize = 12)
ax.set_xlim(17, 83)
ax.set_ylim(30, 310)
plt.title('Step Function Fit for Wage vs. Age', fontsize = 10)
plt.show()

此外还有逻辑回归模型，跟上面多项式回归一个道理。
$$\mathrm{Pr}(y_i>250|x_i)=\frac{\exp ({\beta }_0+{\beta }_1{C}_1(x_i)+\cdots +{\beta }_K{C}_K(x_i))}{1+\exp ({\beta }_0+{\beta }_1{C}_1(x_i)+\cdots +{\beta }_K{C}_K(x_i))}$$

代码如下：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from ISLP
import load_data
Wage = load_data('Wage')
Wage['High_Earner'] = (Wage['wage'] > 250).astype(int)
cutpoints = [33, 40, 50, 60, 65]
bins = [Wage['age'].min()] + cutpoints + [Wage['age'].max()]
Wage['Age_Cut'] = pd.cut(Wage['age'],
    bins = bins,
    include_lowest = True,
    right = True)
y = Wage['High_Earner']
X_dummies = pd.get_dummies(Wage['Age_Cut'], drop_first = True, dtype =
    int)
X_logit = sm.add_constant(X_dummies)
model_logit = sm.Logit(y, X_logit).fit(disp = False)
age_grid_values = np.linspace(Wage['age'].min(), Wage['age'].max(),
    500)
grid_cut = pd.cut(age_grid_values, bins = bins, include_lowest = True,
    right = True)
grid_dummies = pd.get_dummies(grid_cut, drop_first = True, dtype =
    int)
X_grid_logit = sm.add_constant(grid_dummies)
predictions_logit = model_logit.get_prediction(X_grid_logit)
summary_logit = predictions_logit.summary_frame(alpha = 0.05)
pred_probs = summary_logit['predicted']
lower_ci = summary_logit['ci_lower']
upper_ci = summary_logit['ci_upper']
fig, ax = plt.subplots(figsize = (6, 8))
ax.plot(age_grid_values, pred_probs, color = 'darkgreen', linewidth =
    2, linestyle = '-', label = 'Logit Step Function Fit')
ax.plot(age_grid_values, lower_ci, color = 'darkgreen', linewidth = 2,
    linestyle = '--', label = '95% CI')
ax.plot(age_grid_values, upper_ci, color = 'darkgreen', linewidth = 2,
    linestyle = '--')
ax.fill_between(age_grid_values, lower_ci, upper_ci, color =
    'darkgreen', alpha = 0.05)
age_low = Wage.loc[Wage['High_Earner'] == 0, 'age']
age_high = Wage.loc[Wage['High_Earner'] == 1, 'age']
ax.plot(age_low, np.zeros_like(age_low), '|', color = 'gray', alpha =
    0.8, markeredgewidth = 1)
ax.plot(age_high, np.full_like(age_high, 0.20), '|', color = 'gray',
    alpha = 0.8, markeredgewidth = 1)
ax.set_xlabel('Age', fontsize = 12)
ax.set_xlim(17, 83)
ax.set_ylim(-0.01, 0.21)
ax.set_yticks(np.arange(0.0, 0.21, 0.05))
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.title(
    'FIGURE 7.2 (Right): Logit Step Function Fit for Pr(Wage > $250k)',
    fontsize = 10)
plt.show()

很明显，分段函数需要有断点 (breakpoints) 才能发挥作用，在这个数据集不行。

基函数 (Basis Functions)

多项式回归和分段常数回归模型实际上是基函数 (basis function) 方法的特例。
其思想是拥有一系列可以应用于变量 $X$ 的函数或变换：$b_1(X),b_2(X),\dots ,b_K(X)$。我们不再对 $X$ 拟合线性模型，而是拟合如下模型：
$$y_i={\beta }_0+{\beta }_1{b}_1(x_i)+{\beta }_2{b}_2(x_i)+{\beta }_3{b}_3(x_i)+\cdots +{\beta }_K{b}_K(x_i)+{ϵ}_i.$$
基函数 $b_1(\cdot ),b_2(\cdot ),\dots ,b_K(\cdot )$ 是固定的和已知的。对于多项式回归，基函数是 $b_j(x_i)=x_i^j$；对于分段常数函数，它们是 $b_j(x_i)=I(c_j\le x_i<c_{j+1})$。
因此，可以使用最小二乘法来估计未知的回归系数。

回归样条 (Regression Splines)

还是基函数。

分段多项式 (Piecewise Polynomials)

回归样条 (regression spline) 不是在 $X$ 的整个范围上拟合一个高次多项式，而是在 $X$ 的不同区域拟合单独的低次多项式，这被称为分段多项式回归 (piecewise polynomial regression)。例如，一个分段三次多项式 (piecewise cubic polynomial) 的拟合方式是，拟合一个形如下式的三次回归模型：
$$y_i={\beta }_0+{\beta }_1{x}_i+{\beta }_2{x}_i^2+{\beta }_3{x}_i^3+{ϵ}_i,$$
其中系数 ${\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3$ 在 $X$ 范围的不同部分是不同的。系数发生变化的这些点被称为节点 (knots)。

例如，一个没有节点的分段三次多项式就是一个标准的三次多项式。一个在点 $c$ 处有单个节点的分段三次多项式具有以下形式：
$$y_i=\{\begin{array}{ll}{\beta }_{01}+{\beta }_{11}{x}_i+{\beta }_{21}{x}_i^2+{\beta }_{31}{x}_i^3+{ϵ}_i& \text{if }{x}_i<c\\ {\beta }_{02}+{\beta }_{12}{x}_i+{\beta }_{22}{x}_i^2+{\beta }_{32}{x}_i^3+{ϵ}_i& \text{if }{x}_i\ge c.\end{array}$$
约等于把两个模型合一块了，一个用于观测值子集 $x_i<c$，另一个用于观测值子集 $x_i\ge c$。两个都可以最小二乘拟合。
同样可以拟合分段线性 (piecewise linear) 函数。

上图是三次拟合的结果，其中在 age=50 处有一个单节点。
可以发现函数是不连续的 (discontinuous)。
（这里书有一个很明显地断层，我拟合的数据里没有，他甚至接上了）

约束和样条 (Constraints and Splines)

分段多项式的效果不好，在这里就多加一个约束，保证拟合的曲线必须是连续的。

上图追加了约束：在 age=50 处不能有跳跃。
不过函数还是有点不自然。

上图为增加了两个额外的约束：现在分段多项式的一阶导数和二阶导数在 age=50 处都是连续的。这是三个约束，因此图中的曲线被称为三次样条 (cubic spline)。
一般来说，一个有 $K$ 个节点的三次样条总共使用 $4+K$ 个自由度。

上图是一个线性样条 (linear spline)，它在 age=50 处是连续的。
$d$ 次样条的一般定义是：它是一个分段的 $d$ 次多项式，并且在每个节点处，其直到 $d-1$ 阶的导数都是连续的。

四张图的代码如下，我觉得大伙更想要子图的形式。

Wage = load_data('Wage')
knot = 50
age_grid_values = np.linspace(Wage['age'].min(), Wage['age'].max(), 500)
Age_Grid = pd.DataFrame({'age': age_grid_values})
y = Wage['wage']

def plot_setup(ax, title, color='blue'):
    ax.scatter(Wage['age'], Wage['wage'], facecolors='none', edgecolors='gray', alpha=0.6)
    ax.set_title(title, fontsize=12)
    ax.set_xlabel('Age', fontsize=10)
    ax.set_ylabel('Wage', fontsize=10)
    ax.set_xlim(17, 78)
    ax.set_ylim(30, 310)
    ax.axvline(x=knot, color='black', linestyle=':', linewidth=1) 

    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    
def fit_and_plot(formula, data, grid_data, ax, title, color):
    X = dmatrix(formula, data=data, return_type="dataframe")
    X_grid = dmatrix(formula, data=grid_data, return_type="dataframe")
    model = sm.OLS(y, X).fit()
    predictions = model.predict(X_grid)
    plot_setup(ax, title, color)
    ax.plot(age_grid_values, predictions, color=color, linewidth=2)

fig, axes = plt.subplots(2, 2, figsize=(10, 10))
fig.suptitle('FIGURE 7.3: Various Piecewise Polynomials and Splines (Knot at Age=50)', fontsize=14, y=1.02)

formula_tl = f'age + I(age**2) + I(age**3) + I(np.maximum(0, age - {knot})) + I(np.maximum(0, age - {knot})**2) + I(np.maximum(0, age - {knot})**3)'
fit_and_plot(formula_tl, Wage, Age_Grid, axes[0, 0], 'Piecewise Cubic', 'blue')

formula_tr = f'age + I(age**2) + I(age**3) + I(np.maximum(0, age - {knot})**2) + I(np.maximum(0, age - {knot})**3)'
fit_and_plot(formula_tr, Wage, Age_Grid, axes[0, 1], 'Continuous Piecewise Cubic', 'green')

formula_bl = f'bs(age, knots=[{knot}], degree=3)'
fit_and_plot(formula_bl, Wage, Age_Grid, axes[1, 0], 'Cubic Spline', 'red')

formula_br = f'bs(age, knots=[{knot}], degree=1)'
fit_and_plot(formula_br, Wage, Age_Grid, axes[1, 1], 'Linear Spline', 'darkred')

plt.tight_layout()
plt.show()

样条基表示 (The Spline Basis Representation)

这一节表示如何“定义”约束。可以使用基模型来表示一个回归样条。
一个有 $K$ 个节点的三次样条可以被建模为： $$y_i={\beta }_0+{\beta }_1{b}_1(x_i)+{\beta }_2{b}_2(x_i)+\cdots +{\beta }_{K+3}{b}_{K+3}(x_i)+{ϵ}_i,$$
只要为基函数 $b_1,b_2,\dots ,b_{K+3}$ 做出适当的选择即可。

截断幂基函数 (truncated power basis function) 的定义如下：
$$h(x,\xi )=(x-\xi {)}_{+}^3=\{\begin{array}{ll}(x-\xi {)}^3& \text{if }x>\xi \\ 0& \text{otherwise},\end{array}(7.10)$$
其中 $\xi$ 是节点。可以证明，在三次多项式的模型中添加一个形如 ${\beta }_{4}h(x,\xi )$ 的项，将只会导致在 $\xi$ 处的三阶导数出现不连续；而函数本身以及它的一阶和二阶导数在每个节点处都将保持连续。

书里展示了结果，在边界区域的置信区间显得相当不稳定。

我自己复刻之后是下面的。

不知道书里怎么做到的，我用三次样条复刻的，看起来置信区间非常稳定。
代码如下：

Wage = load_data('Wage')
X = Wage['age']
y = Wage['wage']
knots = [25, 40, 60]
age_grid = np.linspace(X.min(), X.max(), 500)
def truncated_power_basis(X, knots):
    X = np.asarray(X)
    T = pd.DataFrame({'age': X})
    for k in knots:
        T[f't{k}'] = np.maximum(0, X - k)**3
    return T
X_spline_basis = truncated_power_basis(X, knots)
X_spline = sm.add_constant(X_spline_basis)
X_spline['age2'] = X**2
X_spline['age3'] = X**3
fit_spline = sm.OLS(y, X_spline).fit()
X_grid_basis = truncated_power_basis(age_grid, knots)
X_grid_spline = sm.add_constant(X_grid_basis)
X_grid_spline['age2'] = age_grid**2
X_grid_spline['age3'] = age_grid**3
pred_spline = fit_spline.get_prediction(X_grid_spline)
y_spline_pred = pred_spline.predicted_mean
ci_spline = pred_spline.conf_int(obs=True) 
X_nat_spline_basis = dmatrix("cr(X, knots=knots)", {"X": X, "knots": knots}, return_type="dataframe")
fit_nat_spline = sm.OLS(y, X_nat_spline_basis).fit()
X_grid_nat_spline_basis = dmatrix("cr(X, knots=knots)", {"X": age_grid, "knots": knots}, return_type="dataframe")
pred_nat_spline = fit_nat_spline.get_prediction(X_grid_nat_spline_basis)
y_nat_spline_pred = pred_nat_spline.predicted_mean
ci_nat_spline = pred_nat_spline.conf_int(obs=True) 
plt.figure(figsize=(10, 6))
plt.scatter(X, y, facecolors='none', edgecolors='gray', alpha=0.8)
for k in knots:
    plt.axvline(k, color='gray', linestyle='--')
plt.plot(age_grid, y_spline_pred, color='blue', linestyle='-', linewidth=2, label='Cubic Spline')
plt.plot(age_grid, ci_spline[:, 0], color='blue', linestyle='--')
plt.plot(age_grid, ci_spline[:, 1], color='blue', linestyle='--')
plt.plot(age_grid, y_nat_spline_pred, color='red', linestyle='-', linewidth=2, label='Natural Cubic Spline')
plt.plot(age_grid, ci_nat_spline[:, 0], color='red', linestyle='--')
plt.plot(age_grid, ci_nat_spline[:, 1], color='red', linestyle='--')
plt.xlabel('Age', fontsize=14)
plt.ylabel('Wage', fontsize=14)
plt.yticks(np.arange(50, 300, 50)) 
plt.xlim(18, 77) 
plt.ylim(30, 290) 
plt.legend(fontsize=12)
plt.show()

选择节点的数量和位置

如何科学地选择节点的数量和位置至关重要。

这里书里提到了两种方法：

节点位置

1. 基于先验知识手动放置（不常用）
   • 思路：在你认为函数变化剧烈的区域（例如，斜率变化大的地方）放置更多节点，在函数平稳的区域放置较少节点。
   • 缺点：这需要你对数据有深入的了解，并且主观性强，不便于自动化。
2. 均匀分布（常用且自动化）
   • 思路：这是最常用的方法。你只需要指定想要的自由度，软件会自动将节点放置在预测变量的均匀分位数上。

上图就是均匀分布，白色虚线是选择的位置。节点的位置是自动选择的，即年龄的第25、50和75百分位数。
代码如下：

Wage = load_data('Wage')
X = Wage['age']
y = Wage['wage']
y_binary = (Wage['wage'] > 250).astype(int)
age_grid = np.linspace(X.min(), X.max(), 500)
df_val = 4 
X_nat_spline_basis_ols = dmatrix(f"cr(X, df={df_val})", {"X": X}, return_type="dataframe")
fit_nat_spline_ols = sm.OLS(y, X_nat_spline_basis_ols).fit()
X_grid_nat_spline_basis_ols = dmatrix(f"cr(X, df={df_val})", {"X": age_grid}, return_type="dataframe")
pred_nat_spline_ols = fit_nat_spline_ols.get_prediction(X_grid_nat_spline_basis_ols)
y_nat_spline_pred_ols = pred_nat_spline_ols.predicted_mean
ci_nat_spline_ols = pred_nat_spline_ols.conf_int(obs=True) 
knots_for_plot = np.percentile(X, [25, 50, 75])
X_nat_spline_basis_log = dmatrix(f"cr(X, df={df_val})", {"X": X}, return_type="dataframe")
fit_nat_spline_log = sm.GLM(y_binary, X_nat_spline_basis_log, family=sm.families.Binomial()).fit()
X_grid_nat_spline_basis_log = dmatrix(f"cr(X, df={df_val})", {"X": age_grid}, return_type="dataframe")
pred_nat_spline_log = fit_nat_spline_log.get_prediction(X_grid_nat_spline_basis_log)
y_nat_spline_pred_log = pred_nat_spline_log.predicted_mean
ci_nat_spline_log = pred_nat_spline_log.conf_int(alpha=0.05)
fig, axes = plt.subplots(1, 2, figsize=(16, 6)) 
fig.suptitle('Natural Cubic Spline', fontsize=16, y=1.02) 
ax0 = axes[0]
ax0.scatter(X, y, facecolors='none', edgecolors='gray', alpha=0.6, s=20) 
for k in knots_for_plot:
    ax0.axvline(k, color='gray', linestyle='--', linewidth=0.8)
ax0.plot(age_grid, y_nat_spline_pred_ols, color='red', linestyle='-', linewidth=2)
ax0.plot(age_grid, ci_nat_spline_ols[:, 0], color='red', linestyle='--', linewidth=1)
ax0.plot(age_grid, ci_nat_spline_ols[:, 1], color='red', linestyle='--', linewidth=1)
ax0.set_xlabel('Age', fontsize=14)
ax0.set_ylabel('Wage', fontsize=14)
ax0.set_yticks(np.arange(50, 350, 50)) 
ax0.set_xlim(18, 82)
ax0.set_ylim(30, 320)
ax1 = axes[1]
ax1.plot(X[y_binary == 0], [0.005] * sum(y_binary == 0), '|', color='gray', alpha=0.5)
ax1.plot(X[y_binary == 1], [0.195] * sum(y_binary == 1), '|', color='gray', alpha=0.5)
for k in knots_for_plot:
    ax1.axvline(k, color='gray', linestyle='--', linewidth=0.8)
ax1.plot(age_grid, y_nat_spline_pred_log, color='red', linestyle='-', linewidth=2)
ax1.plot(age_grid, ci_nat_spline_log[:, 0], color='red', linestyle='--', linewidth=1)
ax1.plot(age_grid, ci_nat_spline_log[:, 1], color='red', linestyle='--', linewidth=1)
ax1.set_xlabel('Age', fontsize=14)
ax1.set_ylabel('Pr(Wage>250 | Age)', fontsize=14)
ax1.set_xlim(18, 82)
ax1.set_ylim(-0.01, 0.21) 
ax1.set_yticks(np.arange(0.00, 0.21, 0.05))
plt.tight_layout(rect=[0, 0, 1, 0.98]) 
plt.show()

节点数量

1. 主观判断：看“样子”
   • 做法：尝试不同自由度的样条，画出拟合曲线，选择一条看起来最平滑、又能捕捉到主要趋势的曲线。
   • 缺点：非常主观，不同的人可能会做出不同的选择。
2. 客观方法：交叉验证（推荐）
   • 目标：找到一个在未见过的数据上表现最好的模型，从而在“欠拟合”和“过拟合”之间取得最佳平衡。
   • 十折交叉验证的具体步骤：
     1. 将数据集随机分成10个大小相似的子集（或“折”）。
     2. 对于每一个候选的自由度 KK：
        • 将第1折数据作为测试集，用剩下的9折数据（训练集）拟合一个具有 KK 个自由度的样条模型。
        • 用这个拟合好的模型对第1折测试集进行预测，并计算预测值与真实值之间的误差（如残差平方和 RSS）。
        • 重复此过程，让每一折数据都轮流充当一次测试集。
     3. 将10次计算出的误差（如RSS）求平均，得到这个自由度 KK 下的交叉验证均方误差。
   • 如何选择：对一系列不同的 KK 值重复上述过程。最终，我们选择那个能产生最小交叉验证均方误差的 KK 值。

书里的图如上。
可以看见3个自由度之后，误差曲线就平摊了。选择3个自由度是一个经济且有效的选择。
左图为自然三次样条，右图为三次样条，都差不多。

但目前python复刻这图，最低要求3个自由度，df=1和df=2的情况对不上。

代码如下：

Wage = load_data('Wage')
dfs = np.arange(3, 11)  
cv = KFold(n_splits=10, shuffle=True, random_state=42)
cv_mse_natural = [] 
cv_mse_cubic = []   
def cv_spline_selection(df_range, data, cv_splitter, spline_type):
    mse_values = []
    for df in df_range:
        fold_mses = []   
        for train_index, test_index in cv_splitter.split(data):
            train_data = data.iloc[train_index]
            test_data = data.iloc[test_index]      
            if spline_type == 'ns':      
                formula = f"wage ~ bs(age, df={df}, degree=3, include_intercept=False)"
            elif spline_type == 'cs':      
                formula = f"wage ~ cr(age, df={df})"
            else:
                raise ValueError("spline_type 必须是 'ns' 或 'cs'")   
            model = ols(formula, data=train_data).fit()  
            predictions = model.predict(test_data)
            rss = np.sum((test_data['wage'] - predictions) ** 2)
            mse = rss / len(test_data)
            fold_mses.append(mse)
        mse_values.append(np.mean(fold_mses))   
    return mse_values

cv_mse_natural = cv_spline_selection(dfs, Wage, cv, spline_type='ns')
cv_mse_cubic = cv_spline_selection(dfs, Wage, cv, spline_type='cs')
best_df_natural = dfs[np.argmin(cv_mse_natural)]
best_df_cubic = dfs[np.argmin(cv_mse_cubic)]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].plot(dfs, cv_mse_natural, color='red', marker='o', markersize=4, 
             linestyle='-')
axes[0].set_title('Degrees of Freedom of Natural Spline', fontsize=12)
axes[0].set_xlabel('Degrees of Freedom of Natural Spline', fontsize=11)
axes[0].set_ylabel('Mean Squared Error', fontsize=11)
axes[0].set_xticks(dfs[::2])
axes[1].plot(dfs, cv_mse_cubic, color='blue', marker='o', markersize=4, 
             linestyle='-')
axes[1].set_title('Degrees of Freedom of Cubic Spline', fontsize=12)
axes[1].set_xlabel('Degrees of Freedom of Cubic Spline', fontsize=11)
axes[1].set_ylabel('Mean Squared Error', fontsize=11)
axes[1].set_xticks(dfs[::2])
y_min = 1580
y_max = 1685
axes[0].set_ylim(y_min, y_max)
axes[1].set_ylim(y_min, y_max)
y_ticks = np.arange(1600, 1681, 20)
axes[0].set_yticks(y_ticks)
axes[1].set_yticks(y_ticks)
plt.tight_layout(pad=3.0) 
plt.show()

与多项式回归的比较

上图是将一个15个自由度的自然三次样条与一个15阶的多项式进行比较。多项式会在尾部表现出剧烈的抖动行为。

回归样条通常比多项式回归能提供更好的结果。

平滑样条 (Smoothing Splines)

概述

拟合的时候我们希望 RSS $={\sum }_{i=1}^n(y_i-g(x_i){)}^2$ 足够小。然而，这种方法存在一个问题。如果我们不对 $g(x_i)$ 施加任何约束，我们总可以通过选择一个能对所有 $y_i$ 进行插值 (interpolates) 的函数 $g$ 来使 RSS 等于零。

我们真正想要的是一个能使 RSS 很小，并且自身也很平滑 (smooth) 的函数 $g$。即最小化：
$$\sum _{i=1}^n(y_i-g(x_i){)}^2+\lambda \int g^{\prime \prime }(t{)}^{2}dt$$
其中 $\lambda$ 是一个非负的调节参数 (tuning parameter)。最小化函数 $g$ 被称为平滑样条 (smoothing spline)。

这个函数包含三个东西：

1. 损失函数 (Loss Function) ：${\sum }_{i=1}^n(y_i-g(x_i){)}^2$
   • 就是 RSS
2. 惩罚项 (Penalty Term) $\lambda \int g^{\prime \prime }(t{)}^{2}dt$
   • 目标：惩罚函数 $g(x)$ 的“不光滑”或“粗糙度”。
   • $\int [g^{\prime \prime }(t){]}^{2}dt$ 的值越大，说明 $g(x)$ 约粗糙，反之越平滑。
3. 调节参数 (Tuning Parameter)$\lambda$
   • $\lambda$ 控制了惩罚项的权重，从而在“拟合数据”和“保持平滑”之间进行权衡。
   • $\lambda \to 0$：惩罚项失效。$g(x)$编程差值曲线。
   • $\lambda \to \infty$：惩罚项占绝对主导。$g(x)$ 退化为普通的现行回归线
   • 适中的 $\lambda$：合理的平滑曲线。

可以证明，$g(x)$ 具有一些特殊的性质：它是一个在 $x_1,...,x_n$ 的唯一值上设有节点的分段三次多项式 (piecewise cubic polynomial)，并且在每个节点处都具有连续的一阶和二阶导数。
换句话说，函数 $g(x)$ 是一个以 $x_1,...,x_n$ 为节点的自然三次样条。

选择平滑参数 $\lambda$

调节参数 $\lambda$ 控制着平滑样条的粗糙度，并因此控制着其有效自由度 (effective degrees of freedom)。可以证明，当 $\lambda$ 从 0 增加到 $\infty$ 时，我们记为 $d{f}_{\lambda }$ 的有效自由度，会从 $n$ 减少到 2。

$${\hat{\mathbf{g}}}_{\lambda }={\mathbf{S}}_{\lambda }\mathbf{y}$$
其中 ：

• ${\hat{\mathbf{g}}}_{\lambda }$ 是训练点 $x_1,...,x_n$ 上的拟合值的 $n$ 维向量。
• ${\mathbf{S}}_{\lambda }$ 是一个 的平滑矩阵。
• $\mathbf{y}$ 是相应向量。
  那么，有效自由度就可以被定义为：
  d f λ = ∑ i = 1 n { S λ } i i df_{\lambda} = \sum_{i=1}^{n} \{\mathbf{S}_{\lambda}\}_{ii} dfλ​=i=1∑n​{Sλ​}ii​ 即矩阵 ${\mathbf{S}}_{\lambda }$ 的迹，对角线元素之和。

对于这类超参数的选择，交叉验证永远是个好办法。
对于平滑样条，可以使用留一法交叉验证 (Leave-One-Out Cross-Validation, LOOCV) 的误差可以非常高效地计算出来，其计算成本与单次拟合几乎相同，使用的公式如下：
RSS c v ( λ ) = ∑ i = 1 n ( y i − g ^ λ ( − i ) ( x i ) ) 2 = ∑ i = 1 n [ y i − g ^ λ ( x i ) 1 − { S λ } i i ] 2 \text{RSS}_{cv}(\lambda) = \sum_{i=1}^{n} (y_i - \hat{g}_{\lambda}^{(-i)}(x_i))^2 = \sum_{i=1}^{n} \left[ \frac{y_i - \hat{g}_{\lambda}(x_i)}{1 - \{\mathbf{S}_{\lambda}\}_{ii}} \right]^2 RSScv​(λ)=i=1∑n​(yi​−g^​λ(−i)​(xi​))2=i=1∑n​[1−{Sλ​}ii​yi​−g^​λ​(xi​)​]2
其中：

• ${\hat{g}}_{\lambda }^{(-i)}(x_i)$：是使用全部数据拟合的模型在第 $i$ 个点上的普通残差。
• { S λ } i i \{\mathbf{S}_{\lambda}\}_{ii} {Sλ​}ii​：是平滑矩阵 ${\mathbf{S}}_{\lambda }$ 的第 $i$ 个对角线元素，被称为该点的杠杆值。
• 1 − { S λ } i i 1-\{\mathbf{S}_{\lambda}\}_{ii} 1−{Sλ​}ii​：起到了调整作用。如果一个点的杠杆值很高（即 { S λ } i i \{\mathbf{S}_{\lambda}\}_{ii} {Sλ​}ii​ 接近1），说明模型在很大程度上“依赖”这个点来塑造自身。当这个点被剔除时，拟合结果会变化很大。因此，它的普通残差会被放大，以更好地估计其真正的预测误差。

这样只需要 1 次拟合就可以计算出 LOOCV 误差。

上图是平滑样条，红色是预先指定 df=16，蓝色是LOOCV自动选择的 $\lambda$ 值得到的平滑样条。
当然书里是 6.8，我这里计算的是 7，我严重怀疑 Wage 数据集被改了。

df 越小模型越简单，选择蓝色曲线自然更好，更平滑。

代码如下：

Wage = load_data('Wage')
age_grid = pd.DataFrame({'age': np.linspace(Wage['age'].min(), Wage['age'].max(), 500)})
DF_HIGH = 16 
ns_high_formula = f"wage ~ bs(age, df={DF_HIGH}, degree=3, include_intercept=False)"
ns_high_model = ols(ns_high_formula, data=Wage).fit()
pred_high = ns_high_model.predict(age_grid)
DF_LOW = 7
ns_low_formula = f"wage ~ bs(age, df={DF_LOW}, degree=3, include_intercept=False)"
ns_low_model = ols(ns_low_formula, data=Wage).fit()
pred_low = ns_low_model.predict(age_grid)
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(Wage['age'], Wage['wage'], facecolors='none', edgecolors='gray', alpha=0.6, s=15, linewidths=0.5)
ax.plot(age_grid['age'], pred_high, color='red', linewidth=2, 
        label=f'{DF_HIGH} Degrees of Freedom')
ax.plot(age_grid['age'], pred_low, color='blue', linewidth=2, 
        label=f'{DF_LOW} Degrees of Freedom')
ax.set_xlabel('Age', fontsize=14)
ax.set_ylabel('Wage', fontsize=14)
ax.set_title('') 
ax.set_xlim(17, 83)
ax.set_ylim(0, 310)
ax.set_xticks(np.arange(20, 90, 10))
ax.set_yticks(np.arange(0, 350, 50))
ax.legend(loc='upper right', fontsize=10)
plt.tight_layout()
plt.show()

局部回归 (Local Regression)

局部回归为每一个预测点 $x_0$ 都单独拟合一个模型，而且这个模型只依赖于 $x_0$ 附近的数据点

蓝色曲线代表生成数据的真实函数 $f(x)$，而浅橙色曲线对应于局部回归的估计值 $\hat{f}(x)$。

在 $X=x_0$ 处的局部回归算法

1. 收集训练点中，其 $x_i$ 值最接近 $x_0$ 的那一部分数据（比例为 $s=k/n$）。
2. 为这个邻域内的每个点分配一个权重 $K_{i0}=K(x_i,x_0)$，使得离 $x_0$ 最远的点权重为零，而最近的点权重最高。除了这 $k$ 个最近邻点外，所有其他点的权重都为零。
3. 使用上述权重，通过找到最小化以下公式的 ${\hat{\beta }}_0$ 和 ${\hat{\beta }}_1$ 来对 $y_i$ 和 $x_i$ 进行加权最小二乘回归 (weighted least squares regression)：
   $$\sum _{i=1}^n{K}_{i0}(y_i-{\beta }_0-{\beta }_1{x}_i{)}^2$$
4. 在 $x_0$ 处的拟合值为 $\hat{f}(x_0)={\hat{\beta }}_0+{\hat{\beta }}_1{x}_0$。

局部回归有时被称为 “基于内存”的程序 (memory-based procedure)，这是因为每一次都要重复算法实时计算。
图代码如下：

np.random.seed(42)
N = 100
X = np.sort(np.random.uniform(0, 1, N))
f_true = 1.0 * np.sin(X * np.pi * 1.5) + 0.5 * X 
Y = f_true + np.random.normal(0, 0.25, N) 
X_grid = np.linspace(0, 1, 500)
f_true_grid = 1.0 * np.sin(X_grid * np.pi * 1.5) + 0.5 * X_grid
loess_fit = lowess(Y, X, frac=0.4, is_sorted=True)
X_loess, Y_loess = loess_fit[:, 0], loess_fit[:, 1]

def plot_local_fit(ax, X, Y, x0, frac, color_main='tab:red', color_weight='gold'):
    bandwidth = np.abs(X - x0).max() * frac
    dist = np.abs(X - x0)
    weights_raw = np.maximum(0, 1 - (dist / bandwidth)**3)**3
    
    idx_local = weights_raw > 0.001
    X_local, Y_local, W_local = X[idx_local], Y[idx_local], weights_raw[idx_local]

    X_mat = np.vstack([np.ones_like(X_local), X_local]).T
    W_diag = np.diag(W_local)
    try:
        beta = np.linalg.inv(X_mat.T @ W_diag @ X_mat) @ X_mat.T @ W_diag @ Y_local
    except np.linalg.LinAlgError:
        return
    x_range = np.array([X_local.min(), X_local.max()])
    y_fit_line = beta[0] + beta[1] * x_range
    ax.plot(x_range, y_fit_line, color=color_main, linewidth=2)
    ax.plot([x0, x0], [-1.0, 1.5], color=color_main, linewidth=1, linestyle='-')
    ax.scatter(X_local, Y_local, facecolors='none', edgecolors=color_main, alpha=0.8, s=40, linewidths=1.5)
    f_hat_x0 = beta[0] + beta[1] * x0
    ax.plot(x0, f_hat_x0, 'o', color=color_main, markersize=8)
    mu = x0
    sigma = bandwidth / 2.5 
    x_curve = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)
    x_curve = x_curve[(x_curve >= 0) & (x_curve <= 1)]
    y_curve = norm.pdf(x_curve, loc=mu, scale=sigma)
    y_curve_scaled = (y_curve / y_curve.max()) * 1.0 
    ax.fill_between(x_curve, 0.0, y_curve_scaled, color=color_weight, alpha=0.7)
    return f_hat_x0
x0_left = 0.05
x0_right = 0.45
fig, axes = plt.subplots(1, 2, figsize=(12, 5.5), sharey=True)
ax = axes[0]
ax.set_title(f'Target Point $x_0 = {x0_left}$ (Boundary)', fontsize=12)
ax.plot(X_grid, f_true_grid, color='tab:blue', linewidth=2, label='True Function $f(x)$')
ax.plot(X_loess, Y_loess, color='orange', alpha=0.8, linewidth=2, label='Local Regression $\\hat{f}(x)$')
ax.scatter(X, Y, facecolors='none', edgecolors='gray', alpha=0.6, s=30, linewidths=0.5)
plot_local_fit(ax, X, Y, x0=x0_left, frac=0.4, color_main='orangered', color_weight='gold')
ax = axes[1]
ax.set_title(f'Target Point $x_0 = {x0_right}$ (Interior)', fontsize=12)
ax.plot(X_grid, f_true_grid, color='tab:blue', linewidth=2)
ax.plot(X_loess, Y_loess, color='orange', alpha=0.8, linewidth=2)
ax.scatter(X, Y, facecolors='none', edgecolors='gray', alpha=0.6, s=30, linewidths=0.5)
plot_local_fit(ax, X, Y, x0=x0_right, frac=0.4, color_main='orangered', color_weight='gold')
for ax in axes:
    ax.set_xlim(-0.05, 1.05)
    ax.set_ylim(-1.2, 1.6)
    ax.set_xlabel('X', fontsize=12)
    ax.set_ylabel('Y', fontsize=12)
plt.tight_layout(pad=1.5)
plt.show()

跨度 (span)

跨度定义了在计算 $x_0$ 处的局部回归时所使用的点的比例，是局部回归中最重要的调节参数，其作用类似于平滑样条中的 $\lambda$。

• 小跨度：
  • 使用的邻域很小，只包含很少的数据点。
  • 结果：拟合的曲线可以捕捉非常局部的、细微的变化，但也会变得非常**“曲折”和“嘈杂”** ，方差高，容易过拟合。
• 大跨度：
  • 使用的邻域很大，包含更多的数据点。
  • 结果：拟合的曲线更加平滑，捕捉的是更宏观的趋势，但可能会忽略一些重要的局部特征，偏差高，可能欠拟合。
• 选择跨度：与之前一样，最佳的跨度 $s$ 可以通过交叉验证 来确定。

推广

局部回归的思想非常灵活，可以推广到更复杂的情况：

1. 变系数模型：在一个多变量场景中，可以让模型对某些变量是全局的，而对另一个关键变量（如时间）是局部的。例如，研究收入和年龄的关系，但这个关系可能随时间变化。我们可以让系数随时间进行局部变化，从而捕捉关系的演变。
2. 多维局部回归：局部回归可以自然地扩展到两个或更多个预测变量。例如，要预测点 $(X_1,X_2)$ 的值，可以找到一个二维邻域，然后在这个邻域内拟合一个加权的多元线性回归模型。
3. 局限性：维数灾难
   尽管局部回归在低维情况下表现优异，但当预测变量 $p$ 很大时（例如超过3或4个），它会遭遇与 k-近邻法 同样的维数灾难 问题。
   • 原因：在高维空间中，数据点变得极其稀疏。为了在目标点 $x_0$ 周围找到一个足够大的邻域以包含有代表性的数据点，这个邻域的物理范围可能必须非常大，以至于它不再是“局部”的。这就失去了局部拟合的意义，模型会退化为一个粗糙的全局模型。

广义可加模型 (Generalized Additive Models)

上面都是基于单个预测变量 $X$ 灵活预测响应变量 $Y$ 的方法。
这一节介绍基于多个预测变量 $X_1,...,X_p$ 灵活预测 $Y$ 的问题。

复刻无果，怎么也对不上，放弃了。

用于回归问题的GAMs

一种扩展多元线性回归模型 $$y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_p{x}_{ip}+{ϵ}_i$$ 以允许每个特征与响应变量之间存在非线性关系的自然方法是，将每个线性部分 ${\beta }_j{x}_{ij}$ 替换为一个（平滑的）非线性函数 $f_j(x_{ij})$。然后我们可以将模型写为
$$\begin{array}{rl}{y}_i& ={\beta }_0+\sum _{j=1}^p{f}_j(x_{ij})+{ϵ}_i\\ & ={\beta }_0+f_1(x_{i1})+f_2(x_{i2})+\cdots +f_p(x_{ip})+{ϵ}_i\end{array}$$
它被称为可加模型 (additive model)，因为我们为每个 $X_j$ 计算一个独立的 $f_j$，然后将它们各自的贡献相加。

GAM可以将我们之前学过的所有单变量平滑方法作为“构建模块”来使用。比如说本章常用例子：
考虑在 Wage 数据上拟合以下模型的任务：
$$\text{wage}={\beta }_0+f_1(\text{year})+f_2(\text{age})+f_3(\text{education})+ϵ$$
year 和 age 是定量变量，而变量 education 是一个有五个水平的定性变量：<HS, HS, <Coll, Coll, >Coll，分别指代个人完成的高中或大学教育水平。

上图就是这么拟合的。

还可以使用平滑样条。需要使用反向拟合 (backfitting)，这个方法每次我们更新一个函数时，我们只需将该变量的拟合方法应用于一个偏残差 (partial residual)。

上图是平滑样条拟合的，$f_1$ 和 $f_2$ 分别是自由度为4和5的平滑样条。

GAMs的优缺点

优点：

1. 自动非线性拟合：无需手动为每个变量尝试不同的变换（如对数、平方根），GAM能自动从数据中学习非线性关系。
2. 预测更准确：通过捕捉非线性，通常能获得比线性模型更好的预测性能。
3. 可解释性强：模型的可加结构允许我们单独可视化并解释每个变量与响应之间的关系，这是GAM最强大的优势之一。
4. 平滑度可控：可以为每个 $f_j$ 指定不同的自由度，以控制其平滑度（灵活性）。

缺点：

1. 主要局限：可加性假设：模型假设变量之间没有交互作用。一个变量的效应不依赖于另一个变量的取值。例如，它可能无法捕捉“高等教育带来的工资溢价在年轻人中更高”这类交互效应。
2. 解决方案：
   • 可以手动添加交互项，如 $X_j\times X_k$。
   • 可以包含二维平滑项 $f_{jk}(X_j,X_k)$ 来直接建模两个变量之间的非线性交互作用（但这会增加模型复杂性）。

GAM在线性模型和完全非参数的“黑箱”模型（如随机森林、梯度提升机）之间提供了一个极佳的折中方案。

• 它比线性模型灵活得多，能捕捉复杂的非线性模式。
• 它比一些复杂的集成方法可解释性强得多，结果可以像线性模型一样被可视化和平滑地解释。

用于分类问题的GAMs

标准的逻辑回归如下：
$$\log \left(\frac{p(X)}{1-p(X)}\right)={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p$$
其中$P(Y=1|X)$ 与 $P(Y=0|X)$ 的对数几率 (log of the odds)。
引入非线性函数 $f_j$ 就可以扩展它。
$$\log \left(\frac{p(X)}{1-p(X)}\right)={\beta }_0+f_1(X_1)+f_2(X_2)+\cdots +f_p(X_p)$$
上式为逻辑回归GAM (logistic regression GAM)。

同样将一个GAM拟合到 Wage 数据，以预测个人年收入超过$250,000的概率。我们拟合的GAM形式如下：
$$\log \left(\frac{p(X)}{1-p(X)}\right)={\beta }_0+{\beta }_1\text{year}+f_2(\text{age})+f_3(\text{education})$$
其中，
$$p(X)=\text{Pr}(\text{wage}>250|\text{year},\text{age},\text{education})$$

上图是拟合结果。这里有个问题：在 <HS（低于高中）这个教育水平上，拟合的曲线有一个非常宽的置信区间。因为数据中没有任何一个受教育程度低于高中的人年收入超过25万美元。

上图是从分析中排除 <HS 这个类别的个体后重新拟合的结果，可以看见教育是阶梯式增长了。

优缺点

• 优点：
  • 自动捕捉非线性：无需手动试探变换，就能发现如 age 与收入概率之间的“倒U型”关系。
  • 可解释性强：可以直观地可视化每个变量如何影响成为“高收入者”的对数几率。
  • 预测性能可能更优：通过更准确地建模变量关系，可能得到更好的分类性能。
• 缺点：
  • 可加性假设：模型假设变量之间没有交互作用。例如，它无法捕捉“高等教育带来的收入溢价在中年阶段最为显著”这类交互效应。
  • 解决方案：同样可以手动添加交互项或使用二维平滑项。