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)
