【视频】R语言广义加性模型GAMs非线性效应、比较分析草种耐寒性实验数据可视化

全文链接：https://tecdat.cn/?p=36979

原文出处：拓端数据部落公众号

广义加法模型（Generalized Additive Models, GAMs）作为一种高度灵活的统计工具，显著扩展了广义线性模型（Generalized Linear Models, GLMs）的框架。GAMs的核心思想在于，将GLM中的一个或多个线性预测变量替换为这些变量的平滑函数，从而允许模型捕捉预测变量与条件响应之间复杂且非线性的关系，而无需事先对这些关系的具体形态做出假设。这一过程通过引入惩罚平滑样条技术实现，该方法在保持模型灵活性的同时，有效防止了过拟合现象。

GAMs的工作原理根植于基展开（basis expansion）的概念之中。简而言之，基展开意味着将协变量（在此语境下，如时间等）映射到一组精心设计的基函数上，这些基函数旨在全面覆盖协变量观测值的范围。一种常用的基函数类型是三次回归基（cubic regression splines），它们能够灵活地适应数据的局部特征。

值得注意的是，除了三次回归基外，还有多种类型的基展开方法可用于构建惩罚平滑模型，包括但不限于多维平滑技术，用于处理具有多个协变量的复杂情况；空间平滑技术，特别适用于具有空间相关性的数据；以及单调平滑技术，确保模型输出随协变量单调变化，以满足特定领域的分析需求。这些多样化的基展开方法为GAMs提供了广泛的应用前景和强大的建模能力。

何时应使用广义加法模型？

GAM之所以备受推崇，主要归因于其在处理复杂数据关系上的独特优势。具体而言，当研究者面临以下情境时，GAM成为强有力的分析工具：

1. 非线性关系探索：当响应变量与预测变量间可能存在未知的非线性关系时，GAM能灵活捕捉这些复杂模式。
2. 形式无关性：无需预设关系形式，GAM通过数据驱动的方式自动学习最佳拟合模型。
3. 高效捕获非线性效应：相比传统高阶多项式，GAM能以更优雅的方式捕获非线性效应，同时规避其潜在的过拟合和解释复杂性。
4. 稳健性：在追求模型拟合精度的同时，GAM通过惩罚平滑技术有效控制过拟合风险。

环境设置和初始 GAM 模型

现在，加载数据。这些数据包含草种耐寒性实验的测量值

检查数据结构和维度

dplyr::glimpse(CO2)

现在对建模进行一些操作

plant <- CO2 |>
  as_tibble() |>
  rename(plant = Plant, type = Type, treatment = Treatment) |>
  mutate(plant = factor(plant, ordered = FALSE))

我为这些数据拟合了一个广义加性模型，指定了每个截距之间的参数交互作用以及随机（分层）截距，并使用二氧化碳吸收 作为非负的实值响应。非线性部分允许二氧化碳浓度的非线性效应随不同水平的冷处理变量而变化。

看看这个模型的总结

这里似乎有很多“重大”影响，但我们到底如何解释这些呢？

标记的系数是控制这些样条形状的基函数权重。但是，如果不能完全理解这些基本函数是什么样子的，以及它们如何共同作用形成样条曲线。

现在，我将介绍 3 个步骤，您可以使用这些步骤来帮助解释和报告来自 GAMS 的非线性效应。

第 1 步：可视化 GAM 输出

到目前为止，视觉效果是我们理解 GAM 的最佳首选工具。

在链路（或链接）尺度上绘制GAM（广义加性模型）的效应时，用户最常采用的可视化手段是部分效应图（Partial Effects Plots）。这些图形在链路尺度上直观地展示了各个平滑函数（Smooth Functions）的单独分量效应，此过程中假设模型中的其他所有项均被固定为零（即，保持其他因素不变，仅关注某一特定预测变量的影响）。由于多数平滑项在作图时会被中心化为零以提高可解释性，因此这些图形易于解读，使得用户能够迅速获得关于变量间关系的初步理解。

具体而言，若要在GAM中查看特定平滑项（如处理因素“nonchilled treatment”）的部分效应，用户可以通过选择该平滑项并观察其在链路尺度上的表现来实现。

该图显示，对于 的较小值，对线性预测变量的影响大多为负值（低于零），但对于 的中间值，它很快就会变为正值（高于零）。然后，它开始在较大的值处趋于稳定。我们还可以看看其他层次的影响

这个图看起来大致相似，但效果并不完全相同。

在探讨广义加性模型（GAMs）的效应时，确实需要注意其可视化表示的局限性和如何更全面地理解模型结果。

GAM效应的可视化局限性

尽管在链路尺度上绘制GAM的部分效应图是用户常用的可视化手段，但这种方法有其内在的局限性。这些图主要展示了在保持其他所有预测变量为零（或基准水平）的情况下，单个平滑函数对响应变量的预期影响。然而，这种“孤立”的展示方式可能无法全面反映预测变量之间的交互作用以及它们对响应变量的综合影响。

复杂模型结构的挑战

当GAM模型包含多个平滑项和交互项时（如y ~ s(rain) + s(temp) + s(year) + ti(rain, temp) + ti(rain, year) + ti(temp, year)），单一预测变量的效应往往分散在多个平滑函数中，这使得直接解释每个平滑项变得困难。此外，如果模型采用对数链接函数等非线性链接，部分效应图在链路尺度上的展示可能与实际尺度上的影响存在显著差异。

超越链路尺度的分析

为了更全面地理解GAM的效应，我们需要超越简单的链路尺度部分效应图。以下是一些建议的方法：

1. 计算并绘制平均平滑效果：利用适当的统计软件包（如R中的mgcv和ggeffects或margins包），可以计算并绘制考虑所有其他预测变量影响的平均平滑效果图。这样的图能够更好地反映预测变量在实际情境下的综合影响。

2. 转换到实际尺度：如果模型使用了非线性链接函数，应尝试将链路尺度上的效应转换为实际尺度（如原始数据尺度或概率尺度），以便更直观地解释模型结果。

3. 比较不同条件下的效应：通过计算和比较不同治疗组或不同协变量水平下的效应，可以更深入地了解预测变量如何影响响应变量，以及这些影响在不同条件下如何变化。

4. 使用更高级的绘图和摘要工具：采用专门的统计绘图和摘要工具（如ggeffects、sjPlot等R包），可以方便地生成各种类型的效应图，包括条件效应图、交互效应图等，从而更全面地展示GAM的复杂结构。

当然，我们可以很容易地将这些图拆分为：

在探索广义加性模型（GAMs）时，了解平滑函数斜率的变化对于深入理解模型行为及其在不同协变量水平下的影响至关重要。

首先，需要注意的是，plot_predictions()函数通常不直接支持绘制斜率。但是，您可以使用与这些包相关或独立的函数来计算平滑函数的一阶导数，并使用图形化工具（如ggplot2）来展示这些斜率。

该图更清楚地表明，在我们达到 260 附近的值之前，斜率是正的，超过该值，函数将趋于平稳。

如何在结果量表上绘制平滑效应？

在响应量表上绘制 GAM 效应

我们现在可以选择在这些图上显示观测值

plot_predictions(model_1, co

这是再次的平均平滑度，但这次是在响应变量的尺度上。

这些绘图增强了我们对拟合模型进行质疑和评估的能力。在解读或报告GAM中的函数时，您可以考虑以下几个基本问题来启动分析：

• 该函数是否在其定义域内达到渐近线？
• 函数图像是剧烈波动还是展现出平滑的趋势？
• 函数是否存在多个峰值或模式？这些模式在实际应用中是否有合理解释？
• 是否存在数据点稀疏的区域，且该区域函数的不确定性相应增加？
• 是否有明显的异常点，导致函数反应异常强烈？

此外，如果您有兴趣进一步探索，确实可以在响应尺度上直接绘制斜率图，以观察函数变化率随协变量的变化。

到目前为止，您可能对这些预测的来源感到好奇，这正是引入我们解释GAM的下一个有力工具——模拟的恰当时机。

第2步：从拟合的GAM模型进行仿真

在深入探讨GAM时，通过模拟数据来加深对其模型及其潜在局限性的理解变得尤为重要。特别是对于GAM，模拟过程涉及到线性预测器（或称设计矩阵）的生成，这是通过将协变量映射到其对应的基函数上而得到的。

理解GAM模型中系数的含义

一个关键步骤是查看线性预测矩阵（我将其简称为(X_{lp})），这个矩阵对于理解GAM中的系数至关重要。在R中，使用mgcv包中的predict.gam()函数，并设置type = 'lpmatrix'，我们可以轻松地生成这个矩阵。无论是针对新数据还是拟合模型时使用的原始数据，这一操作都同样适用。以下是如何针对原始数据生成线性预测矩阵的示例：

Xp <- predict(model_1, type = 'lpmatrix')
dim(Xp)

## [1] 84 28

您可以看到矩阵中有 84 行，对应于我们数据中的 84 个观测值。但是我们有 28 列，其中许多列表示模型中两个平滑项的基函数

这些对应于我们之前从拟合模型中提取的系数 

## [1] TRUE

如果我们使用线性代数将这些系数与设计矩阵 \（（X_{lp}\beta）\） 交叉相乘，我们会得到链接尺度上的预测值：

通过反向链接函数（在我们的对数链接的情况下）运行这些函数，为我们提供了模型中的拟合值exp()

## [1] TRUE

从模型的隐含多元正态后验分布中抽取(\beta)系数来研究预测中的不确定性是一个高级话题，它涉及到贝叶斯统计和MCMC（马尔可夫链蒙特卡洛）方法，这通常用于更复杂的模型评估。不过，对于大多数GAM（广义加性模型）的常规应用，我们通常关注于点预测和预测区间，这些可以通过predict.gam()函数直接获得，而无需显式地抽取(\beta)系数的后验样本。

现在，让我们聚焦于您提到的实际应用场景：当您向GAM模型提供新数据时，如何利用这些数据进行预测。假设您已经有一个拟合好的GAM模型，该模型研究了不同CO₂浓度和温度处理下植物的生长情况。现在，您想要预测在特定CO₂浓度（如278 ppm）和特定处理（如未冷却处理）下植物的某种响应（比如生长率）。

## [1]  1 28

在这里，我们采用了相对简单的场景，并通过所有必要的基础评估来生成设计矩阵。然后，从此方案生成预期值就像以下简单方法一样简单：

exp(newXp %*% beta)

##       [,1]
## 1 27.56059

可能很难理解这种方法的力量，但一旦你掌握了这些步骤，探索目标场景的可能性是无穷无尽的。例如，我们可以使用这些步骤来查看效果的基础基函数是什么样子的。只需模拟一系列假值以覆盖协变量的范围，并通过函数运行这些值（将其他协变量固定到特定值）：

 
newXp <- predict(model_1, type = 'lpmatrix', newdata = newdat)

求哪些系数属于conc

## [1] 17 18 19 20 21 22

现在将 \（X_{lp}\） 矩阵中与这些系数不对应的所有单元格设置为零

在链路尺度上生成预测并绘制函数

 
ggplot(plot_dat, aes(x = conc, y = value, col = basis_func)) +

确实，当处理包含多个协变量的模型时，手动为所有感兴趣的预测场景创建newdata数据框可能会变得既繁琐又容易出错。这就是为什么自动化工具在这种情况下变得极其有价值的原因。通过使用复杂的规则自动设置缺失预测变量的值，可以毫不费力地创建这些方案。

在这里，您会注意到只有目标变量在行与行之间变化，而其余所有变量均保持固定。为了探索并优化这些变量配置，您可利用预设的默认规则（及其自定义选项），这些规则旨在辅助解决潜在的预测挑战。通过巧妙结合marginaleffects的datagrid函数与mgcv的predict.gam()功能（特别是当type = 'lpmatrix'时），您可以高效地利用虚拟数据来构建具有针对性的查询。值得注意的是，在调用plot_predictions()、plot_slopes()或datagrid()等函数时，背后实际上自动调用了必要的逻辑以精确设置设计矩阵，比如conc等变量的处理。

值得一提的是，marginaleffects的强大之处不仅限于GAM，它提供了一个清晰、简洁的框架来探索非线性效应，同时也广泛兼容R中多种模型类（当前已支持超过100种），这一特性极大地促进了模型间的比较与分析。例如，即便是在处理包含复杂多项式交互效应的GLM（尽管这通常不是一个推荐的做法，仅为示例）时，marginaleffects也能游刃有余地助力您将模型拟合至相同数据集，进而深入洞察数据背后的故事。

我们可以从我们之前制作的 GAM 中重新制作浓度效应的图之一

plot_predictions

我们可以生成完全相同的图，但现在有了 GLM。请注意，除了 model 参数之外，调用 to 中的单个字符都不必更改plot_predictions()

plot_predictions(model_2

如何从我们的GAM模型中提炼出更为直接且深刻的问题呢？

第三步：深化GAM的比较分析

多数模型旨在描述而非直接模拟数据背后的因果机制。然而，这恰好赋予了我们利用回归模型探索不同场景间差异、并询问模型何为合理预测的能力。

差异平滑的奥秘：探索各组之间平滑效果的微妙变化

近期关于计算因子-平滑交互中平滑差异的研究，为我们深入理解{}在GAM比较中的强大作用提供了宝贵视角。当研究者希望拟合类似您之前探讨的GAM（含因子-平滑交互），以洞察非线性效应如何随因子水平变化时，问题虽看似简单，但手动处理却颇为棘手。幸运的是，{}提供了一系列专门设计的函数，极大地简化了这一过程。这些函数能够比较不同场景下的预测，无论是计算差异、比率，还是执行用户自定义的函数（如转换对数赔率或弹性），都游刃有余。

具体而言，我们可以通过整合uptake、type、conc和treatment等变量，计算处理间值的平均预期差异，从而得到一种‘全局’效应度量，该度量不受特定type或conc值的影响，是总结模型效果的有力工具。"

在这里，我们可以清晰地观察到，在反应的尺度上，不同治疗之间的平均差异显著且强于某个特定的基准（尽管您在此处未明确提及该基准是什么，可能是指未治疗组或另一种治疗方式）。这种对比不仅提供了估计的效应大小，还附带了相应的p值，以量化这种差异在统计上的显著性。为了获取这些详细的结果，包括效应大小和p值，您可以使用marginaleffects包的comparisons函数，并通过指定type参数来进一步细分或定制比较的类型。

这为我们提供了两个平滑值之间的预期差值。它非常有用，因为它已经考虑了截距的任何变化或模型中可能出现的其他影响。我们可以绘制这些差异：

我们还可以提出诸如非线性斜率增长最快的 conc 值等问题？这种问题在时间序列建模中非常有用，例如，如果我们想知道人口趋势何时增长（或减少）最快。为函数提供自定义函数可以实现这一点：comparisons()

在这里，我们可以看到，对于这两种处理，非线性效应的斜率似乎在conc = 155附近增长最快。

或者我们可以问：在什么 conc 值下，斜率与零有显着差异？通过再次聚合和计算一系列值的一阶导数，我们可以使用默认情况下询问数量是否与零显着不同的函数：

在深入探讨科学报告中如何精准呈现GAM（广义可加模型）的影响时，我们时常会遭遇需要细化假设的情境，而这些假设的设定往往能借助R来轻松实现。

如何在期刊中精准报告GAM的影响？

最终，我将聚焦于解答GAM领域的一个普遍疑问：如何有效地传达这些复杂而精细的分析结果？以下是我提炼的几点关键建议，旨在促进GAM非线性效应在科学交流中的准确传达与深入理解，特别是借助marginaleffects等工具的辅助：

1. 超越p值与统计显著性：应避免过分依赖p值作为效应存在与否的唯一标准。将连续变化的效应简化为二元分类，往往会忽略大量宝贵的信息。相反，应关注效应的实际大小及其对结果预期变化的贡献。

2. 聚焦效应对结果规模的影响：选择GAM与GLM，正是为了捕捉现实世界中非高斯分布的复杂现象，这意味着链接函数常呈现非线性特性。因此，解释效应时，应努力在线性预测变量的尺度上阐明其含义，尽管这可能对部分读者构成挑战。

3. 利用仿真增强理解：通过仿真手段，我们可以深入探究模型在不同情境下的合理预测范围。这不仅有助于自身对模型行为的理解，也需在报告中详尽描述仿真策略及其合理性依据，从而增强结论的可信度。

4. 对比不同模型以评估稳健性：将GAM与其他模型（如多项式回归、线性模型）进行对比分析，是评估结论对函数形式选择敏感性的重要步骤。这种跨模型的比较能够揭示不同假设下结论的稳定性，为科学论断提供更加坚实的支撑。

综上所述，通过避免对p值的过度依赖、关注效应的实际影响、利用仿真深化理解以及进行跨模型比较，我们可以在期刊中更加全面、准确地报告GAM的非线性效应，促进科学知识的有效传播与应用。