預測懷俄明州蛇河流域的水量狡刘，數(shù)據(jù)集snake可以加載alr3包得到扶欣。
install.packages("alr3")
library(alr3)
data("snake")
str(snake)
'data.frame': 17 obs. of 2 variables:
Y: num 10.5 16.7 18.2 17 16.3 10.5 23.1 12.4 24.9 22.8 ...
head(snake)
X Y
1 23.1 10.5
2 32.8 16.7
3 31.8 18.2
4 32.0 17.0
5 30.4 16.3
6 24.0 10.5
更改變量名
names(snake) <- c("content","yield")
str(snake)
'data.frame': 17 obs. of 2 variables:
yield : num 10.5 16.7 18.2 17 16.3 10.5 23.1 12.4 24.9 22.8 ...
with(snake,plot(content,yield,xlab = "water content of snow",ylab = "water yield",las = 1))
snake01.png
散點圖顯示content和yield之間存在線性關系稿静，但首尾疑似存在兩個離群點。
建立線性回歸
yield.fit <- lm(yield~content,data = snake)
summary(yield.fit)
Call:
lm(formula = yield ~ content)
Residuals:
Min 1Q Median 3Q Max
-2.179 -1.515 -0.362 1.628 3.197
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.7254 1.5488 0.47 0.65
content 0.4981 0.0495 10.06 4.6e-08 ***
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.74 on 15 degrees of freedom
Multiple R-squared: 0.871, Adjusted R-squared: 0.862
F-statistic: 101 on 1 and 15 DF, p-value: 4.63e-08
P值高度顯著芯勘，可以拒絕原假設滑频。
回到剛剛的散點圖，為散點圖添加線性回歸模型產(chǎn)生的擬合直線。
with(snake,plot(content,yield))
abline(yield.fit,lwd=3,col="red")
線性回歸必須通過統(tǒng)計假設檢驗帅容。
正態(tài)性：對于固定的自變量值颇象，因變量值呈正態(tài)分布。
獨立性：Y值之間相互獨立并徘。
線性：因變量和自變量之間為線性相關遣钳。
同方差性：因變量的方差不隨自變量的水平不同而變化。
對模型進行回歸診斷
par(mfrow=c(2,2))
plot(yield.fit)
標準方法
正態(tài)性：右上角QQ圖是在正態(tài)分布對應的值下麦乞，標準化殘差的概率圖蕴茴，若滿足正態(tài)假設，那么圖上的點應該落在呈45度角的直線上姐直。
獨立性：從收集的數(shù)據(jù)來驗證倦淀。
線性：左上角殘差與擬合圖中，殘差值和擬合值不存在任何系統(tǒng)的關聯(lián)声畏。
同方差性：左下角位置尺度圖中撞叽，水平線的點應該隨機分布。
改進的方法
正態(tài)性：
(1)car包qqPlot()函數(shù)
library(car)
qqPlot(yield.fit,labels=row.names(snake),id.method="identify",simulate=TRUE,main="Q-Q Plot")
(2)學生化殘差圖
學生化殘差圖
residplot <- function(fit,nbreaks=10){
- z <- rstudent(fit)
- hist(z,breaks = nbreaks,freq = FALSE,
xlab = "Studentized Residual",xlim = c(-3,3),main = "Distribution of Errors")- rug(jitter(z),col = "brown")
- curve(dnorm(x,mean = mean(z)),add = TRUE,col = "blue",lwd = 2)
- lines(density(z)
y,col="red",lwd=2,lty=2)
- legend("topright",
legend = c("Normal Curve","Kernel Density Cruve"),lty = 1:2,col=c("blue","red"),cex=.7)- }
residplot(yield.fit)
誤差的獨立性
car包的Durbin-Watson檢驗砰识。
durbinWatsonTest(yield.fit)
lag Autocorrelation D-W Statistic p-value
1 -0.4152 2.706 0.106
Alternative hypothesis: rho != 0
P值等于0.106不顯著能扒，說明無自相關性佣渴。
線性
car包的crPlots()辫狼，繪制成分殘差圖。
crPlots(yield.fit)
snake06.png
若圖形存在非線性辛润，則說明預測變量的函數(shù)形式建模不夠充分膨处，需要添加一些曲線成分，比如多項式和對數(shù)變換砂竖、指數(shù)變換等真椿。
同方差性
(1)car包的ncvTest()函數(shù)
ncvTest(yield.fit)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.8439 Df = 1 p = 0.3583
原始假設為誤差方差不變，p = 0.3583無法拒接原假設
(2)car包的spreadLevelPlot()函數(shù)
spreadLevelPlot(yield.fit)
Suggested power transformation: 0.6308
如果圖中的點在水平的最佳擬合曲線周圍呈水平隨機分布乎澄，說明滿足方差不變假設突硝，否則建議冪次轉換為0.5，用根號Y代替Y置济，若建議冪次為0解恰，則使用對數(shù)變換。此例中應使用冪次轉換浙于。
yield.fit2 <- lm(sqrt(yield)~content,data = snake)
summary(yield.fit2)
Call:
lm(formula = sqrt(yield) ~ content, data = snake)
Residuals:
Min 1Q Median 3Q Max
-0.329 -0.150 -0.020 0.146 0.365
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.04727 0.19425 10.5 2.5e-08 ***
content 0.06233 0.00621 10.0 4.8e-08 ***
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.219 on 15 degrees of freedom
Multiple R-squared: 0.87, Adjusted R-squared: 0.862
F-statistic: 101 on 1 and 15 DF, p-value: 4.77e-08
冪次變換后線性模型的擬合效果稍微提高了护盈，當然也可能是由于離群點導致，暫不做分析羞酗。
