1. 目的：构建线性回归模型并检验其假设是否成立。
2. 数据来源及背景
2.1 数据来源：数据为本人上课的案例数据，
2.2 数据背景：“玻璃制造公司”主要向新建筑承包商和汽车公司供应产品。该公司认为，他们的年销售额应与新建筑数量以及汽车生产高度相关，因此希望构建线性回归模型来预测其销售额。
glass <- read.csv("glass_mult.csv",header=T)
glass
summary(glass)
3. 应用
3.1 模型构建
查看变量之间相关性，发现自变量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()
3.2 检验模型是否符合线性回归的假设
假设1: 自变量之间是独立的 (independence)
该假设可以通过3.1中探索变量之间的相关性来验证。若自变量之间的相关性小于0.7，则可认为符合假设。
假设2：自变量与因变量之间为线性/可加性的关系 (linearity)
该假设可以通过绘制散点图来判断
假设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：残差满足同方差性 (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： 残差满足独立性 (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)
