R語言實戰讀書筆記(八)回歸

簡單線性：用一個量化驗的解釋變量預測一個量化的響應變量

多項式：用一個量化的解決變量預測一個量化的響應變量，模型的關系是n階多項式

多元線性：用兩個或多個量化的解釋變量預測一個量化的響應變量

多變量：用一個或多個解釋變量預測多個響應變量

Logistic：用一個或多個解釋變量預測一個類別型響應變量

泊松：用一個或多個解釋變量預測一個代表頻數的響應變量

Cox比例風險：用一個或多個解釋變量預測一個事件發生的時間

時間序列：對誤差項相關的時間序列數據建模

非線性：用一個或多個量化的解釋變量預測一個量化的響應變量，不過模型是非線性的

非參數：用一個或多個量化的解釋變量預測一個量化的響應變量，模型的形式源自數據形式，不事先設定

穩健：用一個或多個量化的解釋變量預測一個量化的響應變量，能抵御強影響點的干擾

8.1.1

OLS回歸：普通最小二乘回歸法：包括簡單線性，多項式回歸和多元線性回歸。
OLS回歸是通過預測變量的加權和來預測量化的因變量，其中權重是通過數據估計而得的參數。

例：一個工程師想找到跟橋梁退化有關的最重要的因素，比如使用年限，交通流量，橋梁設計，建造材料和建造方法，建造質量以及天氣情況，並確定它們之間的數學關系。他擬合了一系列模型，檢驗它們是否符合相應的統計假設，探索了所有異常的發現，最終從許多可能的模型中選擇了最佳的模型，如果成功，那么結果將會幫助他完成以下任務。

1.在眾多變量中判斷哪些對預測橋梁退化是有用的，得到它們的相對重要性，從而關注重要的變量。

2.根據回歸所得的等式預測新的橋梁的退化情況，找出那些可能會有麻煩的橋梁

3.利用對異常橋梁的分析，獲得一些意外的信息。比如他發現某些橋梁的退化速度比預測的更快或更慢，那么研究這些“離群點”可能會有重大的發現，能夠幫助理解橋梁退化的機制。

OLS回歸擬合模型的形式：

Y=B0+B1X1+B2X2+...+BnXn

目標是通過減少響應變量的真實值與預測值的差值來獲得模型參數，使得殘差平方和最小。

為了能夠恰當解釋OLS模型的系數，數據必須滿足以下統計假設：

1.正態性：對於固定的自變量值，因變量值成正態分布。

2.獨立性：Y值之間相互獨立

3.線性：因變量與自變量成線性關系

4.同方差性：因變量的方差不隨自變量的水平不同而變化，也可稱為不變方差，但是說同方差性感覺上更犀利。

如果違背以上假設，你的統計顯著性檢驗結果和所得的置信區間很可能就不精確。

8.2.1 用lm擬合回歸模型

summary:

coefficients：列出擬合模型的模型參數（截距和斜率）

confint：提供模型參數的置信區間

fitted：列出擬合模型的預測值

residuals：列出擬合模型的殘差值

anova：生成一個擬合模型的方差分析表，或者比較兩個或更多擬合模型的方差分析表

vcov：列出模型參數的協方差矩陣

AIC：輸出AIC統計量

plot：

predict：

簡單線性回歸：一個因變量和一個自變量

多項式回歸：一個預測變量，同時包含變量的冪

多元線性回歸：不止一個預測變量

8.2.2 線性回歸

fit <- lm(weight ~ height, data = women)
summary(fit)
fitted(fit)  #列出預測值
residuals(fit) #列出殘差

plot(women$height, women$weight, main = "Women Age 30-39", xlab = "Height (in inches)", ylab = "Weight (in pounds)")
abline(fit)

8.2.3 多項式回歸

線性回歸的圖表明，可以通過添加一個二次項X平方來提高回歸的預測精度

fit2 <- lm(weight ~ height + I(height^2), data = women)
summary(fit2)

plot(women$height, women$weight, main = "Women Age 30-39", xlab = "Height (in inches)", ylab = "Weight (in lbs)")

lines(women$height, fitted(fit2))

library(car)
scatterplot(weight ~ height, data = women, spread = FALSE, lty.smooth = 2, pch = 19, main = "Women Age 30-39", xlab = "Height (inches)", ylab = "Weight (lbs.)")

8.2.4 多元線性回歸

states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost")])
cor(states)

library(car)
scatterplotMatrix(states, spread = FALSE, lty.smooth = 2, main = "Scatterplot Matrix")

8.2.5 有交互項的多遠線性回歸

#ht和wt

fit <- lm(mpg ~ hp + wt + hp:wt, data = mtcars)
summary(fit)

8.3 回歸診斷

8.3.1 標准方法

fit <- lm(weight ~ height, data = women)
par(mfrow = c(2, 2))
plot(fit)

OLS的回歸統計假設，

1.正態性：當預測變量值固定時，因變量成正態分布，則殘差值也應該是一個均值為0的正態分布，上圖的右上角。

2.獨立性：滿足

3.線性：若因變量與自變量線性相關，那么殘差值與預測（擬合）值就沒有任何系統關聯。在“殘差圖與擬合圖”（左上）可以清楚地看到一個曲線關系，這說明你可能需要對回歸模型加上一個二次項。

4.同方差性：若滿足不變方差假設，那么在位置尺度圖（左下）中，水平線周圍的點應該隨機分布。該圖滿足此假設。

最后一幅圖提供了可能關注的單個觀測點的信息。從圖形可以鑒別出離群點，高杠桿值點和強影響點。

1.一個觀測點是離群點，表明擬合回歸模型對其預測效果不佳（產生巨大的或正或負的殘差）

2.一個觀測點有很高的杠桿值，表明它是一個異常的預測變量值的組合。也就是說，在預測變量空間中，它是一個離群點。因變量值不參與計算一個觀測點的杠桿值。

3.一個觀測點是強影響點，表明它對模型參數的估計產生的影響過大，非常不成比例。強影響點可以通過Cook距離即Cook's D統計量來鑒別。

以下對模型加了一個平方項

newfit <- lm(weight ~ height + I(height^2), data = women)
par(mfrow = c(2, 2))
plot(newfit)

與上圖的差別，左上圖應該沒有任何關聯，右上圖應該盡可能落在45度角的直線上。左下圖隨機分布，右下識別離群點強影響點。從左上和右上圖看效果比較好。

刪除離群點

newfit <- lm(weight ~ height + I(height^2), data = women[-c(13, 15),])
par(mfrow = c(2, 2))
plot(newfit)

換個例子

fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
par(mfrow = c(2, 2))
plot(fit)

8.3.2 改進的方法

library(car)
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
qqPlot(fit, labels = FALSE, simulate = TRUE, main = "Q-Q Plot")

如下圖，右上角那個點在置信區間之外，表明模型低估了該州的謀殺率

residplot <- function(fit, nbreaks=10) {
z <- rstudent(fit)
hist(z, breaks=nbreaks, freq=FALSE,xlab="Studentized Residual",main="Distribution of Errors")
rug(jitter(z), col="brown")
curve(dnorm(x, mean=mean(z), sd=sd(z)),
add=TRUE, col="blue", lwd=2)
lines(density(z)$x, density(z)$y,
col="red", lwd=2, lty=2)
legend("topright",
legend = c( "Normal Curve", "Kernel Density Curve"),lty=1:2, col=c("blue","red"), cex=.7)
}

residplot(fit)

2.誤差的獨立性

durbinWatsonTest(fit)

#lag Autocorrelation D-W Statistic p-value
#1 -0.2006929 2.317691 0.248
#Alternative hypothesis: rho != 0

p值0.282表明無自相關性，誤差項之間獨立，lag為1表明數據集中每個數據都是與其后一個數據進行比較的。該檢驗適用於時間獨立的數據，對於非聚集型數據並不適用。

3.線性

crPlots(fit, one.page = TRUE, ask = FALSE)

若圖形存在非線性，說明需要添加一些曲線成分，比如多項式項，或對一個或多個變量進行變換。

4.同方差性

library(car)
ncvTest(fit)
spreadLevelPlot(fit)

應該呈隨機分布，滿足方差不變假設，若不滿足，會看到一個非水平的曲線。

8.3.3 線性模型假設的綜合驗證

library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

#Value p-value Decision
#Global Stat 2.7728 0.5965 Assumptions acceptable.
#Skewness 1.5374 0.2150 Assumptions acceptable.
#Kurtosis 0.6376 0.4246 Assumptions acceptable.
#Link Function 0.1154 0.7341 Assumptions acceptable.
#Heteroscedasticity 0.4824 0.4873 Assumptions acceptable.

8.3.4 多重共線性

vif(fit)

sqrt(vif(fit)) > 2 #>2說明有多重共線性

Population Illiteracy Income Frost
FALSE FALSE FALSE FALSE

8.4 異常觀測值

8.4.1 離群點

library(car)
outlierTest(fit)

8.4.2 高杠桿值點

hat.plot <- function(fit){
p <- length(coefficients(fit))
n <- length(fitted(fit))
plot(hatvalues(fit), main = "Index Plot of Hat Values")
abline(h = c(2, 3) * p/n, col = "red", lty = 2)
identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}

hat.plot(fit)

8.4.3 強影響點

cutoff <- 4/(nrow(states) - length(fit$coefficients) - 2)
plot(fit, which = 4, cook.levels = cutoff)
abline(h = cutoff, lty = 2, col = "red")

avPlots(fit, ask = FALSE, onepage = TRUE, id.method = "identify")

influencePlot(fit, id.method = "identify", main = "Influence Plot", sub = "Circle size is proportial to Cook's Distance")

8.5 改進措施

8.5.1 刪除觀測點

8.5.2 變量變換

library(car)
summary(powerTransform(states$Murder))

結果如下，0.6055說明應該用Murder的0.6次方來代替，可以使用0.5次方即平方根，但lambda=1的P值無法拒絕原假設，因此需要變量變換

bcPower Transformation to Normality

　　　　　　　　Est.Power Std.Err. Wald Lower Bound Wald Upper Bound
states$Murder 0.6055 0.2639 0.0884 1.1227

Likelihood ratio tests about transformation parameters
　　　　　　　　　　　　LRT 　　df pval
LR test, lambda = (0) 5.665991 1 0.01729694
LR test, lambda = (1) 2.122763 1 0.14512456

 

library(car)
boxTidwell(Murder ~ Population + Illiteracy, data = states)

結果如下：表明應該是Population的0.87次方，Illiteracy的1.36次方，但是P值>0.05又表明不需要變換

　　　　　　　　Score Statistic 　p-value 　　　　MLE of lambda
Population 　-0.3228003 　　　　0.7468465 　　 0.8693882
Illiteracy 　　 0.6193814 　　　　0.5356651 　　 1.3581188

iterations = 19

8.5.3 增刪變量

如果只是做預測，那么多重共線性不是問題，但是如果還要做解釋，那么就必須解決問題。

8.6 選擇“最佳”的回歸模型

8.6.1 模型比較

fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
fit2 <- lm(Murder ~ Population + Illiteracy, data = states)

anova(fit2, fit1)

結果如下：p值>0.05說明不需要將這兩個變量添加到線性模型中，可以刪除

Model 1: Murder ~ Population + Illiteracy
Model 2: Murder ~ Population + Illiteracy + Income + Frost
Res.Df 　RSS 　　Df 　　Sum of Sq 　　 　　 F Pr(>F)
1 47 　　289.25
2 45 　　289.17  2 　　 0.078505 0.0061 　　0.9939

fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
fit2 <- lm(Murder ~ Population + Illiteracy, data = states)
AIC(fit1, fit2)

結果如下：AIC值越小的模型要優先選擇，所以選模型2

　　df 　　AIC
fit1 6 　　241.6429
fit2 4 　　237.6565

8.6.2 變量選擇

1.逐步回歸

library(MASS)
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
stepAIC(fit, direction = "backward")

2.全子集回歸

library(leaps)
leaps <- regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data = states, nbest = 4)
plot(leaps, scale = "adjr2")

library(car)
subsets(leaps, statistic = "cp", main = "Cp Plot for All Subsets Regression")
abline(1, 1, lty = 2, col = "red")

8.7 深層次分析

8.7.1 交叉驗證：在k重交叉驗證中，樣本被分成k個子樣本，輪流將k-1個子樣本組合作為訓練集，另外1個子樣本作為測試集。記錄k個結果，求平均值。

shrinkage <- function(fit, k = 10) {
require(bootstrap)
# define functions
theta.fit <- function(x, y) {
lsfit(x, y)
}
theta.predict <- function(fit, x) {
cbind(1, x) %*% fit$coef
}

# matrix of predictors
x <- fit$model[, 2:ncol(fit$model)]
# vector of predicted values
y <- fit$model[, 1]

results <- crossval(x, y, theta.fit, theta.predict, ngroup = k)
r2 <- cor(y, fit$fitted.values)^2
r2cv <- cor(y, results$cv.fit)^2
cat("Original R-square =", r2, "\n")
cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
cat("Change =", r2 - r2cv, "\n")
}

# using shrinkage()

fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data = states)
shrinkage(fit)

#Loading required package: bootstrap
#Original R-square = 0.5669502 初始結果
#10 Fold Cross-Validated R-square = 0.4562456 交叉驗證結果
#Change = 0.1107046

fit2 <- lm(Murder ~ Population + Illiteracy, data = states)
shrinkage(fit2)

#Original R-square = 0.5668327
#10 Fold Cross-Validated R-square = 0.5171145
#Change = 0.0497182  比全變量減少的R平方，說明這個比較好

8.7.2相對重要性

zstates <- as.data.frame(scale(states))
zfit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data = zstates)
coef(zfit)

結果如下：說明Illiteracy是最重要的，Frost是最不重要的（最小）

(Intercept) 　　Population 　　Income 　　　　Illiteracy 　　　　Frost
-9.405788e-17 2.705095e-01 1.072372e-02 6.840496e-01 8.185407e-03

 

relweights <- function(fit, ...) {
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))

# correlations between original predictors and new orthogonal variables
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda^2

# regression coefficients of Y on orthogonal variables
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta^2)
rawwgt <- lambdasq %*% beta^2
import <- (rawwgt/rsquare) * 100
lbls <- names(fit$model[2:nvar])
rownames(import) <- lbls
colnames(import) <- "Weights"

# plot results
barplot(t(import), names.arg = lbls, ylab = "% of R-Square",
xlab = "Predictor Variables", main = "Relative Importance of Predictor Variables",
sub = paste("R-Square = ", round(rsquare, digits = 3)),
...)
return(import)
}

# using relweights()

fit <- lm(Murder ~ Population + Illiteracy + Income +
Frost, data = states)
relweights(fit, col = "lightgrey")