Sirius

这里是风(Ventus) 捣鼓各种软件的文科兽控大学汪。

[R]回归分析pt.1——线性回归模型案例

跟着书上的东西做了一个案例。

> year=rep(2008:2010,each=4)
> year
 [1] 2008 2008 2008 2008 2009 2009 2009 2009 2010 2010 2010 2010

> quarter=rep(1:4,3)
> quarter
 [1] 1 2 3 4 1 2 3 4 1 2 3 4

> cpi=c(162.2,164.6,166.5,166.0,
+ 166.2,167.0,168.6,169.5,
+ 171.0,172.1,173.3,174.0)

上面三句很明白~

> plot(cpi,xaxt="n",ylab="CPI",xlab="")

//在图上描出cpi的12个值;xaxt="n"表示x轴不标明刻度,如果不写就是默认的1至12;y轴标题"CPI";x轴标题空
> #draw x-axis

//目前意义不明就背下来了。。
> axis(1,labels=paste(year,quarter,sep="Q"),at=1:12,las=3)

//1下2左3上4右;labels就是指的刻度上写的东西吧。。paste函数把year和quarter的数据连接在一起,中间用“Q”隔开;at指明在哪里放置labels;las=3设置文字为竖直方向(经过测试0和1是横向,2和3是竖向,不知道其区别。。)

这样我们就有了下图~



接下来似乎就是要拟合了> >..

> fit=lm(cpi~year+quarter)
> fit

Call:
lm(formula = cpi ~ year + quarter)

Coefficients:
(Intercept)         year      quarter  
  -7644.487        3.887        1.167  

//先来单词解释(x

coefficients  系数

intercept  截距

线性回归模型的公式:y=c0+c1x1+c2x2+...ckxk

> cpi2011=fit$coefficients[1]+fit$coefficients[2]*2011+fit$coefficients[3]*(1:4)
> cpi2011
[1] 174.4417 175.6083 176.7750 177.9417
//于是代入公式,我们就能求出2011年的cpi的预测值啦

> attributes(fit)
//通过这样可以看到更多细节

$names
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        

$class
[1] "lm"
//只需要fit$xxxx就好啦

如果需要也可以summary(fit)一下owo~

residual  残余(观测值-拟合值)

> plot(fit) //这样就得到了下面四个图w


//同样并没有看懂全部。。

> library(scatterplot3d)

//加载or导入程序包
> s3d=scatterplot3d(year,quarter,cpi,highlight.3d=T,type="h",lab=c(2,3))

//year,quarter,cpi分别用x,y,z标出;highlight.3d=T颜色表示;type="h"画出垂直于x-y平面的直线;lab大概就是指间隔的个数吧。。
> s3d$plane3d(fit)
//似乎是拟合平面
结果如下图~

> data2011=data.frame(year=2011,quarter=1:4)
> data2011
  year quarter
1 2011       1
2 2011       2
3 2011       3
4 2011       4

> cpi2011=predict(fit,newdata=data2011)
> cpi2011
       1        2        3        4
174.4417 175.6083 176.7750 177.9417
//就是代入公式

> style=c(rep(1,12),rep(2,4))
> style
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2

> plot(c(cpi,cpi2011),xaxt="n",ylab="CPI",xlab="",pch=style,col=style)
//1代表圆形,2代表三角形,更多详细见下图pch


col 1代表黑色,2代表红色

> axis(1,at=1:16,las=3,labels=c(paste(year,quarter,sep="Q"),"2011Q1","2011Q2","2011Q3","2011Q4"))

完工~最后成图

碎觉碎觉,快一点了= =

评论
热度(8)
  1. 點滴的記錄Sirius 转载了此文字  到 数学日记

© Sirius | Powered by LOFTER