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

2021-03-02 10:32:28

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

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

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

可悲的是,我觉得我的职业生涯已经随着这个模因的创建而达到顶峰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<-rstanarm :: radon b_radon<-brm(log_radon〜1 +(1 |县),radon,family =高斯,文件=& #34; rad")b_radon#>家庭:高斯#>链接:mu =身份; sigma =身份#>公式:log_radon〜1 +(1 |县)#>数据:ra(观察数:919)样本:4条链,每个链的iter = 2000;预热= 1000;薄= 1; #>预热后样本总数= 4000#> #>组级效果:#> 〜county(级别数:85)#>估计估计错误1-95%CI u-95%CI Rhat Bulk_ESS Tail_ESS#> sd(截距)0.30 0.05 0.22 0.40 1.00 1782 2894# #>人口层面的影响:估计估计错误1-95%CI u-95%CI Rhat Bulk_ESS Tail_ESS#>拦截1.35 0.05 1.26 1.45 1.00 2749 3198 #>家庭特定参数:估计估计错误1-95%CI u-95%CI Rhat Bulk_ESS Tail_ESS#> sigma 0.77 0.02 0.73 0.80 1.00 7374 3105#> #>使用采样(NUTS)抽取样本。对于每个参数,Bulk_ESS#>和Tail_ESS是有效的样本量度量,而Rhat是潜在的#>。拆分链上的比例缩减因子(收敛时,Rhat = 1)。

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

radon_aug<-radon%>%#添加ns并表示group_by(County)%>%mutate(observed_mean =平均值(log_radon),county_n = n())%>%ungroup()%>%#add拟合值tidybayes :: add_fitted_draws(b_radon)%>%变异(observed_minus_model = observed_mean-.value)%&%;%#汇总拟合值ggdist :: mean_qi(.value,observed_minus_model)radon_aug $类型<-"混合模型估算值don $类型<-"观察到的" ggplot(radon_aug)+ aes(x = fct_infreq(县),y = log_radon)+ stat_summary(aes(颜色=类型,shape =类型),data = radon,fun =平均值,geom =" point" )+ geom_point(aes(y = .value,color = type,shape = type))+#要在图形中包含0 geom_blank(aes(y = 0))+ labs(x =" county(in按样本大小递减顺序)",y =" log(radon)")+ geom_hline(yintercept = fixef(b_radon)[1])+ scale_color_manual(值= c(" blue"," grey40"))+实验室(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 =" top",legend.direction =&# 34;水平tal" ,legend.justification =" left" ,)

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

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

radon_aug%>%ungroup()%>%差异(county,county_n,Observed_minus_model)%&%;%ggplot()+ aes(x = county_n,y = observed_minus_model)+ geom_point(alpha = .5)+实验室(x ="县中的观测数,y =&观察到的均值-估计均值)

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

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

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

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

库(mgcv)mcycle<-MASS :: mcycle%>%tibble :: rowid_to_column()#拟合模型gam_20<-gam(accel〜1 + s(times,bs =" cr" ,k = 20),data = mcycle,method =" REML")mcycle $ .fitted< -fitted(gam_20)ggplot(mcycle)+ aes(x = times,y = accel)+ geom_point( alpha = .5)+ geom_line(aes(y = .fitted),颜色=" blue")+实验(x ="撞击后的时间[ms]",y =& #34;加速度[g]")

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

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

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

model.matrix(gam_20)%>%tibble :: trunc_mat(width = 72)#> #说明:dbl [,20] [133 x 20]#> `(拦截)``s(次数).1``s(次数).2`s(次数).3`s(次数).4`#> < dbl> < dbl> < dbl> < dbl> < dbl> #> 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#> 10 1 0.123 0.957 -0.143 -0.0130# #...还有123行,还有15个变量:`s(times).5`< dbl&gt ;,#> #`s(times).6`< dbl&gt ;,`s(times).7`< dbl&gt ;,`s(times).8`< dbl&gt ;,#> #`s(times).9`< dbl&gt ;,`s(times).10`< dbl&gt ;,`s(times).11`< dbl&gt ;,#> #`s(times).12`< dbl&gt ;,`s(times).13`< dbl&gt ;,`s(times).14`< dbl&gt ;,#> #`s(times).15`< dbl&gt ;,`s(times).16`< dbl&gt ;,`s(times).17`< dbl&gt ;,#> #`s(times).18`< dbl&gt ;,`s(times).19`< dbl>

为了可视化矩阵,我使用个人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%。 --

......