数学日记

数学是科学的皇后 - 高斯

 

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

Sirius:

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

> 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
评论
热度(8)

© 数学日记 | Powered by LOFTER