Rcmdr交互操作展示_11300730016

视频连接: http://v.youku.com/v_show/id_XODY4NTMzNDIw.html

数据下载:点击下载

library(Rcmdr)
#使用Rcmdr导入数据

#检查变量间的相关性
> viewHPI <- as.data.frame(HPI[,c("hpi","l","p","wb","gdp","hpy","fp")])
> cor(viewHPI)
hpi l p wb gdp hpy
hpi 1.00000000 0.511156476 0.120064983 0.45105676 -0.07376787 0.50160827
l 0.51115648 1.000000000 0.002601622 0.70621625 0.65617790 0.90806909
p 0.12006498 0.002601622 1.000000000 -0.02522495 -0.06459839 -0.02152737
wb 0.45105676 0.706216247 -0.025224949 1.00000000 0.69582627 0.93471110
gdp -0.07376787 0.656177895 -0.064598394 0.69582627 1.00000000 0.75129816
hpy 0.50160827 0.908069089 -0.021527366 0.93471110 0.75129816 1.00000000
fp -0.23800591 0.610231869 -0.102753947 0.65923370 0.90332033 0.69803395
fp
hpi -0.2380059
l 0.6102319
p -0.1027539
wb 0.6592337
gdp 0.9033203
hpy 0.6980339
fp 1.0000000
#可以大致看出,HPI和FP,GDP负相关;HPI和P,GDP的相关较低。

#散点矩阵图
#非对角线区域绘制变量间的散点图,并添加了平滑和线性拟合曲线。
#对角线显示了是否符合正态分布
Rcmdr> scatterplotMatrix(~fp+gdp+hpi+hpy+l+p+wb, reg.line=lm, smooth=TRUE,
Rcmdr+ spread=FALSE, span=0.5, id.n=0, diagonal = 'qqplot', data=HPI)
#数据并不很好地符合正态分布;直观可以发现HPI和P之间的关系不明显。

#建立线性模型2
Rcmdr> LinearModel.1 <- lm(hpi ~ fp + gdp + hpy + l + p + wb, data=HPI)
Rcmdr> summary(LinearModel.1)

Call:
lm(formula = hpi ~ fp + gdp + hpy + l + p + wb, data = HPI)

Residuals:
Min 1Q Median 3Q Max
-3.5641 -1.6290 -0.2723 1.1049 10.5510

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.382e+01 7.358e+00 -5.955 1.90e-08 ***
fp -4.882e+00 1.955e-01 -24.967 < 2e-16 ***
gdp 4.213e-05 3.046e-05 1.383 0.168841
hpy -1.023e+00 2.924e-01 -3.498 0.000624 ***
l 1.242e+00 1.684e-01 7.377 1.18e-11 ***
p 6.422e-10 1.159e-09 0.554 0.580337
wb 1.117e+01 1.660e+00 6.728 3.77e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.098 on 144 degrees of freedom
Multiple R-squared: 0.9491, Adjusted R-squared: 0.947
F-statistic: 447.8 on 6 and 144 DF, p-value: < 2.2e-16
#该模型可以解释95%的HPI方差,有很高的拟合度。变量GDP和P没有通过显著性检验,在下一个模型中考虑去除。

#学生化残差(残差分位数比较QQ图)
#可交互操作(点击出现label)
Rcmdr> qqPlot(LinearModel.1, simulate=TRUE, id.method="y", id.n=2)
83 114
150 151
#正残差:真实值-预测值>0,即模型低估了114号国家的HPI

> viewHPI["114",] #真实值
hpi l p wb gdp hpy fp
114 25.2 78.4 1759000 6.6 80944 57.7 11.68
> fitted(LinearModel.1)["114"] #预测值
114
14.64903
#比较真实值和预测值发现符合QQ图

#成分+残差图
Rcmdr> crPlots(LinearModel.1, span=0.5)
#结论:WB,L,FP,HPY四个变量符合线性,建模尚可,不必要添加曲线、多项式成分。GDP,P拟合糟糕。

#去除GDP,P变量,建立模型2
Rcmdr> LinearModel.2 <- lm(hpi ~ fp + hpy + l + wb, data=HPI)

Rcmdr> summary(LinearModel.2)

Call:
lm(formula = hpi ~ fp + hpy + l + wb, data = HPI)

Residuals:
Min 1Q Median 3Q Max
-3.7541 -1.5926 -0.3623 0.9260 11.2733

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -40.1325 6.6678 -6.019 1.35e-08 ***
fp -4.6696 0.1128 -41.393 < 2e-16 ***
hpy -0.8399 0.2521 -3.331 0.0011 **
l 1.1451 0.1479 7.745 1.47e-12 ***
wb 10.1956 1.4561 7.002 8.54e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.101 on 146 degrees of freedom
Multiple R-squared: 0.9483, Adjusted R-squared: 0.9469
F-statistic: 669.3 on 4 and 146 DF, p-value: < 2.2e-16
#这四个变量可以解释95%的HPI方差,有很高的拟合度,通过显著性检验
#着四个变量都通过了显著性检验

#观察第二个模型的残差分位数比较QQ图
Rcmdr> qqPlot(LinearModel.2, simulate=TRUE, id.method="identify", id.n=2)
[1] "75" "83" "114"
RcmdrMsg: [4] 注意: 警告: 已经找到了最近的点

#观察第二个模型的成分+残差图
Rcmdr> crPlots(LinearModel.2, span=0.5)
#结论:四个变量的线性拟合都较好

#比较两组模型
Rcmdr> anova(LinearModel.1, LinearModel.2)
Analysis of Variance Table

Model 1: hpi ~ fp + gdp + hpy + l + p + wb
Model 2: hpi ~ fp + hpy + l + wb
Res.Df RSS Df Sum of Sq F Pr(>F)
1 144 634.13
2 146 644.75 -2 -10.623 1.2062 0.3023
#结论是可以选用第二个更为精简的模型

#研究哪个变量对HPI有更大的影响,首先进行标准化
Rcmdr> HPI <- local({
Rcmdr+ .Z <- scale(HPI[,c("fp","hpi","hpy","l","wb")])
Rcmdr+ within(HPI, {
Rcmdr+ Z.wb <- .Z[,5]
Rcmdr+ Z.l <- .Z[,4]
Rcmdr+ Z.hpy <- .Z[,3]
Rcmdr+ Z.hpi <- .Z[,2]
Rcmdr+ Z.fp <- .Z[,1]
Rcmdr+ })
Rcmdr+ })
RcmdrMsg: [5] 注意: 此 HPI 資料集有 151 個列與 16 個欄.

#使用标准化后的变量建立模型3
Rcmdr> LinearModel.3 <- lm(Z.hpi ~ Z.fp + Z.hpy + Z.l + Z.wb, data=HPI)

Rcmdr> summary(LinearModel.3)

Call:
lm(formula = Z.hpi ~ Z.fp + Z.hpy + Z.l + Z.wb, data = HPI)

Residuals:
Min 1Q Median 3Q Max
-0.41178 -0.17469 -0.03974 0.10158 1.23657

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.431e-16 1.876e-02 0.000 1.0000
Z.fp -1.108e+00 2.676e-02 -41.393 < 2e-16 ***
Z.hpy -1.067e+00 3.203e-01 -3.331 0.0011 **
Z.l 1.232e+00 1.590e-01 7.745 1.47e-12 ***
Z.wb 1.309e+00 1.869e-01 7.002 8.54e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2305 on 146 degrees of freedom
Multiple R-squared: 0.9483, Adjusted R-squared: 0.9469
F-statistic: 669.3 on 4 and 146 DF, p-value: < 2.2e-16
#观察相关系数,可以发现WB相对有较高的影响。1标准差的WB变化对带来1.309标准差的HPI增加。

#差异值观测(数据处理或增减务必符合逻辑和理论)
#筛选离散值(Bonferonni离散值检验)-模型预测效果不佳的点
> View(HPI)
> library(car)
> outlierTest(LinearModel.2)
rstudent unadjusted p-value Bonferonni p
114 6.639766 5.8715e-10 8.8659e-08
83 3.825518 1.9335e-04 2.9195e-02

#可以直接在Rcmdr中该操作
Rcmdr> outlierTest(LinearModel.3)
rstudent unadjusted p-value Bonferonni p
114 6.639766 5.8715e-10 8.8659e-08
83 3.825518 1.9335e-04 2.9195e-02

#检测高杠杆值-与其他预测变量有关的离群点
Rcmdr> avPlots(LinearModel.3, id.method="mahal", id.n=2)

#检测强影响点-对模型有很大影响的点
Rcmdr> influencePlot(LinearModel.3, id.method="noteworthy", id.n=2)
StudRes Hat CookD
83 3.825518 0.1052653 0.5611946
114 6.639766 0.1545851 1.1157396
150 -1.188036 0.1099621 0.1864879
#influence plot纵轴可以观察+-2之外属于离散值;横轴>0.2/0.3属于高杠杆指;点大小表示影响力大小

plotEx.data.frame(…)函数编写与示范



由于在Rweb环境下编写函数无法完全展现,因此将所作图形贴于下方,分别为5个数据的不同类型的画图。 pied1 mosaicplot d2 boxplotd3 cdplot1 d4 cdplot2d5

基础R作图&ggplot2简单应用1

#简单绘图示例,使用ggplot2自带数据diamonds
library(ggplot2)
attach(diamonds)
diamond<-diamonds[1:200,]
names(diamonds)
plot(carat,price,type="b",ltyHPI Data=3,pch=4,cex=.5,xlim=c(0,3))
abline(lm(price~carat),col="red")
title(main="pirce and diamond size",col.main="red",xlab="carat",ylab="price",col.lab="green",cex.lab=0.5)
#数据来源Happy Planet Index (HPI)
#变量解释:HPIRank(HPI Rank);C(Country);R(Region);L(life expectancy);
#WB(Well-Being);HPY(happy years);FP(footprint);P(population);GRank(Government Rank)

#转入数据(已在SPSS中作初步调整)

#点击此下载整理后的HPI数据,请直接存放在F盘下方使用
library(foreign)
mydata=read.spss('F:/data.sav')
mydata

#从list变为dataframe
HPIDATA<-as.data.frame(mydata)
HPIDATA

#GDP为x轴,HPI为y轴作出散点和平滑的拟合曲线
library(ggplot2)
attach(HPIDATA)
p1<-ggplot(data=HPIDATA,aes(x=GDP,y=HPI,))
p1+geom_point()+geom_smooth()
detach(HPIDATA)

#将变量名改短一点方便后期分类
names(HPIDATA)[1]<-"HPIR"
names(HPIDATA)[11]<-"GR"

#方便按照类别查看可以按照GR排序
#attach(HPIDATA)
#HPI_GR<-HPIDATA[order(GRank),]
#detach(HPIDATA)
#HPI_GR

#将Government Rank分类为High Government Rank(HGR),Middle~(MGR),Low~(LGR)
#将变量命名为GRCAT(Government Rank Catogory)并区分颜色显示
attach(HPIDATA)
HPIDATA$GRCAT[HPIDATA$GR<51]<-"HGR"
HPIDATA$GRCAT[HPIDATA$GR>=51&HPIDATA$GR<101]<-"MGR"
HPIDATA$GRCAT[HPIDATA$GR>=101]<-"LGR"

#以上分类还可以通过以下代码完成:
#HPIDATA<-within(HPIDATA),{
# GRCAT<-NA
# GRCAT[GR<51] <-"HGR"
# GRCAT[GR>=51&GR<101] <-"MGR"
# GRCAT[GR>=101] <-"LGR"
}

#用color=factor(GRCAT)区分颜色,alpha=0.5表示透明度
library(ggplot2)
as.factor(GRCAT)
p2<-ggplot(data=HPIDATA,aes(x=GDP,y=HPI,color=factor(GRCAT),alpha=0.5))
p2+geom_point()+geom_smooth()

#用size=WB区分点的大小
p3<-p2 + geom_point(aes(colour=factor(GRCAT),size=WB))+geom_smooth()
p3