1. 引言
近年来,在统计学与数据科学迅速发展的背景下,函数型数据分析成为应对高维、连续时间或空间观测数据挑战的重要工具。在实际应用中,研究者经常面临以函数或曲线形式呈现的变量,这类数据具有复杂的结构特征,传统的多元回归模型难以对其进行有效的建模与解释。
为了解决上述问题,1972年由Nelder和Wedderburn [1]首次引入函数型线性模型,分别研究了连续型和离散型响应变量与预测变量之间的关系。但是传统的广义线性模型在处理复杂数据时,模型的预设假设可能难以精确地契合数据的真实分布。针对这一问题,Scalla等人(1984) [2]提出了一种参数化连接函数的方法。该方法事先设定连接函数的形式,并通过数据来估计这些参数,从而调整模型的拟合效果,提高预测精度。然而,在实际应用层面,通过引入参数来优化连接函数的方式存在其固有的局限性。因为我们通常很难事先确定连接函数的具体形式。另外,由于数据的复杂性和独特性,标准的连接函数可能无法准确地拟合数据,从而限制模型的预测能力。而Weisberg和Welsh (1994) [3]结合了Wedderburn (1974) [4]提出的拟似然函数,在连接函数未知的广义线性模型拟合中,使用核平滑估计连接函数,再通过估计出的连接函数估计回归系数,利用交替迭代的方式,有效地解决了连接函数未知时的拟合难题。Chiou和Müller (1998) [5]对未知的方差函数和连接函数进行非参数估计,并将估计代入拟似然函数中,运用局部多项式拟合得到了连接函数估计和方差函数估计的一致性结果,以及回归系数的渐近极限分布。Cardot和Crambes等人[6]将总体最小二乘法推广到函数型线性模型中,提出了模型函数系数的光滑样条估计,并给出了该估计的渐近结果。Goldsmith等人[7]提出了函数型数据广义线性模型的快速拟合方法,通过对大量光滑特征向量进行投影,并利用惩罚样条回归估计系数函数。Fuchs等人[8]将具有函数型预测变量的标量响应的广义模型扩展到包含线性泛函交互项,系数函数的估计使用基展开式和最大似然估计法。
基于以上研究背景,本文在前人研究的基础上,构建了连接函数未知且含有函数型预测变量交互项的广义函数型线性模型,该模型不仅引入了未知的连接函数,还考虑到了函数型预测变量间的交互作用。另外,该模型还具有同时处理函数型预测变量和标量型预测变量的能力,能够更准确的阐述响应变量与预测变量之间的关系,具有重要的理论价值和现实意义。
2. 广义部分函数型线性模型
2.1. 模型介绍
假设有标量型响应变量
,
,
,
是影响
的函数型预测变量,其中
,
是
上的一个有界区间,
维向量
是影响
的标量型预测变量,响应变量
和预测变量
、
有以下关系:
(1)
其中,
是截距,
、
和
是两个函数型预测变量和交互项对应的未知回归系数函数,
是标量型预测变量
对应的未知回归系数,
是未知的连接函数,假设
是独立同分布的,且符合正态分布,均值为0,方差为
。
对于函数型预测变量,我们首先利用函数型主成分分析的方法对其进行降维。假设函数型预测变量
和标量型预测变量
均已中心化。通过Karhunen-Loeve展开和Mercer定理,将函数型预测变量
展开为
(2)
(3)
其中,
、
为函数型主成分得分,
、
为函数型主成分基,且为
空间中的标准正交基。
同理,将回归系数函数
、
展开为
(4)
(5)
(6)
其中,
、
是回归系数函数的系数,
是交互项的回归系数函数的系数。
将上述展开式代入模型(1)中,并将函数型预测变量在
处截断,其中
取决于样本容量
,且随着
而不断增加,得到
(7)
2.2. 回归系数与连接函数的估计
定义参数向量
(8)
(9)
(10)
对于未知的回归系数向量
和连接函数
,我们采取一种迭代交替更新的方法来进行估计。下面简单叙述迭代过程。
步骤一:假设连接函数为恒等映射,利用极大似然估计来估计参数向量
。样本的似然函数为
(11)
对似然函数取对数并对
求偏导,再令偏导为0,即求解下列方程,可求出参数向量的极大似然估计
。
(12)
步骤二:在估计出的参数向量
的基础上,利用局部线性核回归去估计连接函数
及其导数
。设核函数为
,带宽为
,定义与带宽有关的核函数为
,本文选取的核函数为高斯核函数。设函数型预测变量
和标量型预测变量
都属于一个紧集
,则
,为了估计
、
,构造任意
处的加权平方和为
(13)
化为矩阵形式,即
(14)
其中,
,
,
,
对
求导并令其为0,则
(15)
即解出
与
。
步骤三:运用步骤一的方法,将连接函数替换为估计出的连接函数
和
,通过求解含有
的方程来更新
,其中
表示迭代次数,得到
新的估计值。
步骤四:运用步骤二的方法,将参数向量替换为估计后的
,其中
,由此得到连接函数
和
新的估计值
和
。
步骤五:重复上述步骤,直到
收敛,停止迭代更新。
步骤六:得到最终的参数向量的估计值
,连接函数
的估计值
。
3. 数据模拟
我们考虑两个函数型预测变量、三个标量型预测变量以及二元0-1型响应变量的情况。对于函数型预测变量
,
,
,
,
,
可以为任意正整数,这里我们取
为500。
定义两个标准正交基为
,
和
,
。
,
,
生成两个函数型主成分得分
,
,
其中
,
,
,
,
,
,
,
,
。
所以有
, (16)
, (17)
函数型预测变量
和
的图像如图1和图2所示:
Figure 1. Graph of the functional predictor
图1. 函数型预测变量
的图像
Figure 2. Graph of the functional predictor
图2. 函数型预测变量
的图像
对于标量型预测变量
,满足
,
,
。假设标量型预测变量的回归系数理论值为
,
,
。
假设函数型预测变量的回归系数函数为
, (18)
, (19)
其中
,
,
,
,
,
,
,
,
。
对于交互项
,回归系数函数的系数记为
,且有
,
交互项的主成分基为
,满足
, (20)
因此交互项可以表示为
(21)
选择真实的连接函数为
,生成与函数型预测变量和标量型预测变量有关的二元响应变量
,且有
其中
服从伯努利分布。
我们选取样本量
为500。图3和图4展示了样本量为500情况下,估计出的
,
和其相应的95%置信区间带情况,可以看出样本量为500时,估计值已经很接近真实值。
Figure 3. Graph of the regression coefficient function
图3. 回归系数函数
的图像
Figure 4. Graph of the regression coefficient function
图4. 函数型预测变量
的图像
黑色曲线是回归系数函数估计
,
,灰色区域是他们的95%置信区间带,红色曲线是
,
的真实值。
图5展示了
为500时估计出的连接函数
,其中蓝色代表真实的连接函数,红色代表估计出的连接函数,可以看出,估计出的连接函数已经很接近真实值。
Figure 5. Plot of the true link function and the estimated link function
图5. 真实连接函数与估计出的连接函数图像
图6展示了
为500时交互项回归系数函数
的估计值与真实值作差的等高线图,可以看出其差值大部分在[−0.5, 0.5]之间。
Figure 6. Difference between the estimated surface and the true surface
图6. 估计曲面与真实曲面的差值
表1显示了
为500时
和其估计值
的值。通过表中数据可以看出,
为500时估计值接近真实值。
Table 1. Comparison between estimated and true values of scalar predictor
表1. 标量型预测变量的估计值与真实值对比
|
|
|
|
500 |
0.7055 (0.7071) |
0.5671 (0.5774) |
0.4898 (0.4472) |
4. 实例研究
本节所使用的数据源自R语言中的Tecator数据集。选取其中的198个肉类样本,其中蛋白质含量高于18%记为1,低于18%记为0,以此构成二元响应变量
。将198个肉类样本在850~1050波长范围内的光谱曲线数据作为第一个函数型预测变量
,将其导数作为第二个函数型预测变量
;将每个样本光谱曲线对应的MEA (均值)、SKE (偏度)、KUR (峰度)作为三个标量型预测变量,分别记为
、
、
。
图7和图8分别展示了两个函数型预测变量的图像。
Figure 7. Near-infrared spectral curves of meat samples
图7. 肉类近红外光谱曲线
Figure 8. First-order derivative curves of near-infrared spectra of meat samples
图8. 肉类近红外光谱一阶导数曲线
将数据代入模型(1)中,选取累计贡献率85%的函数型主成分作截断,则两个函数型预测变量的主成分基截断个数分别为
,
。标量型预测变量回归系数的结果如表2所示。实验结果说明MEA (均值)、SKE (偏度)、KUR (峰度)均与肉类蛋白质含量呈负相关,即这三个值越大,则对应的肉类蛋白质含量越低。
Table 2. Estimated values of scalar predictors
表2. 标量型预测变量的估计值
|
Estimate |
Std.Error |
t value |
Pr (>|t|) |
|
−0.2528 |
0.1673 |
−1.5111 |
0.13241 |
|
−1.0363 |
0.3350 |
−3.093 |
0.00228** |
|
−1.9254 |
0.6303 |
−3.055 |
0.00257** |
图9和图10分别是两个函数型预测变量回归系数函数的图像。
Figure 9. Regression coefficient of the functional predictor X1
图9. 函数型预测变量X1的回归系数
Figure 10. Regression coefficient of the functional predictor X2
图10. 函数型预测变量X2的回归系数
结论表明肉类的近红外光谱曲线与其蛋白质含量有明显的相关性。由图9可以看出肉类在波长为约850~920 nm时的近红外光吸收与其蛋白质含量呈正相关,在波长为920~1050 nm处的近红外光吸收与其蛋白质含量呈负相关。由图10可以看出在920~990 nm左右的波长段,近红外光谱数据的一阶导数与样本肉类蛋白质含量呈正相关,吸收斜率变化越大,对应肉类样本中的蛋白质含量越高,在其余波长段则呈负相关。因此肉类加工企业、品质检测机构等可利用此模型快速估测肉类样本蛋白质含量的高低,有效降低检测成本。
图11和图12分别展示了近红外光谱曲线及其导数交互项的回归系数的等高线图和三维曲面图。从图中可看出二者之间明显存在交互效应,对肉类蛋白质含量高低判断产生影响。
Figure 11. Contour plot of the interaction term
图11. 交互项的等高线图
Figure 12. 3D surface of the interaction term
图12. 交互项的三维曲面
5. 模型对比
为验证本文所提模型的有效性,我们设计了一项对比实验,具体包含两个目标:(1) 检验在连接函数未知时模型的性能;(2) 评估在模型中考虑函数型预测变量交互项的必要性。
为此,我们将所提模型与以下两个基准模型进行比较。
模型1 [9]:连接函数已知,且包含交互项;
其中
已知;
模型2 [10]:连接函数已知,但不包含交互项。
模型评估采用以下两类指标:首先,在拟合优度方面,选用赤池信息准则(AIC)、残差与决定系数(R-squared)三个指标;其次,在泛化能力方面,通过五折交叉验证计算各模型的预测性能指标。其中AIC的值和残差值越小说明模型的拟合效果和泛化能力越好,而R-squared的值越大说明模型拟合效果越好。如表3所示,本文提出的模型在多项指标上均优于对比模型:其AIC值与残差和更小,R-squared值更高,且预测准确率也更具优势。这些结果一致表明,本文模型具有更优的拟合性能。
Table 3. Model comparison results
表3. 模型对比结果
模型 |
AIC |
R-squared |
残差 |
准确率 |
本文模型 |
271.5259 |
0.9258 |
1.2942 |
97.03% |
模型1 |
279.8883 |
0.9014 |
1.8448 |
94.35% |
模型2 |
292.0326 |
0.6731 |
2.2093 |
90.45% |
6. 结论
本研究旨在建立一类能够刻画函数型预测变量交互作用、且连接函数形式未知的广义函数型线性模型,并应用于肉类蛋白质含量的分析。实证结果表明:肉类样本在850~1050 nm波长区间内的近红外光谱曲线及其一阶导数,均与蛋白质含量显著相关;此外,光谱曲线的均值、偏度与峰度等统计量亦与蛋白质含量呈负相关。因此,在评估肉类蛋白质含量时,可综合依据上述光谱及衍生特征进行判定。
NOTES
*通讯作者。