ChatGPT-4O能跑线性混合模型，但不建议这么用

在完成我的研究生水平的生物统计学课程的作业过程中，我们从简单的模型和方法开始，逐渐转向更复杂的模型和方法，这给ChatGPT-4o带来了越来越多的挑战。上周，我们发现它目前无法正确运行广义加性模型（参考：[原文链接]）。这周，我测试了它运行线性混合模型的能力。

混合模型：（https://en.wikipedia.org/wiki/Mixed_model）是一类特定的回归模型，结合了固定效应和随机效应的模型。固定效应是我们通常假设当一个预测变量对响应变量有影响时的情况，我们通常想要了解这个预测变量对响应变量的具体影响。我们使用随机效应来建模那些我们知道与响应变量相关的变量的影响，但我们并不特别关心这种影响的具体内容，我们只是想避免不考虑它们可能带来的偏差。我在教学中向学生介绍这些模型，因为他们经常需要使用这些模型来分析他们论文项目中收集的数据。

我常举的一个例子是教育研究，我们正在研究一种新的教学法对学生群体的影响。我们可以简单地比较使用新教学法的学生和未使用的学生的学习情况。然而，不同的学生群体由不同的老师教授，我们知道不同老师的教学效果各不相同。我们不能忽略这一点，尤其是因为同一老师教的不同学生群体之间不能假定相互独立。如果我们不考虑老师的影响，我们就会陷入“伪重复”的陷阱。然而，由于各种原因，我们也不希望具体测量和比较不同老师在教学效果上的差异。如果把老师作为固定因素处理，会严重限制我们检测新教学法效果的能力。相反，我们将老师（或教室）作为随机效应。这样假设老师的效果围绕平均值正态分布，每个老师被视为从这个正态分布中随机抽取的样本，我们只需估计这个正态分布的方差。这样我们就可以在不直接比较老师之间差异的情况下，考虑老师之间的差异，并且大大减少了统计所需的自由度。

另一个常见的混合效应模型的例子是重复测量设计，例如治疗前后的测量。不同动物对治疗的反应可能各异，我们只关心所有动物平均反应，无需逐个比较每只动物。然而，忽略每次测量来自哪只动物，又假设来自同一个动物的每次测量都是相互独立的，这当然是不正确的，会导致伪重复。将动物身份作为随机效应纳入，假设响应的分布是正态的，这样就能将动物间的反应差异考虑进去，并估计这种变异，而不需要比较每只动物，同时减少了自由度的损耗。

我在课堂上向学生描述这些问题，并展示不同程度的线性模型，包括有或没有交互项和主效应的模型，以及包括有随机截距和随机斜率或仅有随机截距的混合模型。我的目标是展示将预测变量作为固定效应或随机效应的选择如何影响结果，并说明如何选择最佳的方法。我们还使用R语言的nlme包中的lme命令来拟合这些模型，并解释这些模型的摘要输出、效应图和诊断图。之后，他们被要求完成一项作业，练习运行和解释类似模型的结果。这就是本周我给ChatGPT-4o布置的任务。请观看下面的视频演示。

加载了阿根廷数据集之后，我首先让ChatGPT-4o列出列名，它毫不费力地完成了任务。我还让它创建了三个新的分类变量，以确保即使它们被编码为数字，也被视为分类变量。经过一些抱怨之后，它还是照做了。然后，我让它进行了一次简单的线性回归，以 L. acuta 鸟的数量为响应变量，以中密度沙子的比例作为唯一预测变量。ChatGPT-4o能够正确地拟合线性回归，并提供与学生在R中得到的相同的回归线和汇总统计信息。唯一的例外是AIC值，在ChatGPT-4o/Python中比R低了2个单位，因为R在对数似然计算中包含了额外的自由度（用于离散参数的自由度），而Python没有。值得注意的是，模型发现响应变量和预测变量之间存在显著的正相关性。

然而，该模型的问题在于数据是在分配给三个不同样带的地块中收集的，而未将样带作为预测变量，我们基本上假设了同一样带内的地块彼此独立，这可能导致伪重复取样。因此，我要求它运行一个线性回归，使用相同的响应变量，但将MedSand和fTransect（样带的分类版本）及其交互作用作为预测变量（因此是一个交互模型）。经过与变量名“L.acuta”中的句点的困扰之后，ChatGPT-4o成功拟合了正确的模型，并提供了所有正确的汇总统计信息，除了AIC的2单位差异。重要的是，MedSand和L.acuta之间的关系不再显著（这表明了纳入fTransect变量的重要性）。ChatGPT-4o在被要求提供响应变量由预测变量解释的变异比例时选择了多重R平方，而不是更合适的调整R平方，这也是我的学生常犯的错误。在创建效应图时，我遇到了真正的挑战。首先，它将三个面板分别绘制出来。然后，当我要求它重叠绘制三个回归线时，它分别对三个样带拟合模型。虽然对最佳拟合线来说是可以接受的，但在被要求绘制95%置信区间时这种方法就不适用了。当我尝试将轴的范围设置成与R中的相同，我发现置信区间看起来不同，这时我回头看到，它实际上是按样带分割了数据。当被要求使用交互模型时，它遵从了要求，我可以看到它正在使用Python中的get_prediction命令来获取每个样带的预测值和95%置信区间。这令人担忧，因为如果不仔细检查Python代码，这种差异可能会被忽略。

如图所示，左侧使用每个横断面的数据子集，右侧使用整个数据集。交互式模型生成的结果图显示，最佳拟合曲线相同，但可以看出95%置信区间不同。

接下来，我让它创建一个加性模型，该模型具有相同的响应变量，但仅包含MedSand和fTransect的主要效应，没有相互影响。此模型允许每个地段的截距有所不同，但强制L. acuta和MedSand之间的回归线斜率保持一致。ChatGPT-4o/Python能够运行加性模型，并生成所有正确的摘要统计信息（与一致的2个AIC单位差异）。然而，它使用了之前的交互模型来创建效应图。这是因为它使用了Python中的“sns.lmplot”命令来创建图形，而没有选项不包括相互作用。当我提出这个问题时，它能够使用“get_prediction”命令提供加性模型中的正确线条。然而，当我要求它创建类似于R中“plot(allEffects())”命令产生的“效应”图时，虽然它制作了这些图，但它完全搞砸了置信带和区间。首先，它使用了标准差来绘制误差条：

sns.pointplot(x='fTransect', y='L_acuta', data=df_clean, ci="sd")

绘制fTransect与L_acuta之间的点图，数据来源于df_clean，置信区间为标准差。

然后当我指出来之后，它使用了95%的置信区间。然而，生成的图表仍然与R中的效应图不符。当我询问时，原来它根据“实际观测到的数据”使用简单的点估计来计算95%的置信区间。我不太清楚它是如何做到这一点的，但这些结果并不是基于加法模型。经过进一步提示，它能够重现R中的fTransects效应图。当我要求它为MedSand变量做同样的事情时，它却出错了，于是我决定继续进行下一步。

对于横断变量的加法模型所生成的效应图，可以使用基于实际数据的 sns.pointplot 生成（左），或者使用加法模型中的 get_prediction 生成（右）。请注意左侧正确的估计均值但95%置信区间不正确。

然后，我让它运行一个模型，其中因变量相同，但预测变量是MedSand，另一个预测变量是MedSand和fTransect的交互项，而不包括fTransect的主要效应。这相当于允许不同穿越线的回归线斜率不同，但强制其截距相同。ChatGPT-4o能够设置并运行该模型，汇总统计信息均正确，并且AIC差异一致性为2个单位。然而，此时它仍然展示交互模型的图表。在再次提示之后，它明白了我想要利用模型来预测并绘制95%置信区间，因此能够用R重新创建效应图。

这是我们作业的核心部分。我让它运行了一个混合效应模型，响应变量与MedSand作为固定预测因子相同，但将fTransect作为随机截距。随机截距允许截距在三个不同样带上略有不同，假设这些样带是从潜在样带的随机抽样中来的，假设这些潜在样带有正态分布。在这里，我们没有假设随机斜率，也就是说我们假设每个样带的斜率是相同的。因此，这是之前没有交互项的模型的随机效应版本。幸运的是，ChatGPT-4o能够使用Python中的“statsmodels”库中的“mixed_linear_model”命令来运行这个混合效应模型，因此我们没有像我在上一篇文章中遇到的那种失败。它生成了带有正确最佳拟合线的效果图，并提供了接近正确的p值。它误解了我关于随机成分解释的方差比例的问题，这可能是由于我的措辞不清晰。在尝试为效果图添加95％置信区间时，发现Python中的statsmodels库没有针对混合模型的get_prediction命令，因此无法提供准确的95％置信区间。ChatGPT-4o试图仅使用固定效应提供置信区间带，但这些当然不准确。当被提示时，它试图通过将随机截距方差加入标准误差计算来引入随机效应，但结果仍然不准确，可能是因为使用的Wald近似。它试图通过实施自助法来生成更准确的95％置信区间，但很快就超时了，我继续前进。ChatGPT-4o能够正确计算组内相关系数（随机效应解释的方差比例），简称ICC，并能够正确解释它。它还能够提供边缘和条件R方值（分别仅包括固定效应，或包括固定效应和随机效应）。然而，它没有提供AIC值。

仅包含随机截距的混合效应模型使用ChatGPT-4o生成的效应图（左），以及使用R中的lme模型生成的效应图（右）。请注意，虽然左图的最佳拟合线是正确的，但其95%置信区间实际上是不正确的。

接下来，我让它拟合一个具有相同响应变量和预测变量的模型，并为 fTransect 添加一个随机截距，但同时也为 MedSand 添加了一个随机斜率。这使不同横断面中的回归线有所不同，但斜率的变化将被假设为具有正态分布，这三组横断面被看作是从该正态分布中随机抽取的样本。ChatGPT-4o 能够使用 MixedLM.from_formula 命令拟合并运行该模型，生成了正确的最佳拟合线和接近正确的 p 值。我通过手动计算两个模型的对数似然和参数数量来获取它们的 AIC 值。有趣的是，对于仅包含随机斜率的模型，2 个单位的差异消失了，因此 R 和 Python 计算出的 AIC 值相同。ChatGPT-4o 能够计算出组内相关系数，但在尝试获取效应图的准确 95% 置信区间时又失败了。它建议我使用 R，或作为 R 包装器的 pymer4 库，或实验性的 scikit-mixed-effects 库，但它没有这些库的访问权限。

接下来，我要求比较最后两个模型的ANOVA结果，它使用似然比检验来决定是否有必要在这个模型中加入随机斜率，或者随机截距已经足够。结果发现这是无法完成的，因为statsmodels库不支持混合模型的似然比检验（如果我错了，请纠正我）。相反，它使用AIC值和组内相关系数来解释为什么一个模型比另一个模型更好。这表明，AIC值的差异在R和Python之间是不同的。R中的带随机斜率和截距的模型的AIC值为604，而Python中的同一模型的AIC值为602。仅带随机截距的模型的AIC值为601。现在，在R中，这4个单位的差异使得很容易说带随机斜率的模型表现较差，应该选择带有随机截距的模型。然而，在Python中，1个单位的AIC差异不足以做出这样的判断。根据简约原则，当AIC值相差在2以内时，仍建议使用更简单的模型，但这一点可以讨论。当我要求它比较两个模型的ICC值时，它混淆了。然而，当我提醒ChatGPT-4它之前已经成功计算过这些值时，它能够提取并比较它们。根据现有信息，它建议不要在模型中包含随机斜率。然而，当我指出这将假设三个样带的斜率相同这一不合理假设时，它承认如果这是一个不合理的假设，那么可能最好还是加入随机斜率。我爱这种“客观性”！ :-) 它还说不用担心MedSand同时作为固定效应和随机斜率，在这种情况下，问题只会在每个样带只有一个MedSand值时出现。

最后，我让它绘制了带有随机斜率的模型的QQ图，就像我让学生做的那样，来思考模型假设的有效性。首先，它使用了未处理的残差，结果看起来非常奇怪，似乎不符合预期。在提示下，它将残差标准化，即用其标准差除以残差，得到的QQ图看起来更好，因为点沿着对角线排列，但中心部分的点仍然没有对齐。原始的R语言QQ图看起来完全不同，显示出偏离正态分布的迹象。区别在于R使用的是条件残差，这些残差同时包括固定效应和随机效应。不幸的是，Python中的MixedLM似乎没有生成包含固定效应和随机效应的条件预测的方法，因此无法得到条件残差。这是一个大问题，因为这样我们无法在ChatGPT-4o/Python中制作合适的QQ图，也就无法正确诊断线性混合模型。最好的情况是我们可以在ChatGPT-4o/Python中创建和运行模型，然后在R中重新编码和验证模型，但这听起来并不实际。到了这一步，我们可能还不如全部用R来做。

以下是QQ图，展示了ChatGPT-4o生成的具有随机斜率的混合效应模型，左边是标准化残差，右边是条件残差。注意条件残差图中明显的偏斜，而标准化残差图中却没有这种偏斜。

当我上传了我的学生用R创建的诊断QQ图时，ChatGPT-4能够检测到分布尾部的偏差偏离正态，但它还指出了我不太注意到的中心偏差——这些可能只是假象。它确实得出结论说分布“并非完美地正态分布”，但这其实是一个很弱的结论，因为没有任何事物是完全正态分布的。我在这张图上看到了非常明显的偏离正态性的偏差，这正是统计学展现艺术魅力的地方！

最后，我让它绘制了因变量 L.acuta 的分布情况。它正确地制作了一个直方图（使用的区间与学生在 R 中使用的不同），并能够检测到数据的高度右偏且非正态特性，以及这对之前进行的线性回归可能造成的影响。它还建议了解决方法，例如数据转换和广义线性混合模型（GLMM）。它正确地建议我们需要谨慎使用本次会话中创建的模型，而我会更直接地说这些模型毫无价值，让学生看到这一点。不幸的是，Python 中的 statsmodels 库目前仅支持少数贝叶斯实现的 GLMM 以及 scikit-mixed-effects 库，但这些都不适用于 ChatGPT-4o。

所以，总结一下，我发现今天虽然ChatGPT-4可以通过Python来拟合和运行线性混合效应模型，但它与在R中运行相同模型相比，存在一些实际限制。虽然它可以提供（大部分）正确的汇总统计量，但由于无法生成条件预测，它无法提供正确的95%置信带和区间。此外，它也无法直接提供AIC值或在模型之间进行似然比检验。最后，它无法基于条件残差生成正确的诊断图，这使得即使能够构建模型，也无法诊断它们。这些限制并不是ChatGPT-4本身的局限性，而是Python架构下底层库的局限性，这些问题可能随着时间的推移得到解决。此外，容易被ChatGPT-4生成的图表误导，尤其是在如何创建95%置信带这一点上。基于这些原因，我目前不建议使用ChatGPT-4/Python来创建、运行、诊断或使用线性混合模型。