先回顾以下线性回归的概念:线性回归是利用预测变量的一个线性组合函数来预测响应变量的统计分析方法,该线性回归模型的形式如下:
其中,x1,x2,...,xk为预测变量,y为对预测的响应变量。
下面将在澳大利亚消费者价格指数(CPI)的数据上使用lm()函数做线性回归分析,该数据为2008年到2010年澳大利亚的季度消费者价格指数。
首先,需要创建数据集并绘制散点图,下面的代码中,使用axis()函数手动添加一个横坐标,参数las=3设置文字为垂直方向。
> quarter <- rep(1:4, 3)
> 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="")
> axis(1, labels=paste(year,quarter, seq="Q"), at=1:12, las=3)
接下来,查看CPI与其他变量之间的相关系数,包括自考year(年份)和quarter(季度)这两个变量。
> cor(year,cpi)
[1] 0.9096316
> cor(quarter,cpi)
[1] 0.3738028
接下来,在前面的数据上使用lm()函数建立一个线性回归模型,其中year和quarter为预测变量,CPI为响应变量。
> fit <- lm(cpi~year+quarter)
> fit
Call:
lm(formula=cpi ~ year + quarter)
Coefficients:
(Intercept) year quarter
-7644.488 3.888 1.167
根据上面建立的线性模型,CPI的计算公式为:
cpi=c0 + c1*year+c2*quarter
其中,c0,c1,c2为拟合模型fit的系数。因此,2011年的CPI值可以计算如下。另一种更简单的计算CPI值的方法是使用predict()函数 。
> (cpi2011 <- fit$coefficients[[1]] + fit$coefficients[[2]]*2011 + fit$coefficients[[3]]*(1:4))
[1] 174.4417 175.6083 176.7750 177.9417
这个模型的更多细节可以通过下面的代码获得。
> attributes(fit)
$names
[1] "coefficients" "residuals" "effects" "rank" "fitted.values"
[6] "assign" "qr" "df.residual" "xlevels" "call"
[11] "terms" "model"
$class
[1] "lm"
> fit$coefficients
(Intercept) year quarter
-7644.487500 3.887500 1.166667
观测值与拟合结果的残差使用residuals()函数来计算。
> residuals(fit)
1 2 3 4 5 6 7 8 9 10 11 12
-0.57916667 0.65416667 1.38750000 -0.27916667 -0.46666667 -0.83333333 -0.40000000 -0.66666667 0.44583333 0.37916667 0.41250000 -0.05416667
> summary(fit)
Call:
lm(formula=cpi ~ year + quarter)
Residuals:
Min 1Q Median 3Q Max
-0.8333 -0.4948 -0.1667 0.4208 1.3875
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7644.4875 518.6543 -14.739 1.31e-07 ***
year 3.8875 0.2582 15.058 1.09e-07 ***
quarter 1.1667 0.1885 6.188 0.000161 ***
---
Signif. codes: 0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1
Residual standard error: 0.7302 on 9 degrees of freedom
Multiple R-squared: 0.9672, Adjusted R-squared: 0.9599
F-statistic: 132.5 on 2 and 9 DF, p-value: 2.108e-07
接下来,使用下面的代码绘制拟合模型的图像。
>plot(fit)
我们还可以绘制拟合模型的3D图像,下面的代码中使用函数scatterplot3d()创建一个3D散点图,并使用plane3d()函数绘制拟合平面,参数lab指定x轴和y轴上的刻度。
> library(scatterplot3d)
> s3d <- scatterplot3d(year,quarter,cpi,highlight.3d=T,type="h", lab=c(2,3))
> s3d$plane3d(fit)
基于拟合模型,2011年的CPI可以通过如下方式预测,预测值用小三角表示。
> data2011 <- data.frame(year=2011,quarter=1:4)
> cpi2011 <- predict(fit, newdata=data2011)
> style <- c(rep(1,12),rep(2,4))
> plot(c(cpi,cpi2011),xaxt="n",ylab="CPI",xlab="",pch=style,col=style)
> axis(1,at=1:16,las=3,labels=c(paste(year,quarter,sep="Q"),"2011Q1","2011Q2","2011Q3","2011Q4"))
联系客服