线性回归是利用预测变量的一个线性组合函数来预测响应变量的同济分析方法,该线性回归模型的形式如下:
y=c0+c1x1+c2x2+......ckxk 其中, x1 , x2 , x3 ,…, xk 为预测变量, y 为对预测的响应变量。下面将在澳大利亚消费者价格知识CPI的数据上使用函数lm()做线性回归分析,该数据为2008到2010年澳大利亚的季度消费者价格指数。
首先,需要创建数据集并绘制散布图。在下面的代码中,使用函数axis()手动添加一个横坐标,参数las=3设置文字为垂直方向,如图5-1所示。
> year <- rep(2008:2010,each=4) > quarter<-rep(1:4,3) > cpi<-c(162.2,164.4,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 = "") > #draw x-axis > axis(1,labels=paste(year,quarter,sep="Q"),at=1:12,las=3)接下来,查看CPI与其他变量之间的相关系数,包括year(年份)和quarter(季度)这两个变量。
> quarter [1] 1 2 3 4 1 2 3 4 1 2 3 4 > cor(year,cpi) [1] 0.9106228 > cor(quarter,cpi) [1] 0.3739437接下来,在前面的数据上使用函数lm()建立一个线性回归模型,其中 year和quarter为预测变量,CPI为响应变量
> fit <- lm(cpi~year quarter) > fit Call: lm(formula = cpi ~ year quarter) Coefficients: (Intercept) year quarter -7694.746 3.913 1.173根据上面建立的线性模型,CPI的计算公式为:
cpi=c0 c1∗year c2∗quarter
其中,c0,c1,c2为拟合模型fit的系数。因此,2011年的CPI值可以计算如下,另一种更简单地计算CPI值的方法是使用函数preict(),这种方法将在本节的最后介绍
> (cpi2011 <- fit$coefficients[[1]]+fit$coefficients[[2]]*2011+fit$coefficients[[3]]*(1:4)) [1] 174.4650 175.6383 176.8117 177.9850该模型的更多细节可以通过下面的代码获得
> attributes(fit) $names [1] "coefficients" "residuals" "effects" "rank" "fitted.values" "assign" "qr" [8] "df.residual" "xlevels" "call" "terms" "model" $class [1] "lm" > fit$coefficients (Intercept) year quarter -7694.745833 3.912500 1.173333观测值与拟合结果的残差使用函数residuals()来计算 differences between observed values and fitted values
> residuals(fit) 1 2 3 4 5 6 7 8 9 10 -0.5275000 0.4991667 1.4258333 -0.2475000 -0.4400000 -0.8133333 -0.3866667 -0.6600000 0.4475000 0.3741667 11 12 0.4008333 -0.0725000 > summary(fit) Call: lm(formula = cpi ~ year + quarter) Residuals: Min 1Q Median 3Q Max -0.8133 -0.4619 -0.1600 0.4125 1.4258 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -7694.7458 506.0354 -15.206 1.00e-07 *** year 3.9125 0.2519 15.533 8.33e-08 *** quarter 1.1733 0.1840 6.379 0.000128 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7124 on 9 degrees of freedom Multiple R-squared: 0.9691, Adjusted R-squared: 0.9622 F-statistic: 141 on 2 and 9 DF, p-value: 1.61e-07接下来使用下面的代码绘制拟合模型的图像,如图5-2所示
plot(fit)我们还可以绘制拟合模型3D图像,下面代码中使用函数scatterplt3d()创建一个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"))