随机效果和受罚样条线是一回事

很长一段时间以来，我一直对某些事情感到好奇。我随随便便就读到教科书，包装文档和推文中的一个事实是：随机效果和受罚的平滑样条线是一回事。听起来如此深刻和启发。这是什么意思？它们是一样的吗？我错过了什么深度统计诊断？

我花了几个月的时间来尝试理解这种等效性。我无法给您完整的数学处理方法，但是我有主旨，可以为您指出这些等式。在这篇文章中，我将尝试强调两者之间的联系。

混合效应模型使用部分合并在总体均值（完全合并）和单个均值（不合并）之间取得平衡。

可悲的是，我觉得我的职业生涯已经随着这个模因的创建而达到顶峰pic.twitter.com/5ilRFonsy7

-埃里克·佩德森（@ericJpedersen）2020年8月12日

让我们回顾一下这些含义。混合效应模型显然是多年来博客的主要关注点，用于估计“随机”或“变化”效应。这是经典的等式设置：

\ [\ mathbf {y} = \ mathbf {X \ beta} + \ mathbf {Zb} + \ mathbf {\ epsilon} \\\ mathbf {b} \ sim \ textsf {Normal}（0，\ sigma_b）\\ \ mathbf {\ epsilon} \ sim \ textsf {Normal}（0，\ sigma_y）\\\ mathbf {X}：\ textrm {固定效果模型矩阵} \\\ mathbf {Z}：\ textrm {随机效果模型矩阵} \\\ sigma_b，\ sigma_y：\ textrm {方差组件} \\\ sigma_b：\ textrm {发生魔术的地方} \\\]这里的魔术是σb，因为它把b下所有b的所有效应都联系在一起共同的分布。如果将σb替换为固定数（例如10），则b中的所有影响将是独立的并且彼此不知道：组之间将不会汇集任何信息。

考虑来自Gelman和Hill（2007）的the数据集示例。 Min测量在明尼苏达州进行。我们想估算每个县的平均ra测量值。我们有反复出现的措施情况，有些县比其他县有更多的观察结果。我们使用带有brms的贝叶斯混合效应模型来估计县估计的人口分布，而县级估计是随机变化的影响。它们是从随机分布中得出的，我们从数据中估算出其规模。

库（tidyverse）theme_set（theme_grey（base_size = 14））库（brms）radon＆lt;-rstanarm :: radon b_radon＆lt;-brm（log_radon〜1 +（1 |县），radon，family =高斯，文件=＆ ＃34; rad＆＃34;）b_radon＃＆gt;家庭：高斯＃＆gt;链接：mu =身份； sigma =身份＃＆gt;公式：log_radon〜1 +（1 |县）＃＆gt;数据：ra（观察数：919）样本：4条链，每个链的iter = 2000；预热= 1000;薄= 1; ＃＆gt;预热后样本总数= 4000＃＆gt; ＃＆gt;组级效果：＃＆gt; 〜county（级别数：85）＃＆gt;估计估计错误1-95％CI u-95％CI Rhat Bulk_ESS Tail_ESS＃＆gt; sd（截距）0.30 0.05 0.22 0.40 1.00 1782 2894＃ ＃＆gt;人口层面的影响：估计估计错误1-95％CI u-95％CI Rhat Bulk_ESS Tail_ESS＃＆gt;拦截1.35 0.05 1.26 1.45 1.00 2749 3198 ＃＆gt;家庭特定参数：估计估计错误1-95％CI u-95％CI Rhat Bulk_ESS Tail_ESS＃＆gt; sigma 0.77 0.02 0.73 0.80 1.00 7374 3105＃＆gt; ＃＆gt;使用采样（NUTS）抽取样本。对于每个参数，Bulk_ESS＃＆gt;和Tail_ESS是有效的样本量度量，而Rhat是潜在的＃＆gt;。拆分链上的比例缩减因子（收敛时，Rhat = 1）。

我们可以在模型估计值旁边绘制观察到的县均值。首先，我进行了一些争论，以便计算观察到的均值和估计均值之间的差异，以备后用。

radon_aug＆lt;-radon％＆gt;％＃添加ns并表示group_by（County）％＆gt;％mutate（observed_mean =平均值（log_radon），county_n = n（））％＆gt;％ungroup（）％＆gt;％＃add拟合值tidybayes :: add_fitted_draws（b_radon）％＆gt;％变异（observed_minus_model = observed_mean-.value）％＆％;％＃汇总拟合值ggdist :: mean_qi（.value，observed_minus_model）radon_aug $类型＆lt;-＆＃34;混合模型估算值don $类型＆lt;-＆＃34;观察到的＆＃34; ggplot（radon_aug）+ aes（x = fct_infreq（县），y = log_radon）+ stat_summary（aes（颜色=类型，shape =类型），data = radon，fun =平均值，geom =＆＃34; point＆＃34; ）+ geom_point（aes（y = .value，color = type，shape = type））+＃要在图形中包含0 geom_blank（aes（y = 0））+ labs（x =＆＃34; county（in按样本大小递减顺序）＆＃34;，y =＆＃34; log（radon）＆＃34;）+ geom_hline（yintercept = fixef（b_radon）[1]）+ scale_color_manual（值= c（＆＃34; blue＆＃34;，＆＃34; grey40＆＃34;））+实验室（color = NULL，shape = NULL）+主题（axis.text.x = element_blank（），axis.ticks.x = element_blank（），面板.grid.major.x = element_blank（），panel.grid.minor.x = element_blank（），legend.title = element_blank（），legend.position =＆＃34; top＆＃34;，legend.direction =＆＃ 34;水平tal＆＃34; ，legend.justification =＆＃34; left＆＃34; ，）

我们看到了部分池化的经典示例。对于有很多观测的县，估计均值几乎无法调整。对于没有数据的县，估算值将推向总体平均值（以上概述中的截距）。

下图显示了观察到的均值与估计均值之间的差，从上图中的蓝色正方形中减去了灰色三角形。

radon_aug％＆gt;％ungroup（）％＆gt;％差异（county，county_n，Observed_minus_model）％＆％;％ggplot（）+ aes（x = county_n，y = observed_minus_model）+ geom_point（alpha = .5）+实验室（x =＆＃34;县中的观测数，y =＆观察到的均值-估计均值）

平滑=随机效果主张的争论在于，wejust所做的只是平滑情况。从某种意义上说，这些随机效应是平滑的固定效应。

现在，让我们看一下mgcv中的广义加性模型，以展示受罚的平滑样条曲线。那真是令人mouth目结舌，但是从本质上讲，加性模型就像是标准线性模型的平滑扩展包。我们仍在进行回归，但是我们有了一些新的语法，并且我们的模型现在可以更轻松地进行非线性关系了。

我将介绍一个基本示例，说明样条曲线的基础函数如何加权以近似非线性趋势，但这不会成为完整的教程。其他人已经对附加模型或mgcv软件包进行了视频介绍。我首先从语言学家的教程中学习它们，然后从mgcv教科书中学习，但是还有其他在线资源。

我们使用mcycle数据集来提供模拟摩托车事故中的头部加速度。我们将拟合模型，从中绘制平滑度，然后通过模型所做的工作。

库（mgcv）mcycle＆lt;-MASS :: mcycle％＆gt;％tibble :: rowid_to_column（）＃拟合模型gam_20＆lt;-gam（accel〜1 + s（times，bs =＆＃34; cr＆＃34; ，k = 20），data = mcycle，method =＆＃34; REML＆＃34;）mcycle $ .fitted＆lt; -fitted（gam_20）ggplot（mcycle）+ aes（x = times，y = accel）+ geom_point（ alpha = .5）+ geom_line（aes（y = .fitted），颜色=＆＃34; blue＆＃34;）+实验（x =＆＃34;撞击后的时间[ms]＆＃34;，y =＆ ＃34;加速度[g]＆＃34;）

我们告诉gam（）在时间预测变量（s（times，...））上使用截距项和平滑项来估计加速度。具体来说，我们使用三次回归样条曲线基础（bs =＆＃34; cr＆＃34;）创建了平滑的曲线，其中k = 20-1条曲线。我们的模型通过将称为基础函数的较小组成部分相加来估计函数，并且定义这些组成部分的空间是基础。这些基础函数经过加权和求和以产生称为样条的平滑趋势。样条线的名称是受样条线的启发而产生的，这些样条线是可弯曲的木板条，可以加重并固定在适当的位置以形成优美的曲线。

重申一下，从概念上讲，我们将时间预测器分解为一堆单独的摆动线（基本函数），并对它们进行加权和求和，以近似得出一些非线性函数。关于正交多项式的Mypost说明了相同的原理，但是具有多项式基函数。RichardMcElreath在统计重新思考过程中提供了贝叶斯模型中30分钟的友好样条。我从他的描述中了解到的一条话是，通过样条曲线，我们可以用一组“综合”预测变量来代替预测变量，例如时间。

拉动摆动的一种简单方法是使用模型矩阵。我们有20个用于拦截的列和19个基本函数。

model.matrix（gam_20）％＆gt;％tibble :: trunc_mat（width = 72）＃＆gt; ＃说明：dbl [，20] [133 x 20]＃＆gt; `（拦截）``s（次数）.1``s（次数）.2`s（次数）.3`s（次数）.4`＃＆gt; ＆lt; dbl＆gt; ＆lt; dbl＆gt; ＆lt; dbl＆gt; ＆lt; dbl＆gt; ＆lt; dbl＆gt; ＃＆gt; 1 1 -0.143 -0.165 -0.154 -0.311＃ 2 1 -0.0317 -0.193 -0.133 -0.293＃ 3 1 0.292 -0.270 -0.0716 -0.242＃ 4 1 0.492 -0.309 -0.0357 -0.208＃ 5 1 0.669 -0.332 -0.00603 -0.175＃ 6 1 0.933 0.0792 -0.0455 -0.0280＃ 7 1 0.808 0.278 -0.0969 -0.0106＃ 8 1 0.730 0.383 -0.122 -0.00387＃ 9 1 0.286 0.851 -0.182 0.00169＃＆gt; 10 1 0.123 0.957 -0.143 -0.0130＃ ＃...还有123行，还有15个变量：`s（times）.5`＆lt; dbl＆gt ;,＃＆gt; ＃`s（times）.6`＆lt; dbl＆gt ;,`s（times）.7`＆lt; dbl＆gt ;,`s（times）.8`＆lt; dbl＆gt ;,＃＆gt; ＃`s（times）.9`＆lt; dbl＆gt ;,`s（times）.10`＆lt; dbl＆gt ;,`s（times）.11`＆lt; dbl＆gt ;,＃＆gt; ＃`s（times）.12`＆lt; dbl＆gt ;,`s（times）.13`＆lt; dbl＆gt ;,`s（times）.14`＆lt; dbl＆gt ;,＃＆gt; ＃`s（times）.15`＆lt; dbl＆gt ;,`s（times）.16`＆lt; dbl＆gt ;,`s（times）.17`＆lt; dbl＆gt ;,＃＆gt; ＃`s（times）.18`＆lt; dbl＆gt ;,`s（times）.19`＆lt; dbl＆gt;

为了可视化矩阵，我使用个人Rpackage中的辅助函数在ggplot2中绘制矩阵。我们看到的是在x轴上的时间，在截距和每个基本函数的一行上。

＃帮助函数绘制矩阵ggmatplot的行<-tjmisc :: ggmatplot＃帮助函数标记在theme_grey（）图上的annotate_grey＆lt;-tjmisc :: annotate_label_grey ggmatplot（cbind（mcycle $ times，model.matrix（gam_20 ）），x_axis_column = 1）+ annotate_grey（＆＃34; intercept＆＃34;，0，1.02，size = 5）+ annotate_grey（＆＃34; individual \ nbasis \ nfunctions＆＃34;，0，.16，size = 5）+ expand_limits（y = 1.2）+实验（x =＆＃34; time [ms]＆＃34;，title = NULL）

现在，我们可以通过乘以模型系数来对其加权。在此，我们使用diag（coef（））技巧来防止将加权的预测变量求和。

weighted_coefs＆lt;-model.matrix（gam_20）％*％diag（coef（gam_20））ggmatplot（cbind（mcycle $ times，weighted_coefs），x_axis_column = 1）+ annotate_grey（＆＃34;加权拦截＆＃34;，35， -40，size = 5）+ annotate_grey（加权基函数＆＃34;，0，26，size = 5）+实验室（x =＆＃34; time [ms]＆＃34;，title = NULL ）

现在，我们可以看到数据集中的两个主要拐点。大约20毫秒和30毫秒的基函数变得非常活跃，以便在那时将样条曲线推离0。

如果将这些线加在一起，则得到回归线（intercept加上平滑样条线）。

ggmatplot（cbind（mcycle $ times，weighted_coefs），x_axis_column = 1，n_colors = 1）+ stat_summary（aes（group = 1），color =＆＃34; maroon＆＃34;，fun = sum，geom =＆＃34; line＆＃34;，size = 1.5）+ annotate_grey（＆＃34; sum＆＃34;，10，-70，size = 5，color =＆＃34; maroon＆＃34;）+实验室（x =＆＃34;时间[ms]＆＃34;，标题= NULL）

到目前为止，我们的图显示了具有基函数的回归，但是平滑样条曲线又向前延伸了一步：它们惩罚了摆动性，以防止过度拟合。想法如下：对于上面的演示，选择20维样条曲线的基础（19条曲线，因为显然删除了1条，以便于识别）。但是那20个数字是从哪里来的呢？稀薄的空气。如果我们指定尺寸为50怎么办？那是50个预测变量（1个拦截函数和49个基本函数）。用这种方法过拟合数据真的很容易吗？

gam_50＆lt;-gam（accel〜s（times，bs =＆＃34; cr＆＃34;，k = 50），data = mcycle，method =＆＃34; REML＆＃34;）mcycle $ .fitted50＆lt;- fit.values（gam_50）ggplot（mcycle）+ aes（x =次，y = accel）+ geom_point（alpha = .5）+ geom_line（aes（y = .fitted，color =＆＃34; 20＆＃34;） ，size = 1）+ geom_line（aes（y = .fitted50，color =＆＃34; 50＆＃34;），size = 1）+ labs（x =＆＃34;撞击后的时间[ms]＆＃34; ，y =＆＃34;加速度[g]＆＃34;，颜色=＆＃34;尺寸＆＃34;）

嗯，它们看起来几乎没有什么不同。没有过度拟合。这是怎么回事？我已经说过了摇摆不定正在受到惩罚。

在后台，该模型试图平衡两个相互竞争的目标：一方面，我们希望最大限度地拟合数据。在线性回归中，此目标等于使平方误差的总和最小化。另一方面，我们要使摆动性（过度拟合）最小化。钝化的平滑样条线，这是通过首先指定一个惩罚矩阵来完成的，该惩罚矩阵定义了该样条线的摆动性。这两个功能在以下等式中相互抵触：

\ [\ begin {align *} \ mathbf {\ hat {β}}＆amp; = \ operatorname {arg min} _ \ mathbf {β} \ \ | \ mathbf {y} − \ mathbf {Xβ} \ | ^ 2 + \ lambda \ mathbf {β} ^ \ intercal \ mathbf {Sβ} \\＆amp; \ \ text {（want} \ mathbf {\ hat {β}} \ text {可使拟合误差和摆动最小化）}} \\ \ | \ mathbf {y} − \ mathbf {Xβ} \ | ^ 2＆amp ;： \ text {平方误差之和（最小化误差以提高拟合度）} \\ \ lambda \ beta ^ \ intercal \ mathbf {Sβ}＆amp; ：\ text {wiggliness惩罚} \\ \ mathbf {Xβ}＆amp ;： \ text {样条基础乘以权重} \\\ mathbf {S}＆amp ;： \ text {惩罚矩阵（定义样条的wiggliness）}} \\ \ lambda＆amp ;： \ text {平滑度参数（增加罚则以增加强度）} \\\ end {align *} \]不必担心这里的确切数学原理：只需要知道错误现在与摇摆不定配对即可，而摇摆不定是由惩罚矩阵S和平滑度参数λ控制。是的，摆动是技术术语。我基于西蒙·伍德的幻灯片中的方程式，该方程式使用“ fit-wiggliness权衡”一词。

出于我们的目的，我们不必担心惩罚矩阵。我会简单地描述它。对于该模型，通过使用估计的样条曲线的二阶导数来定义摆动。一阶导数测量曲线沿x的斜率/陡度，二阶导数测量斜率/陡度随x的变化量。因此，摆动是斜率的变化，惩罚矩阵为与此相关的每个模型系数提供了惩罚摆动性。出色的免费提供将将惩罚矩阵绘制为热图。 x和y代表模型系数（基本函数的权重），因此沿主对角线我们看到对每个系数都施加了惩罚。在两个非对角线中，我们看到相邻的基函数的权重共同受到惩罚或惩罚。

我认为最严厉的惩罚是第5和第6个基本函数，因为这两个基本函数覆盖了数据集中的大多数行，但是我对这种解释并不是100％充满信心。

要查看实际的过度拟合，我们可以通过使用固定的（fx = TRUE）回归样条曲线在上述样条比较图中禁用平滑惩罚。现在，模型的主要目标是使误差最小化，而50维样条曲线为模型提供了许多很多自由度。

gam_20_fx＆lt;-gam（accel〜s（times，bs =＆＃34; cr＆＃34;，k = 20，fx = TRUE），data = mcycle，method =＆＃34; REML＆＃34;）gam_50_fx＆lt; -gam（accel〜s（times，bs =＆＃34; cr＆＃34;，k = 50，fx = TRUE），data = mcycle，method =＆＃34; REML＆＃34;）mcycle $ .fitted20_fx＆lt; -fit.values（gam_20_fx）mcycle $ .fitted50_fx＆lt;--fit.values（gam_50_fx）ggplot（mcycle）+ aes（x =次，y =加速度）+ geom_point（alpha = .5）+ geom_line（aes（y = .fitted20_fx，color =＆＃34; 20＆＃34;），size = 1）+ geom_line（aes（y = .fitted50_fx，color =＆＃34; 50＆＃34;），size = 1）+实验室（x = ＆＃34;撞击后时间[ms]＆＃34;，y =＆＃34;加速度[g]＆＃34;，颜色=＆＃34;尺寸＆＃34;）

因此，当我们禁用惩罚时，50维样条线可以在整个位置自由摆动。

平滑参数λ是超参数。它控制样条系数（基函数权重），并根据数据进行估算。我们可以手动设置λ，然后向上摇动它。在这种情况下，模型将尝试找到减少建模误差的最小摆动曲线：一条直线。

gam_20_sp＆lt;-gam（accel〜s（times，bs =＆＃34; cr＆＃34;，k = 20，sp = 10000000），data = mcycle，method =＆＃34; REML＆＃34;）mcycle $。 fit20_sp＆lt;-拟合值（gam_20_sp）ggplot（mcycle）+ aes（x =次，y =加速度）+ geom_point（alpha = .5）+ geom_line（aes（y = .fitted，color =＆＃34; estimated＆ ＃34;），size = 1）+ geom_line（aes（y（。= fitted20_fx，color =＆＃34; no smoothing＆＃34;），size = 1）+ geom_line（aes（y（y = .fitted20_sp，color =＆＃ 34; 10,000,000＆＃34;），大小= 1）+实验（x =＆＃34;撞击后的时间[ms]＆＃34;，y =＆＃34;加速度[g]＆＃34;，颜色= ＆＃34; lambda＆＃34;）

我们需要某种方式来讨论发生了多少平滑。一方面，我们可以将每个基函数视为独立的预测变量，在拟合曲线时会用尽所有的自由度。另一方面，我们可能会对基函数权重进行过多的惩罚，以至于它们产生一条直线，因此，这组预测变量的作用就像单个预测变量一样。也就是说，他们有效地估算了其中只有1个自由度的作用值的曲线。实际上，这就是mgcv描述模型平滑度的方式：它报告了每个平滑度背后的有效（或估计）自由度（EDF）。

如果我们看一下模型摘要，就会发现我们的20维基准平滑度的EDF为11.78（请参见平滑项的近似重要性下的edf）。

摘要（gam_20）＃＆gt; ＃＆gt;家庭：高斯＃＆gt;链接功能：身份＃＆gt; ＃＆gt;公式：＃＆gt;加速〜1 + s（次，bs =＆＃34; cr＆＃34 ;, k = 20）＃＆gt; ＃＆gt;参数系数：估计标准误差t值Pr（＆gt; | t |）＃＆gt; （截取）-25.546 1.956 -13.06＆lt; 2e-16 ***＃＆gt; --＃＆gt;签名。代码：0＆＃39; ***＆＃39; 0.001＆＃39; **＆＃39; 0.01＆＃39; *＆＃39; 0.05＆＃39;。＆＃39; 0.1＆＃39; ＆＃39; 1＃＆gt; ＃＆gt;平滑项的近似含义： edf Ref.df F p值＃＆gt; s（次）11.78 13.97 33.93＆lt; 2e-16 ***＃＆gt; --＃＆gt;签名。代码：0＆＃39; ***＆＃39; 0.001＆＃39; **＆＃39; 0.01＆＃39; *＆＃39; 0.05＆＃39;。＆＃39; 0.1＆＃39; ＆＃39; 1＃＆gt; ＃＆gt; R-sq。（adj）＝ 0.782，说明的偏差＝ 80.1％。 --

......