r语言线性回归分析中如何剔除极值?

R语言进行数据挖掘决策树和随机森林
1、使用包party建立决策树
这一节学习使用包party里面的函数ctree()为数据集iris建立一个决策树。属性Sepal.Length(萼片长度)、Sepal.Width(萼片宽度)、Petal.Length(花瓣长度)以及Petal.Width(花瓣宽度)被用来预测鸢尾花的Species(种类)。在这个包里面,函数ctree()建立了一个决策树,predict()预测另外一个数据集。在建立模型之前,iris(鸢尾花)数据集被分为两个子集:训练集(70%)和测试集(30%)。使用随机种子设置固定的随机数,可以使得随机选取的数据是可重复利用的。
观察鸢尾花数据集的结构
设置随机数起点为1234
set.seed(1234)
使用sample函数抽取样本,将数据集中观测值分为两个子集
ind &- sample(2, nrow(iris), replace=TRUE, prob=c(0.7, 0.3))
样本的第一部分为训练集
trainData &- iris[ind==1,]
样本的第二部分为测试集
testData &- iris[ind==2,]加载包party建立一个决策树,并检测预测见过。函数ctree()提供一些参数例如MinSplit, MinBusket, MaxSurrogate 和 MaxDepth用来控制决策树的训练。下面我们将会使用默认的参数设置去建立决策树,至于具体的参数设置可以通过?party查看函数文档。下面的代码中,myFormula公式中的Species(种类)是目标变量,其他变量是独立变量。library(party)
符号'~'是连接方程或公式左右两边的符号
myFormula &- Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
建立决策树
iris_ctree &- ctree(myFormula, data=trainData)
检测预测值
table(predict(iris_ctree), trainData$Species)显示结果如下:
由上图可知,setosa(山鸢尾)40观测值全部正确预测,而versicolor(变色鸢尾)有一个观测值被误判为virginica(维吉尼亚鸢尾),virginica(维吉尼亚鸢尾)有3个观测值被误判为versicolor(变色鸢尾)。
打印决策树
print(iris_ctree)plot(iris_ctree)plot(iris_ctree, type=&simple&)
在图1中,每一个叶子的节点的条形图都显示了观测值落入三个品种的概率。在图2中,这些概率以每个叶子结点中的y值表示。例如:结点2里面的标签是“n=40 y=(1,0,0)”,指的是这一类中一共有40个观测值,并且所有的观测值的类别都属于第一类setosa(山鸢尾)。
接下来,需要使用测试集测试决策树。
在测试集上测试决策树
testPred &- predict(iris_ctree, newdata = testData)table(testPred, testData$Species)结果如下:
从上图的结果可知,决策树对变色鸢尾和维吉尼亚鸢尾的识别仍然有误判。因此ctree()现在的版本并不能很好的处理部分属性不明确的值,在实例中既有可能被判到左子树,有时候也会被判到右子树上。2、使用包rpart建立决策树rpart这个包在本节中被用来在'bodyfat'这个数据集的基础上建立决策树。函数raprt()可以建立一个决策树,并且可以选择最小误差的预测。然后利用该决策树使用predict()预测另外一个数据集。首先,加载bodyfat这个数据集,并查看它的一些属性。
data(&bodyfat&, package = &TH.data&)dim(bodyfat)attributes(bodyfat)bodyfat[1:5,]跟第1节一样,将数据集分为训练集和测试集,并根据训练集建立决策树。set.seed(1234)ind &- sample(2, nrow(bodyfat), replace=TRUE, prob=c(0.7, 0.3))bodyfat.train &- bodyfat[ind==1,]bodyfat.test &- bodyfat[ind==2,]library(rpart)
编写公式myFormula
myFormula &- DEXfat ~ age + waistcirc + hipcirc + elbowbreadth + kneebreadth
训练决策树
bodyfat_rpart &- rpart(myFormula, data = bodyfat.train,
control = rpart.control(minsplit = 10))
plot(bodyfat_rpart)
添加文本标签
text(bodyfat_rpart, use.n=T)结果如下图所示:
选择预测误差最小值的预测树,从而优化模型。
opt &- which.min(bodyfat_rpart$cptable[,&xerror&])cp &- bodyfat_rpart$cptable[opt, &CP&]bodyfat_prune &- prune(bodyfat_rpart, cp = cp)plot(bodyfat_rpart)text(bodyfat_rpart, use.n=T)优化后的决策树如下:
对比结果就会发现,优化模型后,就是将hipcirc&99.5这个分层给去掉了,也许是因为这个分层没有必要,那么大家可以思考一下选择预测误差最小的结果的决策树的分层反而没有那么细。之后,优化后的决策树将会用来预测,预测的结果会与实际的值进行对比。下面的代码中,使用函数abline()绘制一条斜线。一个好的模型的预测值应该是约接近真实值越好,也就是说大部分的点应该落在斜线上面或者在斜线附近。
根据测试集预测
DEXfat_pred &- predict(bodyfat_prune, newdata=bodyfat.test)
预测值的极值
xlim &- range(bodyfat$DEXfat)plot(DEXfat_pred ~ DEXfat, data=bodyfat.test, xlab=&Observed&,
ylab=&Predicted&, ylim=xlim, xlim=xlim)
abline(a=0, b=1)绘制结果如下:
3、随机森林我们使用包randomForest并利用鸢尾花数据建立一个预测模型。包里面的randomForest()函数有两点不足:第一,它不能处理缺失值,使得用户必须在使用该函数之前填补这些缺失值;第二,每个分类属性的最大数量不能超过32个,如果属性超过32个,那么在使用randomForest()之前那些属性必须被转化。也可以通过另外一个包'cforest'建立随机森林,并且这个包里面的函数并不受属性的最大数量约束,尽管如此,高维的分类属性会使得它在建立随机森林的时候消耗大量的内存和时间。
ind &- sample(2, nrow(iris), replace=TRUE, prob=c(0.7, 0.3))trainData &- iris[ind==1,]testData &- iris[ind==2,]library(randomForest)
Species ~ .指的是Species与其他所有属性之间的等式
rf &- randomForest(Species ~ ., data=trainData, ntree=100, proximity=TRUE)table(predict(rf), trainData$Species)结果如下:
由上图的结果可知,即使在决策树中,仍然有误差,第二类和第三类话仍然会被误判,可以通过输入print(rf)知道误判率为2.88%,也可以通过输入plot(rf)绘制每一棵树的误判率的图。最后,在测试集上测试训练集上建立的随机森林,并使用table()和margin()函数检测预测结果。
irisPred &- predict(rf, newdata=testData)table(irisPred, testData$Species)
绘制每一个观测值被判断正确的概率图
plot(margin(rf, testData$Species))显示结果如下:
> 本站内容系网友提交或本网编辑转载,其目的在于传递更多信息,并不代表本网赞同其观点和对其真实性负责。如涉及作品内容、版权和其它问题,请及时与本网联系,我们将在第一时间删除内容!
今天发现一个很不错的博客(), 博主致力于研究R语言在数据挖掘方面的应用,正好近期很想系统的学习一下R语言和数据挖掘的整个流程,看了这个博客的内容,心里久久不能平静.决定从今天 开始,只要晚上能在11点之前把碗洗好,就花一个小时的时间学习博客上的内容,并把学习过程中记不住的信息记录下来,顺便把离英语四级的差 ...
摘要: 今天发现一个很不错的博客(),博主致力于研究R语言在数据挖掘方面的应用,正好近期很想系统的学习一下R语言和数据挖掘的整个流程,看了这个博客的内容,心里久久不能平静.决定从今天开始 ... 今天发现一个很不错的博客(),博主致力于研究R语言在数据挖掘 ...
(1):数据挖掘相关包的介绍 今天发现一个很不错的博客(),博主致力于研究R语言在数据挖掘方面的应用,正好近期很想系统的学习一下R语言和数据挖掘的整个流程,看了这个博客的内容,心里久久不能平静.决定从今天开始,只要晚上能在11点之前把碗洗好,就花一个小时的时间学习博客上的内容,并把学习过程中记不住的信息记录 ...
今天发现一个很不错的博客(),博主致力于研究R语言在数据挖掘方面的应用,正好近期很想系统的学习一下R语言和数据挖掘的整个流程,看了这个博客的内容,心里久久不能平静.决定从今天开始 ... 今天发现一个很不错的博客(), 博主致力于研究R语言在数据挖掘方面的 ...
今天开了第二次的各个比赛小组之间的交流会议,会上我说明了一下我们小组的进度,并得到了老师对我们小组进度的肯定.尽管老师认为我们的速度算快了,但是我觉得我们还是有很多问题没有解决,并没有想象的那么乐观.出现的问题主要由以下几点:
1.对R语言以 ...
第二次进度汇报:
加伟:搭建了服务器客户端环境,正在调试数据包的互传
志杰:摒弃了Rserve
赖威:确定使用JRI进行调用R.
下一个三天的计划:
加伟:实现服务器客户端互传.
志杰:实现一种分类算法在R软件上的实现,在Java中调用的效果
赖威:同上, ...
1.线性回归线性回归就是使用下面的预测函数预测未来观测量: 其中,x1,x2,...,xk都是预测变量(影响预测的因素),y是需要预测的目标变量(被预测变量). 线性回归模型的数据来源于澳大利亚的CPI数据,选取的是2008年到2011年的季度数据. rep函数里面的第一个参数是向量的起始时间,从,第二个参数表示向量里面的每个元素都被4个小 ...
原文来源:/examples/decision-tree 并加上个人整理R语言基础入门之五:简单线性回归 - 生物信息 - 生物秀
标题: R语言基础入门之五:简单线性回归
摘要: [R语言基础入门之五:简单线性回归]线性回归可能是数据分析中最为常用的工具了,如果你认为手上的数据存在着线性定量关系,不妨先画个散点图观察一下,然后用线性回归加以分析。下面简单介绍一下如何在R中进行线性回归。一、回归建模
我们利用R语言中内置的trees数据,其中包含了Volume(体积)、Girth(树围)、Height(树高)这三个变量,我们希望以体积为…… [关键词:回归 变量变换 标准化 散点图 正态分布 统计量 自变量]……
线性回归可能是数据分析中最为常用的工具了,如果你认为手上的数据存在着线性定量关系,不妨先画个散点图观察一下,然后用线性回归加以分析。下面简单介绍一下如何在R中进行线性回归。
一、回归建模
我们利用R语言中内置的trees数据,其中包含了Volume(体积)、Girth(树围)、Height(树高)这三个变量,我们希望以体积为因变量,树围为自变量进行线性回归。
plot(Volume~Girth,data=trees,pch=16,col="red")
model=lm(Volume~Girth,data=trees)
abline(model,lty=2)
summary(model)
首先绘制了两变量的散点图,然后用lm函数建立线性回归模型,并将回归直线加在原图上,最后用summary将模型结果进行了展示,从变量P值和F统计量可得回归模型是显著的。但截距项不应该为负数,所以也可以用下面方法将截距强制为0。
model2=lm(Volume~Girth-1,data=trees)
二、模型诊断
在模型建立后会利用各种方式来检验模型的正确性,对残差进行分析是常见的方法,下面我们来生成四种用于模型诊断的图形。
par(mfrow=c(2,2))
plot(model)
par(mfrow=c(1,1))
这里左上图是残差对拟合值作图,整体呈现出一种先下降后下升的模式,显示残差中可能还存在未提炼出来的影响因素。右上图残差QQ图,用以观察残差是否符合正态分布。左下图是标准化残差对拟合值,用于判断模型残差是否等方差。右下图是标准化残差对杠杆值,虚线表示的cooks距离等高线。我们发现31号样本有较大的影响。
三、变量变换
因为31号样本有着高影响力,为了降低其影响,一种方法就是将变量进行开方变换来改善回归结果,从残差标准误到残差图,各项观察都说明变换是有效的。
plot(sqrt(Volume)~Girth,data=trees,pch=16,col="red")
model2=lm(sqrt(Volume)~Girth,data=trees)
abline(model2,lty=2)
summary(model2)
四、模型预测
下面根据上述模型计算预测值以及置信区间,predict函数可以获得模型的预测值,加入参数可以得到预测区间
plot(sqrt(Volume)~Girth,data=trees,pch=16,col="red")
model2=lm(sqrt(Volume)~Girth,data=trees)
data.pre=data.frame(predict(model2,interval="prediction"))
lines(data.pre$lwr~trees$Girth,col="blue",lty=2)
lines(data.pre$upr~trees$Girth,col="blue",lty=2)
我们还可以将树围和树高都加入到模型中去,进行多元回归。如果要考虑的变量很多,可以用step函数进行变量筛选,它是以AIC作为评价指标来判断一个变量是否应该加入模型,建议使用这种自动判断函数时要谨慎。对于嵌套模型,还可以使用anova建立方差分析表来比较模型。对于变量变换的形式,则可以使用MASS扩展包中的boxcox函数来进行COX变换。
本文来自:http://www./lang/chinese/546
相关热词:
..........
生物秀是目前国内最具影响力的生物医药门户网站之一,致力于IT技术和BT的跨界融合以及生物医药领域前沿技术和成功商业模式的传播。为生物医药领域研究人员和企业提供最具价值的行业资讯、专业技术、学术交流平台、会议会展、电子商务和求职招聘等一站式服务。
官方微信号:shengwuxiu
电话:021-苹果/安卓/wp
苹果/安卓/wp
积分 675, 距离下一级还需 125 积分
权限: 自定义头衔, 签名中使用图片
道具: 彩虹炫, 雷达卡, 热点灯, 雷鸣之声, 涂鸦板, 金钱卡, 显身卡, 匿名卡, 抢沙发下一级可获得
权限: 隐身
购买后可立即获得
权限: 隐身
道具: 金钱卡, 雷鸣之声, 彩虹炫, 雷达卡, 涂鸦板, 热点灯
开心签到天数: 149 天连续签到: 1 天[LV.7]常住居民III
书接上回如果存在着严重的多重共线性,则需要使用合适的方法尽量地降低多重共线性,有两种比较常用的方法:一、逐步回归逐步回归主要分为向前逐步回归(forward)、向后逐步回归(backward)和向后向前逐步回归(both)。逐步回归本身并不是一种新的回归或者参数的估计方法,所用到的参数估计方法都是原来的,是从众多的变量中选出最优模型的变量的一套方法。即假如因变量Y ,4 个自变量分别是X1 ,X2 , X3 ,X4 。当所有自变量都进入到模型中时,Y =α +β1X1+β2X2 +β3X3 +β4X4 +μ 。现在是如何利用逐步回归方法从中选取最优的模型?向前逐步回归的思路是逐个引入变量。具体来讲是,先用因变量与每个自变量都进行回归,选取最优的模型,假如第一步选取的最优模型是Y =α +β1X1 +μ ;接着在第一步的最优模型的基础上,从剩余的变量X2,X3 ,X4 中每个分别加入到第一步的最优模型中,得Y =α +β1X1 +βj Xj +μ , j = 0,2,3,4, j = 0即为Y =α +β1X1 +μ ,比较这四个模型,如果发现模型Y =α +β1X1 +β3X3+μ 最优;接着再在第二步的最优模型Y =α +β1X1 +β3X3 +μ 上,从剩余的变量X2, X4中每个分别加入到第二步的最优模型中,得Y =α +β1X1+β3X3+βjXj+μ , j = 0,2,4,比较这三个模型,如果 j = 0时,模型最优,则最终选取的最优模型是Y =α +β1X1 +β3X3+μ向后逐步回归的思路是先引入全部自变量,然后逐个剔除不重要的变量,其剔除变量的思路和向前逐步回归的思路类似。向后向前逐步回归先逐步剔除变量,但可以后面的步骤中重新引入原先被剔除的变量,其方向是双向的,而向后逐步回归的自变量一旦被剔除后,在后面的步骤中就不会被重新引入,是单向的。注意,上文所指的最优模型一般通过一些准则来确定,比如F 值、可决系数R2、AIC 等。继续上篇提到的财政收入影响因素的例子:首先介绍一下step函数的用法,它是属于stats包,使用之前需先加载。step(object, scope,scale = 0,& &&&direction = c(&both&,&backward&, &forward&),& &&&trace = 1, keep = NULL, steps = 1000, k =2, ...)
14:06:40 上传
向前逐步回归的最优模型是把所有自变量都引入模型,没有剔除任何变量。
14:06:44 上传
向后逐步回归中,从AIC最小的变量依次逐步剔除了农业,建筑业,受灾三个变量,第四步不剔除变量时最优,即最终模型中包含工业,人口,消费三个变量。二、岭回归当解释变量之间存在多重共线性时,即X′X ≈ 0,则Var(βˆ) =σ 2 (X′X)−1将会增大,原因是X′X接近奇异。如果将X′X加上一个正常数对角阵λ I (λ & 0,I 为单位矩阵)即X′X +λ I,使得 X′X+λI ≈ 0的可能性比 X′X ≈ 0的可能性更小,那么X′X +λ I接近奇异的程度就会比X′X小的多,这就是岭回归的最初想法。R里MASS包的lm.ridge()函数可以用来做岭估计,其用法与lm()用法类似。可以证明β 的岭回归估计为βˆ (λ) = (X’X+λI)-1 X’Yλ 称为岭参数.岭估计的关键是选取岭参数λ,岭迹法是选取岭参数λ的常用方法之一。若记βˆ (λ)为βiˆ (λ )的第i个分量,它是λ 的一元函数。当λ 在[0,∞)上变化时,βˆ (λ)的图形称为岭迹(ridge trace)。βˆ (λ )的每个分量βj ˆ(λ )
的岭迹画在同一个图上,根据岭迹的变化趋势选择λ值,使得各个回归系数的岭估计大体上稳定,并且各个回归系数岭估计值的符号比较合理并符合实际。lm.r是属于MASS包的,用法和lm类似& lm.r&-lm.ridge(revenue~industry+agriculture+construction+consumption+pop+disaster,data=dat)& lm.r& && && && && && & industry& &agriculture&&construction& &consumption
6.&&1. -7.&&4.&&6. & && && & pop& && &disaster -7.&&4. 不指定λ值时,默认为0,结果和OLS一致。下面生成一个lambda序列,从0到0.3,间隔0.001,。同时把不同参数的估计值βˆ (λ )估计出来,画出岭迹图。如下:
14:06:46 上传
当λ取0.25-0.3之间时,参数的估计大致趋于稳定。& select(lm.ridge(revenue~industry+agriculture+construction+consumption+pop+disaster,data=dat,lambda=seq(0,0.3,0.001)))modified HKB estimator is 0. modified L-W estimator is 0. smallest value of GCV&&at 0.004 通过select函数可以选取更为精确的岭参数,本例中我们取λ=0.004& lm.ridge(revenue~industry+agriculture+construction+consumption+pop+disaster,data=dat,lambda=0.004)& && && && && && & industry& &agriculture&&construction& &consumption& && && &&&pop& && &disaster
5.&&1. -3.&&1.&&5. -5.&&4. 再代入到lm.ridge()函数中,就可以估计出相应的岭估计结果。本节完,下节开始讲异方差性问题。、
支持楼主:、
购买后,论坛将把您花费的资金全部奖励给楼主,以表示您对TA发好贴的支持
载入中......
本帖被以下文库推荐
& |主题: 90, 订阅: 11
[img]http://pic./forum//184005lbzbae6212af7azk.
如果是面板数据的 一组自变量 如何判断他们间的共线性
论坛好贴推荐
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
为做大做强论坛,本站接受风险投资商咨询,请联系(010-)
邮箱:service@pinggu.org
合作咨询电话:(010)
广告合作电话:(刘老师)
投诉电话:(010)
不良信息处理电话:(010)
京ICP证090565号
京公网安备号
论坛法律顾问:王进律师- Database Error
Discuz! Database Error
已经将此出错信息详细记录, 由此给您带来的访问不便我们深感歉意.

我要回帖

更多关于 r语言多元线性回归 的文章

 

随机推荐