Boston数据回归
这个问题使用的变量dis（到五个波士顿就业中心的距离的加权平均值）和nox（每百万人口中一氧化氮的浓度，单位为百万）。我们将dis作为预测变量，将nox作为因变量。
1. rm(list = ls()) 2. set.seed(1) 3. library(MASS) 4. attach(Boston)
1. ## The following objects are masked from Boston (pos = 14): 2. ## 3. ## age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, 4. ## rad, rm, tax, zn
- 使用
poly()函数拟合三次多项式回归来预测nox使用dis。报告回归输出，并绘制结果数据和多项式拟合。
1. m1 <- lm(nox ~ poly(dis, 3)) 2. summary(m1)
1. ## 2. ## Call: 3. ## lm(formula = nox ~ poly(dis, 3)) 4. ## 5. ## Residuals: 6. ## Min 1Q Median 3Q Max 7. ## -0.121130 -0.040619 -0.009738 0.023385 0.194904 8. ## 9. ## Coefficients: 10. ## Estimate Std. Error t value Pr(>|t|) 11. ## (Intercept) 0.554695 0.002759 201.021 < 2e-16 *** 12. ## poly(dis, 3)1 -2.003096 0.062071 -32.271 < 2e-16 *** 13. ## poly(dis, 3)2 0.856330 0.062071 13.796 < 2e-16 *** 14. ## poly(dis, 3)3 -0.318049 0.062071 -5.124 4.27e-07 *** 15. ## --- 16. ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 17. ## 18. ## Residual standard error: 0.06207 on 502 degrees of freedom 19. ## Multiple R-squared: 0.7148, Adjusted R-squared: 0.7131 20. ## F-statistic: 419.3 on 3 and 502 DF, p-value: < 2.2e-16
1. dislim <- range(dis) 2. ... 3. lines(x = dis.grid, y = lm.pred$fit, col = "red", lwd = 2) 4. matlines(x = dis.grid, y = cbind(lm.pred$fit + 2* lm.pred$se.fit, 5. lm.pred$fit - 2* lm.pred$se.fit) 6. , col = "red", lwd = 1.5, lty = "dashed")
摘要显示，在nox使用进行预测时，所有多项式项都是有效的dis。该图显示了一条平滑的曲线，很好地拟合了数据。
- 绘制多项式适合不同多项式次数的范围（例如，从1到10），并报告相关的残差平方和。
我们绘制1到10度的多项式并保存RSS。
1. # 2. train.rss <- NA 3. ... 4. # 在训练集中显示模型拟合 5. train.rss
1. ## [1] 2.768563 2.035262 1.934107 1.932981 1.915290 1.878257 1.849484 2. ## [8] 1.835630 1.833331 1.832171
正如预期的那样，RSS随多项式次数单调递减。
- 执行交叉验证或其他方法来选择多项式的最佳次数，并解释您的结果。
我们执行LLOCV并手工编码：
1. cv.error <- matrix(NA, nrow = nrow(Boston), ncol = 10) 3. ... 4. names(result) <- paste( "^", 1:10, sep= "" ) 5. result
1. ## ^1 ^2 ^3 ^4 ^5 ^6 2. ## 0.005471468 0.004022257 0.003822345 0.003820121 0.003785158 0.003711971 3. ## ^7 ^8 ^9 ^10 4. ## 0.003655106 0.003627727 0.003623183 0.003620892
1. plot(result ~ seq(1,10), type = "b", pch = 19, bty = "n", xlab = "powers of dist to empl. center", 2. ylab = "cv error") 3. abline(h = min(cv.error) + sd(cv.error), col = "red", lty = "dashed")
基于交叉验证，我们将选择dis平方。
- 使用
bs()函数拟合回归样条曲线使用nox进行预测dis。
我们以[3,6,9]大致相等的4个区间划分此范围
1. library(splines) 2. m4 <- lm(nox ~ bs(dis, knots = c(3, 6, 9))) 3. summary(m4)
1. ## 2. ## Call: 3. ## lm(formula = nox ~ bs(dis, knots = c(3, 6, 9))) 4. ## 5. ## Residuals: 6. ## Min 1Q Median 3Q Max 7. ## -0.132134 -0.039466 -0.009042 0.025344 0.187258 8. ## 9. ## Coefficients: 10. ## Estimate Std. Error t value Pr(>|t|) 11. ## (Intercept) 0.709144 0.016099 44.049 < 2e-16 *** 12. ## bs(dis, knots = c(3, 6, 9))1 0.006631 0.025467 0.260 0.795 13. ## bs(dis, knots = c(3, 6, 9))2 -0.258296 0.017759 -14.544 < 2e-16 *** 14. ## bs(dis, knots = c(3, 6, 9))3 -0.233326 0.027248 -8.563 < 2e-16 *** 15. ## bs(dis, knots = c(3, 6, 9))4 -0.336530 0.032140 -10.471 < 2e-16 *** 16. ## bs(dis, knots = c(3, 6, 9))5 -0.269575 0.058799 -4.585 5.75e-06 *** 17. ## bs(dis, knots = c(3, 6, 9))6 -0.303386 0.062631 -4.844 1.70e-06 *** 18. ## --- 19. ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 20. ## 21. ## Residual standard error: 0.0612 on 499 degrees of freedom 22. ## Multiple R-squared: 0.7244, Adjusted R-squared: 0.7211 23. ## F-statistic: 218.6 on 6 and 499 DF, p-value: < 2.2e-16
1. # 绘图结果 2. ... 3. # 4. matlines( dis.grid, 5. ... 6. col = "black", lwd = 2, lty = c("solid", "dashed", "dashed"))
- 现在针对一定范围的自由度拟合样条回归，并绘制结果拟合并报告结果RSS。描述获得的结果。
我们使用3到16之间的dfs拟合回归样条曲线。
1. box <- NA 3. for (i in 3:16) { 4. ... 5. } 7. box[-c(1, 2)]
1. ## [1] 1.934107 1.922775 1.840173 1.833966 1.829884 1.816995 1.825653 2. ## [8] 1.792535 1.796992 1.788999 1.782350 1.781838 1.782798 1.783546
ISLR包中的College数据集。
- 将数据分为训练集和测试集。使用学费作为因变量，使用其他变量作为预测变量，对训练集执行前向逐步选择，确定仅使用预测变量子集的令人满意的模型。
1. rm(list = ls()) 2. set.seed(1) 3. library(leaps) 4. attach(College)
1. ## The following objects are masked from College (pos = 14): 2. ## 3. ## Accept, Apps, Books, Enroll, Expend, F.Undergrad, Grad.Rate, 4. ## Outstate, P.Undergrad, perc.alumni, Personal, PhD, Private, 5. ## Room.Board, S.F.Ratio, Terminal, Top10perc, Top25perc
1. # 训练/测试拆分行的索引号 2. train <- sample( length(Outstate), length(Outstate)/2) 3. test <- -train 4. ... 5. abline(h=max.adjr2 - std.adjr2, col="red", lty=2)
所有cp，BIC和adjr2得分均显示6是该子集的最小大小。但是，根据1个标准误差规则，我们将选择具有4个预测变量的模型。
1. ... 2. coefi <- coef(m5, id = 4) 3. names(coefi)
## [1] "(Intercept)" "PrivateYes" "Room.Board" "perc.alumni" "Expend"
- 将GAM拟合到训练数据上，使用学费作为响应，并使用在上一步中选择的函数作为预测变量。绘制结果，并解释您的发现。
1. library(gam) 2. ... 3. plot(gam.fit, se=TRUE, col="blue")
- 评估在测试集上获得的模型，并解释获得的结果。
1. gam.pred <- predict(gam.fit, College.test) 2. gam.err <- mean((College.test$Outstate - gam.pred)^2) 3. gam.err
## [1] 3745460
1. gam.tss <- mean((College.test$Outstate - mean(College.test$Outstate))^2) 2. test.rss <- 1 - gam.err / gam.tss 3. test.rss
## [1] 0.7696916
- 对于哪些变量（如果有），是否存在与因变量呈非线性关系的证据？
summary(gam.fit)
1. ## 2. ## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 2) + s(PhD, 3. ## df = 2) + s(perc.alumni, df = 2) + s(Expend, df = 5) + s(Grad.Rate, 4. ## df = 2), data = College.train) 5. ## Deviance Residuals: 6. ## Min 1Q Median 3Q Max 7. ## -4977.74 -1184.52 58.33 1220.04 7688.30 8. ## 9. ## (Dispersion Parameter for gaussian family taken to be 3300711) 10. ## 11. ## Null Deviance: 6221998532 on 387 degrees of freedom 12. ## Residual Deviance: 1231165118 on 373 degrees of freedom 13. ## AIC: 6941.542 14. ## 15. ## Number of Local Scoring Iterations: 2 16. ## 17. ## Anova for Parametric Effects 18. ## Df Sum Sq Mean Sq F value Pr(>F) 19. ## Private 1 1779433688 1779433688 539.106 < 2.2e-16 *** 20. ## s(Room.Board, df = 2) 1 1221825562 1221825562 370.171 < 2.2e-16 *** 21. ## s(PhD, df = 2) 1 382472137 382472137 115.876 < 2.2e-16 *** 22. ## s(perc.alumni, df = 2) 1 328493313 328493313 99.522 < 2.2e-16 *** 23. ## s(Expend, df = 5) 1 416585875 416585875 126.211 < 2.2e-16 *** 24. ## s(Grad.Rate, df = 2) 1 55284580 55284580 16.749 5.232e-05 *** 25. ## Residuals 373 1231165118 3300711 26. ## --- 27. ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 28. ## 29. ## Anova for Nonparametric Effects 30. ## Npar Df Npar F Pr(F) 31. ## (Intercept) 32. ## Private 33. ## s(Room.Board, df = 2) 1 3.5562 0.06010 . 34. ## s(PhD, df = 2) 1 4.3421 0.03786 * 35. ## s(perc.alumni, df = 2) 1 1.9158 0.16715 36. ## s(Expend, df = 5) 4 16.8636 1.016e-12 *** 37. ## s(Grad.Rate, df = 2) 1 3.7208 0.05450 . 38. ## --- 39. ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
非参数Anova检验显示了因变量与支出之间存在非线性关系的有力证据，以及因变量与Grad.Rate或PhD之间具有中等强度的非线性关系（使用p值为0.05）。
