R语言线性回归分析实战：模型构建与假设检验

R语言线性回归分析实战：模型构建与假设检验

本文将介绍如何使用R语言进行线性回归分析，包括模型构建、参数估计、模型预测以及模型假设的验证。通过一个具体的案例，帮助读者理解线性回归分析的方法和技巧。

数据来源及背景

数据来源于一个玻璃制造公司的案例，该公司主要向新建筑承包商和汽车公司供应产品。公司认为其年销售额应与新建筑数量以及汽车生产高度相关，因此希望通过构建线性回归模型来预测销售额。

glass <- read.csv("glass_mult.csv",header=T)
glass
summary(glass)

应用

模型构建

首先查看变量之间的相关性，发现自变量BLDG（建筑数量）与因变量SALES（销售额）存在高达0.948的正相关性。

# explore the correlation between variables
cor(glass)

绘制任意两个变量之间的散点图来探索其相关关系。

# get all scatter plots between 2 variables
library(ggplot2)
library(GGally)
ggpairs(glass)
## alternative
## pairs( ~ SALES + BLDG + AUTO, data = glass)

由于BLDG与SALES高度相关，先构建二者之间的一元线性回归模型。

# build linear regression model 1
glass.lm1 <- lm(SALES ~ BLDG, data=glass)
summary(glass.lm1)
anova(glass.lm1)

根据模型结果，BLDG高度显著，且R-squared为0.8993，Adjusted R-squared为0.8926，说明该模型解释了将近90%的variance。

绘制BLDG与SALES的散点图，回归曲线，以及置信区间。

# plot the model as well as the real data  
ggplot(glass, aes(x = BLDG, y = SALES)) + geom_point() + geom_smooth(method = 'lm')   
# could add se = 'F' to geom_smooth to remove standard error lines   
# alternative1  
# ggplot(glass, aes(x = BLDG, y = SALES)) + geom_point() +   
# geom_abline(intercept = coef(glass.lm1)[1], slope = coef(glass.lm1)[2])

# alternative2
## plot(glass$BLDG, glass$SALES)
## abline(glass.lm1)

为了进一步理解一元线性回归模型，接下来手动构建一元线性回归模型并计算其参数（slope and intercept）。

# manually fit the simple linear regression model
library(dplyr)
glass %>%
 summarise(
 N = n(),
 r = cor(SALES, BLDG),
 mean_x = mean(BLDG),
 sd_x = sd(BLDG),
 mean_y = mean(SALES),
 sd_y = sd(SALES),
 slope = r * sd_y / sd_x,
 intercept = mean_y - slope * mean_x
 )

接下来，构建BLDG与AUTO（汽车生产数量）的SALES模型。

# build linear regression model 2
glass.lm2 <- lm(SALES ~ BLDG+AUTO, data=glass)
summary(glass.lm2)
anova(glass.lm2)

根据模型结果，两个自变量均高度显著，且R-squred及Adjusted R-squred分别提高至0.9446和0.939。

探索模型二是否在模型一的基础上具有显著的提升。由于P-value小于0.05，我们可以认为模型二与模型一相比有显著提升。

# check model improvement
anova(glass.lm1, glass.lm2) # significant p-value means significant improvement

基于已有数据，运用模型二对其进行预测，并绘制真实值与预测值的散点图来观察预测准确性。

# model prediction
glass$lm2.pred <- predict(glass.lm2)
# compare the model prediction and the real data points
ggplot(glass, aes(x = lm2.pred, y = SALES)) + geom_point() + geom_abline()

检验模型是否符合线性回归的假设

1. 假设1：自变量之间是独立的 (independence)

该假设可以通过3.1中探索变量之间的相关性来验证。若自变量之间的相关性小于0.7，则可认为符合假设。

2. 假设2：自变量与因变量之间为线性/可加性的关系 (linearity)

该假设可以通过绘制散点图来判断

3. 假设3：残差符合正态分布 (normality)

# test for normality  
hist(glass.lm2$residuals, main = 'Histogram of Residual')
qqnorm(glass.lm2$residuals)
qqline(glass.lm2$residuals)  
# alternative  
## dat1 <- as.data.frame(glass.lm2$residuals)  
## names(dat1) <- 'res'  
## theme_set(  
## theme_minimal() +  
## theme(legend.position = "top"))  
## ggplot(dat1,aes(sample = res)) + stat_qq()

# Shapiro–Wilk test
# H0: 样本数据与正态分布没有显著区别。
# HA: 样本数据与正态分布存在显著区别。
shapiro.test(glass.lm2$residuals)

# roughly achieve qqplot by hand 
par(mfrow=c(1,1))
t <- rank(glass.lm2$residuals)/length(glass.lm2$residuals) #求观察累积概率
q <- qnorm(t) #用累积概率求分位数值
plot(q,glass.lm2$residuals)

4. 假设4：残差满足同方差性 (homoscedasticity)

# test for homoscedasticity
par(mfrow=c(1,2))
plot(glass.lm2$fitted.values, glass.lm2$residuals)
zres <- rstandard(glass.lm2)
plot(glass.lm2$fitted.values, zres)

5. 假设5： 残差满足独立性 (independence)

# test for independence
par(mfrow=c(1,1))
data <- data.frame(YEAR=c(1:17))
newglassdata <- cbind(glass,data)
#newglassdata
plot(newglassdata$YEAR, glass.lm2$residuals)

# check model assumption in one step
library(ggfortify)
autoplot(glass.lm2)
# alternative
# par(mfrow = c(2,2))
# plot(glass.lm2)