《第四课回归分析.doc》由会员分享,可在线阅读,更多相关《第四课回归分析.doc(17页珍藏版)》请在三一文库上搜索。
1、回归分析一、回归分析:是建立因变量Y与自变量X之间关系的模型,一元线性回归使用一个自变量X,多元线性回归使用超过一个自变量,如:.1 线性回归模型 ,进行n次独立观测:, i=1,2,n称为残差向量2 回归系数显著性检验:判断某一自变量的系数是否为0,如果统计量(表示)的值(通常取0.05),认为。3 回归方程显著性检验:检查是否可用线性方程来处理数据。即检验方程系数是否全为0,如果统计量(表示)的值,认为可以用线性方程来处理问题。4 相关性检验:检验因变量与自变量相关性程度。用相关系数的平方来衡量,的值接近1表示相关,接近0表示不相关。二、线性回归模型的计算 1、 lm()函数:完成多元线性
2、回归系数的估计、回归系数的检验、回归方程的检验等工作,返回值称为拟合结果的对象,本质上是一个具有类属性值lm的列表.格式:lm(formula,data,subset,weights,na.action,method=”qr”,model=TRUE,x=FALSE,y=FALSE,qr=TRUE,singular.ok=TRUE,contrasts=NULL,offset,) formula 为模型公式,形如y1+x1+x2的形式,表示常数项、的系数和的系数,去掉 公式中的1,其意义不变。拟合成齐次线性模型:y0+x1+x2或y -1+x1+x2的形式;data 为数据框,由样本数据构成;su
3、bset 为可选项,表示所使用的样本子集;weights为可选向量,表示对应样本的权重;na.action 为函数,表示当数据中出现缺失数据(NA)的处理方法;method 为估计回归系数的计算方法,默认值为“qr”; model、x、y、qr为逻辑变量,如果取TRUE,函数的返回值将给出模型的框架、模型矩阵、响应变量,及QR分解; singular.ok 为逻辑变量,取FALSE表示奇异值拟合是错误的; contrasts 为可选列表;offset 为NULL,或者是数值向量。为附加参数。 2、 summary()函数 lm()函数的返回值称为拟合结果的对象,是一个有类属性值lm的列表,有m
4、odel, coefficients, residuals(残差)等成员,为了获得更多信息,summary()常与lm( )函数一起使用,格式:summary(object, correlation=FALSE, symbolic.cor=FALSE,) 参数意义:object为lm()函数生成的对象,correlation为逻辑变量,取TRUE表示给出估计参数的相关矩阵,symbolic.cor为逻辑变量,取TRUE表示用符号形式给出估计参数的相关矩阵,此参数只有当correlation=TRUE时才有效。例1 年龄相等情况下,建立血压的收缩压与体重(kg),年龄(岁数)有关,建立与,的线性
5、回归方程。解blood-data.frame( X1=c(76.0, 91.5, 85.5, 82.5, 79.0, 80.5, 74.5, 79.0, 85.0, 76.5, 82.0, 95.0, 92.5), X2=c(50, 20, 20, 30, 30, 50, 60, 50, 40, 55, 40, 40, 20), Y= c(120, 141, 124, 126, 117, 125, 123, 125,132, 123, 132, 155, 147)lm.sol|t|) (Intercept) -62.96336 16.99976 -3.704 0.004083 * X1 2.1
6、3656 0.17534 12.185 2.53e-07 *X2 0.40022 0.08321 4.810 0.000713 *- 极为显著*,高度显著*,显著*Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 2.854 on 10 degrees of freedom #残差的标准差Multiple R-squared: 0.9461, Adjusted R-squared: 0.9354 #相关系数的平方,修正相关系数的平方F-statistic: 87.84 on 2 and 10 DF, p
7、-value: 4.531e-07 #F统计量,F统计量对应的P值回归方程与回归系数的检验都是显著的,回归方程为: 3 预测区间与置信区间 如果回归效果显著,就可以利用回归方程进行预测,即对给定的回归自变量值,预测因变 量可能的取值范围,是一个区间估计问题。 假设=, 回归方程的真实值(未知):相应的估计值:predict( )函数计算的估计值、预测区间和置信区间。格式:predict(object,newdata,se.fit=FALSE,scale=NULL, df=Inf,interval=c(“none”,”confidence”,”prediction”),level=0.95,ty
8、pe=c(“response”,”terms”),terms=NULL,na.action=na.pass,pred.var=res.var/weights,weights=1,) 参数意义: object 为lm()函数得到的对象; newdata 为数据框,由预测点构成,如取默认值,将计算已知数据的回归值; se.fit 为逻辑变量,取TRUE表示输出预测值的标准差、自由度和残差尺度信息; scale 为计算标准差的尺度参数; df 为尺度参数的自由度; interval 为计算的区间类型,取“none”(默认值)表示不计算,取”confidence”表示计算置信区间,取”predicti
9、on”表示计算预测区间; level 为置信水平,默认值为0.95,在计算置信区间和预测区间时用到; type 为预测类型,或者取“response”(默认值), 或者取”terms”,当type=”terms”时,为全部选择项; terms 为选择项,只能取NULL,1,2等; na.action 为函数,表示处理newdata中有缺失数据(NA)时的处理方法,默认预测为NA.pred.var预测区间的方差; weights 为数值向量或单侧模型公式,用于预测时方差的权; 附加参数。例2 继续上面的例1,设,求的估计值、预测区间和置信区间(置信水平为0.95)。解blood-data.fra
10、me( X1=c(76.0, 91.5, 85.5, 82.5, 79.0, 80.5, 74.5, 79.0, 85.0, 76.5, 82.0, 95.0, 92.5), X2=c(50, 20, 20, 30, 30, 50, 60, 50, 40, 55, 40, 40, 20), Y= c(120, 141, 124, 126, 117, 125, 123, 125,132, 123, 132, 155, 147)lm.sol-lm(Y 1+X1+X2, data=blood)newdata - data.frame(X1 = 80, X2 = 40)predict(lm.sol,
11、newdata, interval=prediction) predict(lm.sol, newdata, interval=confidence)例 分析合金强度与合金中含碳量之间的关系:(它们之间的一组数据如下)(1)完成一元线性回归的计算;(2)计算自变量在区间0.10,0.23内的回归方程的预测估计值、预测区间和置信区间(取),并将数据点、预测估计曲线、预测区间曲线和置信区间曲线画在同一张图上。解 # 回归估计x-c(0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.20, 0.21, 0.23)y-c(42.0, 43.
12、5, 45.0, 45.5, 45.0, 47.5, 49.0, 53.0, 50.0, 55.0, 55.0, 60.0)lm.sol-lm(y 1+x)summary(lm.sol)# 计算预测值, 并绘图new - data.frame(x = seq(0.10, 0.24, by=0.01) #取x的范围值pp-predict(lm.sol, new, interval=prediction)pc|t|)(Intercept) 3.0000909 1.1247468 2.667348 0.025734051x1 0.5000909 0.1179055 4.241455 0.002169
13、6292 Estimate Std. Error t value Pr(|t|)(Intercept) 3.000909 1.1253024 2.666758 0.025758941x2 0.500000 0.1179637 4.238590 0.0021788163 Estimate Std. Error t value Pr(|t|)(Intercept) 3.0024545 1.1244812 2.670080 0.025619109x3 0.4997273 0.1178777 4.239372 0.0021763054 Estimate Std. Error t value Pr(|t
14、|)(Intercept) 3.0017273 1.1239211 2.670763 0.025590425x4 0.4999091 0.1178189 4.243028 0.002164602 并且进一步检验4组数据的、F值检验全合格,但绘出图分析4组数据做线性回归模型不是全合格: 只有第一组数据适合做线性回归模型。1 残差检验 检验模型误差是否满足正态性和方差齐性。1)用残差图检验,残差为纵坐标,拟合值或数据观测序号或数据观测时间做横坐标,做出的散点图称为残差图。 (a) (b) (c) 正常情况 异方差情况 非线性情况 有异常值点的情况检查1)有无有异常值点;2)画出残差的Q-Q散点图,
15、若这些点位于一条直线上,则残差服从正态分布,否则不然。注意:上面图中只有(a)图是正常的,(b)图是喇叭口,回归值的大小与残差的波动有关,方差的假设有问题,(c)图表明线性模型不合适。residuals()(或resid()函数)计算回归模型的残差,rstandard()计算回归模型的标准化残差,rstudent()计算回归模型的化学生化残差(删除第i个样本数据后得到的标准化残差)格式:residuals(object, .)resid(object, .)rstandard(model, infl = lm.influence(model, do.coef = FALSE), sd = sq
16、rt(deviance(model)/df.residual(model), .)rstudent(model, infl = lm.influence(model, do.coef = FALSE), res = infl$wt.res, .)参数object或model为lm生成的对象;infl为lm.influence返回值得到的影响结构;sd 为模型的标准差;res为模型残差。plot()函数绘出各种残差的散点图,格式:plot(x,which=c(1:3,5),.) ,参数x由lm生成的对象,which是开关变量,其值16,默认子集(1,2,3,5),即绘出1、2、3和5号散点图,该
17、函数可以绘出6张诊断图,第一张残差与预测值散点图,第二张残差的正态Q-Q图,第三张标准差的平方根与预测值的散点图等。例 2060成年女性的血压与年龄之间的回归关系,并作残差分析解# (1) 回归rt-read.table(blood.dat, header=TRUE)lm.sol-lm(YX, data=rt); lm.solsummary(lm.sol)# (2) 残差图pre-fitted.values(lm.sol) #拟合值res-residuals(lm.sol) #残差rst-rstandard(lm.sol) #标准差par(mai=c(0.9, 0.9, 0.2, 0.1)pl
18、ot(pre, res, xlab=Fitted Values, ylab=Residuals)savePlot(resid-1, type=eps)plot(pre, rst, xlab=Fitted Values, ylab=Standardized Residuals)savePlot(resid-2, type=eps)# (3) 对残差作回归rt-read.table(blood.dat, header=TRUE)rt$res-res #给rt数据框增加一列reslm.res-lm(abs(res)X, data=rt); #利用残差的绝对值与自变量做回归lm.ressummary(
19、lm.res)# (4) 修正回归方程s-lm.res$coefficients1+lm.res$coefficients2*rt$X #计算残差的标准差lm.weg-lm(YX, data=rt, weights=1/s2); #用标准差平方的倒数作权重 #权重的选取可减少非齐性方差带来的影响lm.wegsummary(lm.weg) 2 影响分析 影响分析就是检查是否有对估计有异常大数据(高杠杆点)情况。影响分析的方法有:DFFITS准则、Cook距离、COVRATIO准则、帽子值和帽子矩阵。 相关的R函数有:influence.measures(), dffits(), dfbeta()
20、, dfbetas(), cooks.distance(), covratio(), hatvalues(), hat();格式:influence.measures(model)dffits(model, infl = , res = )dfbeta(model, infl = lm.influence(model, do.coef = TRUE), .)dfbetas(model, infl = lm.influence(model, do.coef = TRUE), .)covratio(model, infl = lm.influence(model, do.coef = FALSE)
21、, res = weighted.residuals(model)cooks.distance(model, infl = lm.influence(model, do.coef = FALSE), res = weighted.residuals(model), sd = sqrt(deviance(model)/df.residual(model), hat = infl$hat, .)hatvalues(model, infl = lm.influence(model, do.coef = FALSE), .)hat(x, intercept = TRUE) 参数model为lm生成的对
22、象,infl为lm.influence返回值得到的影响结构,res为模型残差,x为设计矩阵。例 通过智力测试数据建立智力随年龄变化的关系(在exam0606.R文件中)。解 第一步:运行,计算回归系数,并做回归系数与回归方程检验# 1. 回归分析intellect-data.frame( x=c(15, 26, 10, 9, 15, 20, 18, 11, 8, 20, 7, 9, 10, 11, 11, 10, 12, 42, 17, 11, 10), #年龄(月)y=c(95, 71, 83, 91, 102, 87, 93, 100, 104, 94, 113, 96, 83, 84,
23、102, 100, 105, 57, 121, 86, 100) #智力 ) lm.sol-lm(y1+x, data=intellect)summary(lm.sol)第二步 做回归诊断# 2. 回归诊断op - par(mfrow=c(2,2), mar=0.4+c(4,4,1,1), oma= c(0,0,2,0) plot(lm.sol, 1:4) #通过绘图分析有无异常值点par(op)savePlot(diagnoses1, type=eps)influence.measures(lm.sol) #通过函数计算分析有无异常值点,下面是运行结果: dfb.1_ dfb.x dffit
24、 cov.r cook.d hat inf1 0.01664 0.00328 0.04127 1.166 8.97e-04 0.0479 2 0.18862 -0.33480 -0.40252 1.197 8.15e-02 0.1545 . 18 0.83112 -1.11275 -1.15578 2.959 6.78e-01 0.6516 *19 0.14348 0.27317 0.85374 0.396 2.23e-01 0.0531 *20 -0.20761 0.10544 -0.26385 1.043 3.45e-02 0.0567 21 0.02791 -0.01622 0.0329
25、8 1.187 5.74e-04 0.0628 influence.measures()函数得到回归诊断结果共7列,1,2列(dfbetas)对应常数和变量x值, 3列是DFFITS准则值,4列是COVRATIO准则值,5列是Cook距离,6列高杠杆值,7列是影响点记号。(18,19号样本是强影响点记号)。第三步 处理强影响点# 3. 回归修正n-length(intellect$x)weights-rep(1, n); weights18-0.5lm.correct-lm(y1+x, data=intellect, subset=-19, weights=weights)summary(lm
26、.correct)第四步 验证# 4. 验证 #在脚本里做intellect-data.frame( x=c(15, 26, 10, 9, 15, 20, 18, 11, 8, 20, 7, 9, 10, 11, 11, 10, 12, 42, 17, 11, 10), #年龄(月)y=c(95, 71, 83, 91, 102, 87, 93, 100, 104, 94, 113, 96, 83, 84, 102, 100, 105, 57, 121, 86, 100) #智力 ) lm.sol-lm(y1+x, data=intellect)attach(intellect) #把inte
27、llect调入内存,可操作它的元素par(mai=c(0.8, 0.8, 0.2, 0.2)plot(x, y, cex=1.2, pch=21, col=red, bg=orange) #真值abline(lm.sol, col=blue, lwd=2) #修正前,lm.sol是拟合值text(xc(19, 18), yc(19, 18), label=c(19, 18), adj=c(1.5, 0.3) # adj()调整符号坐标位置detach()abline(lm.correct, col=red, lwd=2, lty=5) #绘一条修正后线legend(30, 120, c(Poi
28、nts, Regression, Correct Reg), #点Points、回归Regression、修正回归Correct Reg pch=c(19, NA, NA), lty=c(NA, 1,5), # NA代表缺省 col=c(orange, blue, red)savePlot(diagnoses2, type=eps)# 5. 检验op - par(mfrow=c(2,2), mar=0.4+c(4,4,1,1), oma= c(0,0,2,0) #修改图边空plot(lm.correct, 1:4)par(op) #退出修改参数savePlot(diagnoses3, type
29、=eps) 调整高杠杆点后.完.聚类分析思想:将相似性(或相异性)数据看成是对象之间的“距离”远近的一种度量,将距离近的对象归入一类,不同类之间的对象距离较远。根据分类对象不同分为Q型聚类分析和R型聚类分析,Q型聚类分析是指:对样本进行聚类;R型聚类分析是指:对变量进行聚类。1. 距离和相似系数:变量大致可分为两类1)定量变量,是连续值,具有数值特征,数量上的变化;2)定性变量,这些量并非具有数量上的变化,而只有性质上的差异,大致分两类,一为有序变量,没有数量关系,只有次序关系,如一等品、二等品等;二为名义变量,既无等级,也无次序,如天气(阴,晴)等。1.1. 距离:若数据集中有n个样本,则这
30、n个样本就是数据集中的n个点,如每个样本有k个指标(分量),则为第i个样本的第k个指标,第i个样本的第j个样本的距离记为,聚类过程中,距离近的归为一类,距离远的归属不同类。在R中,dist( )函数计算各样本间的距离,格式: dist(x, method = euclidean, diag = FALSE, upper = FALSE, p = 2)参数x为数值矩阵,或者为数据框,或者是“dist”对象;method为定义距离的方法,默认值是euclidean,表示欧几里德距离(即),还可选绝对值距离、无穷模距离等;diag为逻辑变量,表示输出上三角矩阵值,默认为FALSE,表示仅输出下三角阵
31、值;p表示距离的参数q,默认值为2,即欧几里德距离。1.2 数据变换做聚类分析时,需要将数据做中心化或标准化处理,称:,i=1,2,.,n;j=1,2,.,p为中心化变换,称:,i=1,2,.,n;j=1,2,.,p为标准化变换在R中,scale( )函数做数据的中心化或标准化,格式: scale(x, center = TRUE, scale = TRUE)参数x为样本构成的数值矩阵,center 或者为逻辑变量,表示是否对数据做中心变换,默认值为 TRUE;或者为数值向量,其维数等于矩阵x的列数,表示以center为中心做中心化变换;scale或者为逻辑变量,表示对数据是否做标准变换,默认
32、值为 TRUE,或者为数值向量,其维数等于矩阵x的列数,表示以scale为尺度做标准化变换。 1.3 相似系数 聚类分析不仅对样本分类,而且对变量进行分类,对变量分类时,常用相似系数度量变量之间的相似程度。设表示变量和间的相似系数,要求:(1) =1,当且仅当;(2) ,对一切i,j成立;(3) ,对一切i,j成立。越接近1,则表示和间的关系越密切,越接近0两者间关系越疏远。用cor( )函数计算相关系数。变量间的距离可用相关系数求解: (往往用在变量聚类中)2 系统聚类法最普遍方法,思想:将距离最近类合并成一个类,重复进行两最近类合并,直至所有样本合并为一个类,方法很多:最短距离法:“sin
33、gle”最长距离法:“complete”中间距离法:“median”类平均法:“average”重心法:“centroid”离差平方和法:“Ward”Mcquitty相似法:“mcquitty”2.1 hclust( )函数计算系统聚类,plot( )函数画出系统聚类的树形图(谱系图), plclust( )函数也可画谱系图;格式:hclust(d, method = complete, members = NULL)参数d为dist( )函数生成的对象,即距离。method为系统聚类的方法,默认值为“complete”代表最长距离法,“single”代表最短距离法,“median”代表中间距
34、离法,“average”代表类平均法等;members或者为NULL(默认值),或者为与d有相同变量长度的向量。格式:plot(x, labels=NULL,hang=0.1, .)参数x为hclust( )函数生成的对象;labels为树叶的标记,默认值为NULL; hang为数值,表明谱系图中各类所在的位置,默认值为0.1,取负值表示谱系图中的类从底部画起;其他参数意义与plot( )函数相同。格式:plclust(tree, hang = 0.1, unit = FALSE, level = FALSE, hmin = 0, square = TRUE, labels = NULL, p
35、lot. = TRUE, axes = TRUE, frame.plot = FALSE, ann = TRUE, main = , sub = NULL, xlab = NULL, ylab = Height)参数tree为hclust( )函数生成的对象;unit为逻辑变量,取TRUE表示分叉画在等空间高度,而不是在对象的高度。其他参数意义见帮助文档。例1 设有5个样本,每个样本只有一个指标,分别是1,2,6,8,11,样本间的距离选用欧几里得距离,试用最短距离法、最长距离法等方法进行聚类分析,并画出相应的谱系图。#% 输入数据, 生成距离结构x-c(1,2,6,8,11); dim(x)
36、-c(5,1); d-dist(x)#% 生成系统聚类hc1-hclust(d, single); hc2-hclust(d, complete)hc3-hclust(d, median); hc4-hclust(d, average)#% 绘出所有树形结构图, 并以/2* 2/的形式绘在一张图上opar - par(mfrow = c(2, 2), omi=rep(0,4)plot(hc1,hang=-1); plot(hc2,hang=-1)plot(hc3,hang=-1); plot(hc4,hang=-1)par(opar)#% 保存图形savePlot(hclust_1, type
37、=eps)#% 用plclust函数绘图opar - par(mfrow = c(2, 2), omi=rep(0,4), mar=c(5,4,1,1)plclust(hc1,hang=-1); plclust(hc2,hang=-1) #绘图结果一样plclust(hc3,hang=-1); plclust(hc4,hang=-1)par(opar)3 cophenetic函数 计算系统聚类的cophenetic距离;格式:cophenetic(x), 参数x为hclust( )函数或as.hclust( )函数生成的对象。用途:计算Cophenetic距离与dist( )函数的距离的相关系数,用以评价聚类方法好坏,越接近1,方法越好。例2 上例1中4种聚类方法,分别计算cophenetic距离和各自与距离d的相关系数: #% 使用cophenetic函数, 作聚类方法的比较x-c(1,2,6,8,11); dim(x)-c(5,1); d-dist(x)method=c(single, complete, median, average)cc-numeric(0)for (m in method) dc-cophenetic(hclust(d, m) ccm-cor(d, dc)cc执行结
链接地址:https://www.31doc.com/p-2715696.html