年 级:
成 绩:
指导教师:
多元回归分析法下的航班风险数据评估
应用统计学
指导教师
依托课程 统计软件
摘要:航班起飞风险评估能够对影响航班起飞的因素进行有效地识别、分析、与评估,
是保证航班起飞安全的重要手段与方法。本文通过收集可能影响航班总的风险值的数据,首
先通过散点图对数据进行可视化,结合相关系数选出对航班运行风险总值有着重要影响的因
素,进行一元回归建模,并探索是否存在多项式的关系。接下来对所有模型进行全模型多元
回归,通过所有子集法变量选择选出决定系数比较大的模型。最终确定多元线性模型来评估
各个因素对航班总影响值的影响因素。接下来我们对最终,模型进行回归诊断,确定最终模
型。
关键词:多元回归分析 所有子集法 回归诊断
背景:
航空公司的安全隐患广泛存在于各个运行过程中,例如机组配合程度、机长
经验能力、起飞机场天气等各个环节,而每个环节又都具有自己所属专业特点的
生产要素、组织结构和业务流程,任何一个环节或流程上的失误都有可能引起不
安全事件的发生。鉴于航班起飞风险因素的复杂多样且关系不明,仅凭管理经验
已经无法做出正确的风险评估,这就需要运用专门的方法根据航班起飞数据做出
综合评估,为风险管理提供决策支持。
表格 1 是每种风险因素的风险评级以及总的航班风险级别所对应的符号,表
格 2 是对风险等级重要程度的一些补充说说明。共分为七的等级,从低到高依次
为无风险、风险程度低、风险程度一般、有一定风险、风险程度中等,风险程度
高、不能接受的风险的风险程度。我们将对所得数据进行建模预测,最终对模型
进行分析,找出对总的风险模型影响最大的模型,最终提出改进建议。
数据介绍:
变量
y
1x
2x
3x
表 1 基本变量符号说明
变量含义
航班运行总风险值
机组配合程度
机长经验能力
副驾驶经验能力
变量
4x
5x
变量含义
机组疲劳程度
机组压力程度
变量
8x
9x
6x
机龄
10x
7x
MEL/CDL 风险
11x
变量含义
起飞机场天气风险
起飞机场条件风险
着陆机场天气风险
着陆机场条件风险
变量
12x
13x
14x
15x
变量含义
航路风险
航路备降机场天气风
航路备降机场条件风
特殊运行种类(二放、
险
险
极地、etops 等)
评分尺度
重要程度(ri)
说明
表 2 评分尺度(表)
1-5 表示风险较低
5-8 表示风险中等
8-10 表示风险等高
1
2
4
5
6
8
10
无风险
风险程度很低
风险程度一般
有一定风险
风险程度中等
风险程度高
此种风险程度不能接受
由于变量较多,自编函数对数据进行可视化,如图 1:
图 1 各个变量与因变量的散点图
1、模型的初步建立
初步画出相关系数矩阵并分析
图 2 Y 与 15 个自变量的相关系数矩阵
从图 2 的相关系数矩阵中我们能够看出航班运行总风险值(y)与其他其他
十五个自变量的相关性,显然的我们可以看出航班运行总风险值(y)的值与一
半上以上的自变量的相关性不到 0.5,仅有 x1、x2、x5、x9、x11、x13 六个变
量的相关性超过 0.5,且与起飞机场条件风险(x13)的相关系数最大达到 0.99,
下面我们画出其散点图观察其与航班运行总风险值(y)的具体趋势,以此来确
定初步的拟合模型。
可能相关共线性:x5 和 x6(0.94)、x3 和 x11(0.9)、x1 和 x3(0.84)、
x1 和 x11(0.83)。
图 3 单因素散点图
从上述相关系数中我们知道航班运行总风险值(y)与航路备降机场天气风
险(x13)的相关系数比较大为 0.99,同时在图 3 中我们也看出在画出散点图之
后航班运行总风险值(y)与航路备降机场天气风险(x13)呈正相关的线性关系,
故我们用线性模型拟合数据,把航路备降机场天气风险(x13)做自新变量对 y
做一元回归模型,并试图探索其中的多项式因素对模型进行改善。
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.002865 0.034106 0.084 0.933
x13 0.989576 0.005496 180.040 <2e-16 ***
计算结果如下:
表 3 模型求解的主要参数
模型
1、仅仅有 x13
2、在模型 1 的基础加二次项
2、在模型 1 的基础加三次项
Intercept
X13
加入高次项
0.002865
0.989576
R2
分别得出以下三个模型:
0.9862
0.16934
0.909661
0.007210
0.9866
0.1417634
0.9388649
0.000489
0.9867
y
=
0.003 0.99
+
x
13
,
2
R =
0.9862.
y
=
x
0.17 0.91
13
+
+
0.007
2
x
13
,
2
R =
0.9866.
y
=
0.15 0.94
+
x
13
+
0.004
3
x
13
,
2
R =
0.9867.
可见原有数据的相关性非常高,我们模型所得出的决定系数均达到了 0.98
以上,有一个非常好的效果。当然我们也看到在原始线性模型中假如二次项,三
次项之后,决定系数 2R 均有所提高,加入二次项后提高了 0.004,加入了三次
项之后提高了 0.005。通过以上分析我们虽然看到模型的拟合效果有所提高,但
是考虑到原有线性的拟合效果已经非常满意,为了降低模型的复杂程度,我们最
终选用一元线性的模型来拟合航班运行总风险值(y)与航路备降机场天气风险
(x13)。
下面对原一元线性与加入高次项所得模型的预测曲线进行可视化比较。然后
对模型进行线性模型的基本假设验证。从而对模型的进一步修正提供建议。
图 4 三个模型的拟合图像
图 5 三个模型的简单验证
以上 12 个图形分别初步对所得三个模型进行检验,每一行的第一个图描述
的是估计值与残差的关系,可以看出所得残差基本稳定,符合基本线性假设。第
二幅图是对残差的正态性进行验证,我们可以看出有大量的点落在 45 度直线之
外,故应该判断不符合基本的正态性假设。第三个图是看同方差性,我们可以看
出,散点都分布在红线两侧,我们可以认为是同方差。第四幅图,是观察杠杆值
点的,我们仅能看出该数据存在杠杆值点,下面我们将对模型进一步假设检验,
并找出相应的异常值。
图 6 三模型直方图
在图 6 中我们看到三个模型的残差拟合效果都不是很好,但是二三幅图中可
以看出在原一元线性模型中加上高次项后核密度曲线的峰度明显降低。说明残差
的正态效果稍有改善,但是还明显达不到正态的效果。
> test(fit_1)
正态检验
Shapiro-Wilk normality test
data: residuals(fit)
W = 0.51439, p-value < 2.2e-16
独立性检验
lag Autocorrelation D-W Statistic p-value
1 -0.0976107 2.195219 0.026
Alternative hypothesis: rho != 0
同方差
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.05713015, Df = 1, p = 0.81109
> test(fit_2)
正态检验
Shapiro-Wilk normality test
data: residuals(fit)
W = 0.58661, p-value < 2.2e-16
独立性检验
lag Autocorrelation D-W Statistic p-value
1 -0.0956992 2.191076 0.024
Alternative hypothesis: rho != 0
同方差
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1.286504, Df = 1, p = 0.25669
> test(fit_3)
正态检验
Shapiro-Wilk normality test
data: residuals(fit)
W = 0.60649, p-value < 2.2e-16
独立性检验
lag Autocorrelation D-W Statistic p-value
1 -0.09804999 2.195813 0.05
Alternative hypothesis: rho != 0
同方差
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1.369829, Df = 1, p = 0.24184
上述分别对残差进行了正态性、独立性、同方差性的检验。我们看到通过上
面的检验,三个模型的正态性假设,均未通过。同方差性通过,独立性勉强不是
很差。下面针对这几个缺点对模型一进行探索,进而改善模型。
先考察数据的异常值:
离群值:
rstudent unadjusted p-value Bonferonni p
21 3.284062 0.0011029 0.50071
强影响点:
StudRes Hat CookD
1 0.0234525 0.008205138 2.280204e-06
6 0.0234525 0.008205138 2.280204e-06
12 3.1897827 0.005859903 2.939057e-02
21 3.2840622 0.002292397 1.212768e-02
28 3.2840622 0.002292397 1.212768e-02
130 3.1897827 0.005859903 2.939057e-02
图 7 异常值观测
从上图我们可以看到 12、130、248 三个点是异常值。
为了使模型更加的满足残差的假设,我们使用函数对求因变量的变换:
> summary(powerTransform(Data$y))
bcPower Transformation to Normality
Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
Data$y 0.7226 0.72 0.5639 0.8813
> boxTidwell(y ~ x13, data = Data)
MLE of lambda Score Statistic (z) Pr(>|z|)
1.0575 2.5643 0.01034 *
iterations = 2
可以看出推荐的变换参数为 0.7226,对模型修改后的到:
我们删除 1、6、12、21、28、130、248 后然后对 y 做 0.7226 的变换。x13
推荐的是 1.0575 接近 1,故不对其进行变换。
新模型信息:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.758691 0.017567 43.19 <2e-16 ***
x13 0.461781 0.002814 164.08 <2e-16 ***
Residual standard error: 0.1636 on 445 degrees of freedom
Multiple R-squared: 0.9837, Adjusted R-squared: 0.9837
F-statistic: 2.692e+04 on 1 and 445 DF, p-value: < 2.2e-16
图 8 新模型的正态性检验图
在改善的模型中 2R =0.9837 ,残差也没有增大,回归方程依然显著,但是常数项
系数通过了改善。在图 8 中我们也看到和密度曲线的峰度也明显减低。