• Posts tagged "T检验"

Blog Archives

用R语言解读统计检验-T检验

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-test-t

前言

做统计分析R语言是最好的,R语言对于统计检验有非常好的支持。我会分7篇文章,分别介绍用R语言进行统计量和统计检验的计算过程,包括T检验F检验卡方检验P值KS检验AICBIC等常用的统计检验方法和统计量。

本文先从T检验开始说起!

目录

  1. T检验介绍
  2. 数据集
  3. 单总体T检验
  4. 双总体T检验

1. T检验介绍

T检验,也称 student t检验(Student’s t test),是戈斯特为了观测酿酒质量发明的。戈斯特于1908年在Biometrika上公布T检验,但因其老板认为其为商业机密而被迫使用笔名(学生)。

T检验,是用于检验两个小样本的平均值差异程度的检验方法,样本量在3-30个左右,要求样本为正态分布,但总体标准差可未知。T检验,是用T分布理论来推断差异发生的概率,从而判定两个样本平均值的差异是否显著。T检验可分为单总体T检验,双总体非配对T检验,和双总体配对T检验。下面将分别进行介绍。

2. 数据集

T检验,对于数据有比较严格的要求,所以我们需要先找到一个合适的数据集,作为测试数据集。我发现了R语言自带的一个数据集sleep,是很好的测试数据集,本文接下来的内容,将以这个数据集进行测试,来介绍T检验。

开发环境所使用的系统环境

  • Win10 64bit
  • R: 3.4.2 x86_64-w64-mingw32/x64 b4bit

数据集,记录了10名患者使用两种催眠药物影响的对比数据,共3列,20条记录。

  • extra列,为睡眠时间的变化
  • group列,为药物编号
  • ID列,为患者编号

查看数据集。


> data(sleep)
> sleep
   extra group ID
1    0.7     1  1
2   -1.6     1  2
3   -0.2     1  3
4   -1.2     1  4
5   -0.1     1  5
6    3.4     1  6
7    3.7     1  7
8    0.8     1  8
9    0.0     1  9
10   2.0     1 10
11   1.9     2  1
12   0.8     2  2
13   1.1     2  3
14   0.1     2  4
15  -0.1     2  5
16   4.4     2  6
17   5.5     2  7
18   1.6     2  8
19   4.6     2  9
20   3.4     2 10

由于T检验的前提条件,总体要符合正态分布,只是总体方差未知。所以,我们需要先对选定数据集进行进行正态分布检验。使用Shapiro-Wilk和qq图,作为正态分布检验的方法。

Shapiro-Wilk正态分布检验: 用来检验是否数据符合正态分布,类似于线性回归的方法一样,是检验其于回归曲线的残差。该方法推荐在样本量很小的时候使用,样本在3到5000之间。该检验原假设为H0:数据集符合正态分布,统计量为W。统计量W最大值是1,越接近1,表示样本与正态分布匹配。p值,当p-value小于显著性水平α(0.05),则拒绝H0。


> shapiro.test(sleep$extra) # 总体正态分布

	Shapiro-Wilk normality test

data:  sleep$extra
W = 0.94607, p-value = 0.3114

检验结果,W接近1,p-value>0.05,不能拒绝原假设,所以数据集sleep是符合正态分布的。

接下来,我们再画出qq图,直观的看一下数据符合正态分布的情况。

图中,对角线基本能穿过数据点,也说明数据符合正态分布。

3. 单总体的T检验

目的:单总体T检验,是检验一个样本平均值与一个已知的总体平均值的差异是否显著。当总体分布是正态分布,如总体标准差未知且样本容量较小,那么样本平均数与总体平均数的离差统计量呈t分布。

适用条件:

  1. 已知总体均值Um
  2. 可得到样本均值Xm,及该样本标准误差Xs
  3. 样本来自正态或近似正态总体

计算公式:T统计量

公式解释:

  • Xm:样本均值
  • Um:总体均值
  • Xs:样本标准差
  • n:样本个数
  • 自由度:df=n-1

用R语言进行T检验,取sleep数据集做为总体,符合正态分布,总体平均值为1.54,以group=1为样本组,共10条记录,样本均值为0.75,样本标准差为1.78901。判断group=1分组的药物,是否有显著性改进睡眠?


#  总体均值
> Um<-mean(sleep$extra) 

# 取group=1的样本组
> g1<-sleep$extra[which(sleep$group==1)]

# 计算样本均值,方差,数量
> Xm<-mean(g1)
> Xs<-sd(g1)
> n<-length(g1)

# 进行T检验
> t.test(g1,mu=Um)

	One Sample t-test

data:  g1
t = -1.3964, df = 9, p-value = 0.1961
alternative hypothesis: true mean is not equal to 1.54
95 percent confidence interval:
 -0.5297804  2.0297804
sample estimates:
mean of x 
     0.75 

指标解释:

  • H0:原假设总体和样本均值无显著差异
  • t统计量:-1.3964
  • df,自由度,10-1=9
  • p-value值:0.1961
  • 95 percent confidence interval:95%的置信区间
  • mean:样本均值0.75

结果解读,以0.05为显著性水平,t=-1.3964小于临界值2.228(查表),t值不显著,不能拒绝原假设。以0.05为显著性水平,p-value=0.19大于0.05,不能拒绝原假设,group=1的药物对睡眠没有明显改进。

我们可以用t值进行显著性差异判断,也可以用p值进行显著性差异判断,他们的作用是一样的。t值判断时,需要用计算所得的t值,与显著性水平查表对比。p值相当于是把t值,进行一种标准化的变型,只和已经定义好的显著性水平比就行了,比如0.05, 0.01, 0.001等几个固定值。

手动计算T值和P值,关于P值的详细解释,请查看文章R语言实现统计检验-P值


# 手动计算T值
> t_stat<-(Xm-Um) / (Xs/(sqrt(n)));t_stat 
[1] -1.396415

# 手动计算P值
> p_value<-2*pt(-abs(t_stat),df=n-1);p_value    
[1] 0.1960699

4. 双总体T检验

双总体T检验,是检验两个样本平均值,与其各自所代表的总体的差异是否显著。双总体T检验又分为两种情况,一种是配对的样本T检验,用于检验两种同质对象,在不同条件下所产生的数据差异性;另一种是独立样本非配对T检验,用于检验两组独立的样本的平均数差异性。

4.1 配对T检验

目标:检验两组同质样本,在不同的处理下的样本平均值,是否有显著的差异性。

配对设计:将2组样本的某些重要特征按相近的原则配成对子,消除混杂因素的影响,观察样本之间的处理因素和研究因素的差异,其它因素基本相同,把配对两组样本个体随机处理。

配对过程如下3种情况:

  • 两种同质样本,分别接受两种不同的处理,如性别、年龄、体重、病情程度相同进行配对。
  • 同一样本或同一样本的两个部分,分别接受两种不同的处理。
  • 同一样本自身对比,把同一组样本处理前后的结果进行比较。

计算公式:

公式解释:

  • D:两个样本差值
  • εD:求和
  • ε(D^2):平方和
  • n:样本个数
  • 自由度:df=n-1

使用sleep数据集,按group分成2个组,形成配对的数据集。H0原假设,g1和g2两个样本组均值没有显著性差异,对治疗睡眠的效果是一致的。


# 配对T检验
> t.test(extra ~ group, data = sleep, paired = TRUE)

	Paired t-test

data:  extra by group
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean of the differences 
                  -1.58 

结果解读,以0.05为显著性水平,p-value=0.002833<0.05,说明g1和g2有显著性的差异,可以拒绝原假设。

接下来,我们动手自己算一下T值


# 生成数据集
> g1<-sleep$extra[which(sleep$group==1)]
> g2<-sleep$extra[which(sleep$group==2)]

# 配对T检验,与上面的计算结果是一致的
> t.test(g1,g2,paired=TRUE)

	Paired t-test

data:  g1 and g2
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean of the differences 
                  -1.58 

# 动手计算
> n1<-length(g1)
> n2<-length(g2)
> D<-g1-g2
> Ds<-sum(D)

# 计算T值
> Ds/sqrt((n1*sum(D^2)-sum(D)^2)/(n1-1))
[1] -4.062128

4.2 非配对T检验

目标:用于检验两组独立的样本的平均数差异性。

计算公式:

  • n1 和n2 为两样本容量
  • Xm1和Xm2为两样本均值
  • Xs1和Xs2为两样本的标准差

继续使用sleep数据集,按group分成2个组,形成独立的样本,病人按随机进行配对,去掉关联关系。H0原假设,随机选取g1和g2两个样本组均值没有显著性差异,对治疗睡眠的效果是一致的。


> t.test(extra ~ group, data = sleep)

	Welch Two Sample t-test

data:  extra by group
t = -1.8608, df = 17.776, p-value = 0.07939
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.3654832  0.2054832
sample estimates:
mean in group 1 mean in group 2 
           0.75            2.33 

结果解读,以0.05为显著性水平,p-value=0.07939>0.05,说明g1和g2有没显著性的差异,不能拒绝原假设,两种方法对于治疗效果是一致的。

接下来,我们动手自己算一下T值。


# 计算数量,均值,标准差
> n1<-length(g1)
> n2<-length(g2)
> Xm1<-mean(g1)
> Xm2<-mean(g2)
> Xs1<-sd(g1)
> Xs2<-sd(g2)
 
# 计算T值
> (Xm1-Xm2)/sqrt((Xs1^2/n1)+(Xs2^2/n2))
[1] -1.860813

画出箱线图,观察两组数据的均值。


> boxplot(extra~group,data=sleep)


图中,每个箱的最粗的黑色线,表示两个样本的均值。

以双总体进行T检验,当是样本是配对时,g1和g2有显著性差异;而当样本是独立时,g1和g2没有显著性差异。由于我们变化了初始的假设,是会很大的程度影响统计结果的,所以在使用统计学模型时,要做非常严格的条件判断,从而保证结果的可靠性,可解释性。

转载请注明出处:
http://blog.fens.me/r-test-t

打赏作者

R语言解读一元线性回归模型

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员R,Nodejs,Java
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-linear-regression/

reg-liner

前言

在我们的日常生活中,存在大量的具有相关性的事件,比如大气压和海拔高度,海拔越高大气压强越小;人的身高和体重,普遍来看越高的人体重也越重。还有一些可能存在相关性的事件,比如知识水平越高的人,收入水平越高;市场化的国家经济越好,则货币越强势,反而全球经济危机,黄金等避险资产越走强。

如果我们要研究这些事件,找到不同变量之间的关系,我们就会用到回归分析。一元线性回归分析是处理两个变量之间关系的最简单模型,是两个变量之间的线性相关关系。让我们一起发现生活中的规律吧。

由于本文为非统计的专业文章,所以当出现与教课书不符的描述,请以教课书为准。本文力求用简化的语言,来介绍一元线性回归的知识,同时配合R语言的实现。

目录

  1. 一元线性回归介绍
  2. 数据集和数学模型
  3. 回归参数估计
  4. 回归方程的显著性检验
  5. 残差分析和异常点检测
  6. 模型预测

1. 一元线性回归介绍

回归分析(Regression Analysis)是用来确定2个或2个以上变量间关系的一种统计分析方法。如果回归分析中,只包括一个自变量X和一个因变量Y时,且它们的关系是线性的,那么这种回归分析称为一元线性回归分析。

回归分析属于统计学的基本模型,涉及统计学基础,就会有一大堆的名词和知识点需要介绍。

在回归分析中,变量有2类:因变量 和 自变量。因变量通常是指实际问题中所关心的指标,用Y表示。而自变量是影响因变量取值的一个变量,用X表示,如果有多个自变量则表示为X1, X2, …, Xn。

回归分析研究的主要步骤:

  1. 确定因变量Y 与 自变量X1, X2, …, Xn 之间的定量关系表达式,即回归方程。
  2. 对回归方程的置信度检查。
  3. 判断自变量Xn(n=1,2,…,m)对因变量的影响。
  4. 利用回归方程进行预测。

本文会根据回归分析的的主要步骤,进行结构梳理,介绍一元线性回归模型的使用方法。

reg

2. 数据集和数学模型

先让我们通过一个例子开始吧,用一组简单的数据来说明一元线性回归分析的数学模型的原理和公式。找出下面数据集中Y与X的定量关系。

数据集为2016年3月1日,白天开盘的交易数据,为锌的2个期货合约的分钟线的价格数据。数据集包括有3列,索引列为时间,zn1.Close为ZN1604合约的1分钟线的报价数据,zn2.Close为ZN1605合约的1分钟线的报价数据。

数据集如下:


                    zn1.Close zn2.Close
2016-03-01 09:01:00     14075     14145
2016-03-01 09:02:00     14095     14160
2016-03-01 09:03:00     14095     14160
2016-03-01 09:04:00     14095     14165
2016-03-01 09:05:00     14120     14190
2016-03-01 09:06:00     14115     14180
2016-03-01 09:07:00     14110     14170
2016-03-01 09:08:00     14110     14175
2016-03-01 09:09:00     14105     14170
2016-03-01 09:10:00     14105     14170
2016-03-01 09:11:00     14120     14180
2016-03-01 09:12:00     14105     14170
2016-03-01 09:13:00     14105     14170
2016-03-01 09:14:00     14110     14175
2016-03-01 09:15:00     14105     14175
2016-03-01 09:16:00     14120     14185
2016-03-01 09:17:00     14125     14190
2016-03-01 09:18:00     14115     14185
2016-03-01 09:19:00     14135     14195
2016-03-01 09:20:00     14125     14190
2016-03-01 09:21:00     14135     14205
2016-03-01 09:22:00     14140     14210
2016-03-01 09:23:00     14140     14200
2016-03-01 09:24:00     14135     14205
2016-03-01 09:25:00     14140     14205
2016-03-01 09:26:00     14135     14205
2016-03-01 09:27:00     14130     14205

我们以zn1.Close列的价格为X,zn2.Close列的价格为Y,那么试试找到自变量X和因变量Y的关系的表达式。

为了直观起见,我们可以先画出一张散点图,以X为横坐标,Y为纵坐标,每个点对应一个X和一个Y。


# 数据集已存在df变量中
> head(df)
                    zn1.Close zn2.Close
2016-03-01 09:01:00     14075     14145
2016-03-01 09:02:00     14095     14160
2016-03-01 09:03:00     14095     14160
2016-03-01 09:04:00     14095     14165
2016-03-01 09:05:00     14120     14190
2016-03-01 09:06:00     14115     14180

# 分别给x,y赋值
> x<-as.numeric(df[,1])
> y<-as.numeric(df[,2])

# 画图
> plot(y~x+1)

01

从散点图上发现 X和Y 的排列基本是在一条直线附近,那么我们可以假设X和Y的关系是线性,可以用公式表式为。

Y = a + b * X + c
  • Y,为因变量
  • X,为自变量
  • a,为截距
  • b,为自变量系数
  • a+b*X, 表示Y随X的变化而线性变化的部分
  • c, 为残差或随机误差,是其他一切不确定因素影响的总和,其值不可观测。假定c是符合均值为0方差为σ^2的正态分布 ,记作c~N(0,σ^2)

对于上面的公式,称函数f(X) = a + b * X 为一元线性回归函数,a为回归常数,b为回归系数,统称回归参数。X 为回归自变量或回归因子,Y 为回归因变量或响应变量。如果(X1,Y1),(X2,Y2)…(Xn,Yn)是(X,Y)的一组观测值,则一元线性回归模型可表示为


Yi = a + b * X + ci,     i= 1,2,...n

其中E(ci)=0, var(ci)=σ^2, i=1,2,...n

通过对一元线性回归模型的数学定义,接下来让我们利用数据集做回归模型的参数估计。

3. 回归参数估计

对于上面的公式,回归参数a,b是我们不知道的,我们需要用参数估计的方法来计算出a,b的值,而从得到数据集的X和Y的定量关系。我们的目标是要计算出一条直线,使直接线上每个点的Y值和实际数据的Y值之差的平方和最小,即(Y1实际-Y1预测)^2+(Y2实际-Y2预测)^2+ …… +(Yn实际-Yn预测)^2 的值最小。参数估计时,我们只考虑Y随X的线性变化的部分,而残差c是不可观测的,参数估计法并不需要考虑残差,对于残差的分析在后文中介绍。

令公式变形为a和b的函数Q(a,b), 即 (Y实际-Y测试)的平方和,变成到(Y实际 – (a+b*X))的平方和。

reg2

公式一 回归参数变形公式

通过最小二乘估计推导出a和b的求解公式,详细的推导过程请参考文章一元线性回归的细节

reg3

公式二 回归参数计算公式

其中 x和y的均值,计算方法如下

reg4

公式三 均值计算公式

有了这个公式,我们就可以求出a和b两个的回归参数的解了。

接下来,我们用R语言来实现对上面数据的回归模型的参数估计,R语言中可以用lm()函数来实现一元线性回归的建模过程。


# 建立线性回归模型
> lm.ab<-lm(y ~ 1+x)

# 打印参数估计的结果
> lm.ab

Call:
lm(formula = y ~ 1 + x)

Coefficients:
(Intercept)            x  
   -349.493        1.029  

如果你想动手来计算也可以自己实现公式。


# x均值
> Xm<-mean(x);Xm 
[1] 14034.82

# y均值
> Ym<-mean(y);Ym
[1] 14096.76

# 计算回归系数
> b <- sum((x-Xm)*(y-Ym)) / sum((x-Xm)^2) ;b
[1] 1.029315

# 计算回归常数
> a <- Ym - b * Xm;a
[1] -349.493

回归参数a和b的计算结果,与lm()函数的计算结果是相同的。有了a和b的值,我们就可以画出这条近似的直接线。

计算公式为:

Y= a + b * X = -349.493 + 1.029315 * X 

画出回归线。


> plot(y~x+1)
> abline(lm.ab)

02

这条直线是我们用数据拟合出来的,是一个近似的值。我们看到有些点在线上,有些点不在线上。那么要评价这条回归线拟合的好坏,我们就需要对回归模型进行显著性检验。

4. 回归方程的显著性检验

从回归参数的公式二可知,在计算过程中并不一定要知道Y和X是否有线性相关的关系。如果不存相关关系,那么回归方程就没有任何意义了,如果Y和X是有相关关系的,即Y会随着X的变化而线性变化,这个时候一元线性回归方程才有意义。所以,我们需要用假设检验的方法,来验证相关性的有效性。

通常会采用三种显著性检验的方法。

  • T检验法:T检验是检验模型某个自变量Xi对于Y的显著性,通常用P-value判断显著性,小于0.01更小时说明这个自变量Xi与Y相关关系显著。
  • F检验法:F检验用于对所有的自变量X在整体上看对于Y的线性显著性,也是用P-value判断显著性,小于0.01更小时说明整体上自变量与Y相关关系显著。
  • R^2(R平方)相关系统检验法:用来判断回归方程的拟合程度,R^2的取值在0,1之间,越接近1说明拟合程度越好。

在R语言中,上面列出的三种检验的方法都已被实现,我们只需要把结果解读。上文中,我们已经通过lm()函数构建一元线性回归模型,然后可以summary()函数来提取模型的计算结果。


> summary(lm.ab)      # 计算结果

Call:
lm(formula = y ~ 1 + x)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.9385  -2.2317  -0.1797   3.3546  10.2766 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.495e+02  7.173e+01  -4.872 2.09e-06 ***
x            1.029e+00  5.111e-03 201.390  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.232 on 223 degrees of freedom
Multiple R-squared:  0.9945,	Adjusted R-squared:  0.9945 
F-statistic: 4.056e+04 on 1 and 223 DF,  p-value: < 2.2e-16

模型解读:

  • Call,列出了回归模型的公式。
  • Residuals,列出了残差的最小值点,1/4分位点,中位数点,3/4分位点,最大值点。
  • Coefficients,表示参数估计的计算结果。
  • Estimate,为参数估计列。Intercept行表示常数参数a的估计值 ,x行表示自变量x的参数b的估计值。
  • Std. Error,为参数的标准差,sd(a), sd(b)
  • t value,为t值,为T检验的值
  • Pr(>|t|) ,表示P-value值,用于T检验判定,匹配显著性标记
  • 显著性标记,***为非常显著,**为高度显著, **为显著,·为不太显著,没有记号为不显著。
  • Residual standard error,表示残差的标准差,自由度为n-2。
  • Multiple R-squared,为相关系数R^2的检验,越接近1则越显著。
  • Adjusted R-squared,为相关系数的修正系数,解决多元回归自变量越多,判定系数R^2越大的问题。
  • F-statistic,表示F统计量,自由度为(1,n-2),p-value:用于F检验判定,匹配显著性标记。

通过查看模型的结果数据,我们可以发现通过T检验的截距和自变量x都是非常显著,通过F检验判断出整个模型的自变量是非常显著,同时R^2的相关系数检验可以判断自变量和因变量是高度相关的。

最后,我们通过的回归参数的检验与回归方程的检验,得到最后一元线性回归方程为:

Y = -349.493 + 1.029315 * X 

5. 残差分析和异常点检测

在得到的回归模型进行显著性检验后,还要在做残差分析(预测值和实际值之间的差),检验模型的正确性,残差必须服从正态分布N(0,σ^2)。

我们可以自己计算数据残差,并进行正态分布检验。


# 残差
> y.res<-residuals(lm.ab)

# 打印前6条数据
> head(y.res)
        1         2         3         4         5         6 
6.8888680 1.3025744 1.3025744 6.3025744 5.5697074 0.7162808 

# 正态分布检验
> shapiro.test(y.res)

	Shapiro-Wilk normality test

data:  y.res
W = 0.98987, p-value = 0.1164

# 画出残差散点图
> plot(y.res)

09

对残差进行Shapiro-Wilk正态分布检验,W接近1,p-value>0.05,证明数据集符合正态分布!关于正态分布的介绍,请参考文章常用连续型分布介绍及R语言实现

同时,我们也可以用R语言的工具生成4种用于模型诊断的图形,简化自己写代码计算的操作。


# 画图,回车展示下一张
> plot(lm.ab)    
Hit  to see next plot:   # 残差拟合图
Hit  to see next plot:   # 残差QQ图
Hit  to see next plot:   # 标准化的残差对拟合值 
Hit  to see next plot:   # 标准化残差对杠杆值

04

图1,残差和拟合值对比图

对残差和拟合值作图,横坐标是拟合值,纵坐标是残差。残差和拟合值之间,数据点均匀分布在y=0两侧,呈现出随机的分布,红色线呈现出一条平稳的曲线并没有明显的形状特征,说明残差数据表现非常好。

05

图2,残差QQ图

残差QQ图,用来描述残差是否符合正态分布。图中的数据点按对角直线排列,趋于一条直线,并被对角直接穿过,直观上符合正态分布。对于近似服从正态分布的标准化残差,应该有 95% 的样本点落在 [-2,2] 区间内。

06

图3,标准化残差平方根和拟合值对比图

对标准化残差平方根和拟合值作图,横坐标是拟合值,纵坐标是标准化后的残差平方根。与残差和拟合值对比图(图1)的判断方法类似,数据随机分布,红色线呈现出一条平稳的曲线,无明显的形状特征。

07

图4,标准残差和杠杆值对比图

对标准化残差和杠杆值作图,虚线表示的cooks距离等高线,通常用Cook距离度量的回归影响点。本图中没有出现红色的等高线,则说明数据中没有特别影响回归结果的异常点。

如果想把把4张图画在一起进行展示,可以改变画布布局。


> par(mfrow=c(2,2))
> plot(lm.ab)

08

看到上面4幅中,每幅图上都有一些点被特别的标记出来了,这些点是可能存在的异常值点,如果要对模型进行优化,我们可以从这些来入手。但终于本次残差分析的结果已经很好了,所以对于异常点的优化,可能并不能明显的提升模型的效果。

从图中发现,索引编号为27和192的2个点在多幅图中出现。我们假设这2个点为异常点,从数据中去掉这2个点,再进行显著性检验和残差分析。


# 查看27和192
> df[c(27,192),]
                    zn1.Close zn2.Close
2016-03-01 09:27:00     14130     14205
2016-03-01 14:27:00     14035     14085

# 新建数据集,去掉27和192
> df2<-df[-c(27,192),]

回归建模和显著性检验。


> x2<-as.numeric(df2[,1])
> y2<-as.numeric(df2[,2])
> lm.ab2<-lm(y2 ~ 1+x2)
> summary(lm.ab2)

Call:
lm(formula = y2 ~ 1 + x2)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0356 -2.1542 -0.2727  3.3336  9.5879 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.293e+02  7.024e+01  -4.688 4.83e-06 ***
x2           1.028e+00  5.004e-03 205.391  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.117 on 221 degrees of freedom
Multiple R-squared:  0.9948,	Adjusted R-squared:  0.9948 
F-statistic: 4.219e+04 on 1 and 221 DF,  p-value: < 2.2e-16

对比这次的显著性检验结果和之前结果,T检验,F检验 和 R^2检验,并没有明显的效果提升,结果和我预想的是一样的。所以,通过残差分析和异常点分析,我认为模型是有效的。

6. 模型预测

最后,我们获得了一元线性回归方程的公式,就可以对数据进行预测了。比如,对给定X=x0时,计算出y0=a+b*x0的值,并计算出置信度为1-α的预测区间。

当X=x0,Y=y0时,置信度为1-α的预测区间为

reg5

reg6

我们可以用R语言的predict()函数来计算预测值y0,和相应的预测区间。程序算法如下。


> new<-data.frame(x=14040)
> lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95)

# 预测结果
> lm.pred
       fit      lwr      upr
1 14102.09 14093.73 14110.44

当x0=14040时,在预测区间为0.95的概率时,y0的值为 14102,预测区间为[14093.73,14110.44]。

我们通过图形来表示。


> plot(y~x+1)
> abline(lm.ab,col='red')
> points(rep(newX$x,3),y=lm.pred,pch=19,col=c('red','blue','green'))

03

其中,红色点为y0的值,蓝色点为预测区间最小值,绿色点为预测区间最大值。

对于统计模型中最核心部分就在结果解读,本文介绍了一元回归模型的基本的建模过程和模型的详细解读方法。在我们掌握了这种方法以后,就可以更容易地理解和学习 多元回归,非线性回归 等更多的模型,并把这些模型应用到实际的工作中了。下一篇文章R语言解读多元线性回归模型,请继续阅读。

转载请注明出处:
http://blog.fens.me/r-linear-regression/

打赏作者