预测怀俄明州蛇河流域的水量，数据集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值高度显著，可以拒绝原假设。
回到刚刚的散点图，为散点图添加线性回归模型产生的拟合直线。
with(snake,plot(content,yield))
abline(yield.fit,lwd=3,col="red")
线性回归必须通过统计假设检验。
正态性：对于固定的自变量值，因变量值呈正态分布。
独立性：Y值之间相互独立。
线性：因变量和自变量之间为线性相关。
同方差性：因变量的方差不随自变量的水平不同而变化。
对模型进行回归诊断
par(mfrow=c(2,2))
plot(yield.fit)
标准方法
正态性：右上角QQ图是在正态分布对应的值下，标准化残差的概率图，若满足正态假设，那么图上的点应该落在呈45度角的直线上。
独立性：从收集的数据来验证。
线性：左上角残差与拟合图中，残差值和拟合值不存在任何系统的关联。
同方差性：左下角位置尺度图中，水平线的点应该随机分布。
改进的方法
正态性：
(1)car包qqPlot()函数
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
若图形存在非线性，则说明预测变量的函数形式建模不够充分，需要添加一些曲线成分，比如多项式和对数变换、指数变换等。
同方差性
(1)car包的ncvTest()函数
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()函数
spreadLevelPlot(yield.fit)
Suggested power transformation: 0.6308
如果图中的点在水平的最佳拟合曲线周围呈水平随机分布，说明满足方差不变假设，否则建议幂次转换为0.5，用根号Y代替Y，若建议幂次为0，则使用对数变换。此例中应使用幂次转换。
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
幂次变换后线性模型的拟合效果稍微提高了，当然也可能是由于离群点导致，暂不做分析。
