- 执行多项式回归使用
age预测wage。使用交叉验证为多项式选择最佳次数。选择了什么程度，这与使用ANOVA进行假设检验的结果相比如何？对所得多项式拟合数据进行绘图。
加载工资数据集。保留所有交叉验证误差的数组。我们执行K=10 K倍交叉验证。
1. rm(list = ls()) 2. set.seed(1) 5. # 测试误差 6. cv.MSE <- NA 8. # 循环遍历年龄 9. for (i in 1:15) { 10. glm.fit <- glm(wage ~ poly(age, i), data = Wage) 11. # 我们使用cv.glm的交叉验证并保留cv误差 12. cv.MSE[i] <- cv.glm(Wage, glm.fit, K = 10)$delta[1] 13. } 14. # 检查结果对象 15. cv.MSE
1. ## [1] 1675.837 1601.012 1598.801 1594.217 1594.625 1594.888 1595.500 2. ## [8] 1595.436 1596.335 1595.835 1595.970 1597.971 1598.713 1599.253 3. ## [15] 1595.332
我们通过绘制type = "b"点与线之间的关系图来说明结果。
1. # 用线图说明结果 2. plot( x = 1:15, y = cv.MSE, xlab = "power of age", ylab = "CV error", 3. type = "b", pch = 19, lwd = 2, bty = "n", 4. ylim = c( min(cv.MSE) - sd(cv.MSE), max(cv.MSE) + sd(cv.MSE) ) ) 6. # 水平线 7. abline(h = min(cv.MSE) + sd(cv.MSE) , lty = "dotted") 9. # 最小值 10. points( x = which.min(cv.MSE), y = min(cv.MSE), col = "red", pch = "X", cex = 1.5 )
我们再次以较高的年龄权重对模型进行拟合以进行方差分析。
1. ## Analysis of Deviance Table 2. ## 3. ## Model 1: wage ~ poly(age, a) 4. ## Model 2: wage ~ poly(age, a) 5. ## Model 3: wage ~ poly(age, a) 6. ## Model 4: wage ~ poly(age, a) 7. ## Model 5: wage ~ poly(age, a) 8. ## Model 6: wage ~ poly(age, a) 9. ## Model 7: wage ~ poly(age, a) 10. ## Model 8: wage ~ poly(age, a) 11. ## Model 9: wage ~ poly(age, a) 12. ## Model 10: wage ~ poly(age, a) 13. ## Model 11: wage ~ poly(age, a) 14. ## Model 12: wage ~ poly(age, a) 15. ## Model 13: wage ~ poly(age, a) 16. ## Model 14: wage ~ poly(age, a) 17. ## Model 15: wage ~ poly(age, a) 18. ## Resid. Df Resid. Dev Df Deviance F Pr(>F) 19. ## 1 2998 5022216 20. ## 2 2997 4793430 1 228786 143.5637 < 2.2e-16 *** 21. ## 3 2996 4777674 1 15756 9.8867 0.001681 ** 22. ## 4 2995 4771604 1 6070 3.8090 0.051070 . 23. ## 5 2994 4770322 1 1283 0.8048 0.369731 24. ## 6 2993 4766389 1 3932 2.4675 0.116329 25. ## 7 2992 4763834 1 2555 1.6034 0.205515 26. ## 8 2991 4763707 1 127 0.0795 0.778016 27. ## 9 2990 4756703 1 7004 4.3952 0.036124 * 28. ## 10 2989 4756701 1 3 0.0017 0.967552 29. ## 11 2988 4756597 1 103 0.0648 0.799144 30. ## 12 2987 4756591 1 7 0.0043 0.947923 31. ## 13 2986 4756401 1 190 0.1189 0.730224 32. ## 14 2985 4756158 1 243 0.1522 0.696488 33. ## 15 2984 4755364 1 795 0.4986 0.480151 34. ## --- 35. ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
根据F检验，我们应该选择年龄提高到3次的模型，通过交叉验证 。
现在，我们绘制多项式拟合的结果。
1. plot(wage ~ age, data = Wage, col = "darkgrey", bty = "n") 2. ...
- 拟合step函数以
wage使用进行预测age。绘制获得的拟合图。
1. cv.error <- NA 2. ... 4. # 突出显示最小值 5. points( x = which.min(cv.error), y = min(cv.error, na.rm = TRUE),
交叉验证表明，在k = 8的情况下，测试误差最小。最小值的1sd之内的最简约模型具有k = 4，因此将数据分为5个不同的区域。
44
1. lm.fit <- glm(wage ~ cut(age, 4), data = Wage) 3. matlines(age.grid, cbind( lm.pred$fit + 2* lm.pred$se.fit, 4. lm.pred$fit - 2* lm.pred$se.fit), 5. col = "red", lty ="dashed")
该Wage数据集包含了一些其他的变量，如婚姻状况（maritl），工作级别（jobclass），等等。探索其中一些其他预测变量与的关系wage，并使用非线性拟合技术将模型拟合到数据中。
1. ... 2. summary(Wage[, c("maritl", "jobclass")] )
1. ## maritl jobclass 2. ## 1. Never Married: 648 1. Industrial :1544 3. ## 2. Married :2074 2. Information:1456 4. ## 3. Widowed : 19 5. ## 4. Divorced : 204 6. ## 5. Separated : 55
1. # 关系箱线图 2. par(mfrow=c(1,2)) 3. plot(Wage$maritl, Wage$wage, frame.plot = "FALSE")
看来一对已婚夫妇平均比其他群体挣更多的钱。信息类工作的工资平均高于工业类工作。
多项式和step函数
1. m1 <- lm(wage ~ maritl, data = Wage) 2. deviance(m1) # fit (RSS in linear; -2*logLik)
## [1] 4858941
1. m2 <- lm(wage ~ jobclass, data = Wage) 2. deviance(m2)
## [1] 4998547
1. m3 <- lm(wage ~ maritl + jobclass, data = Wage) 2. deviance(m3)
## [1] 4654752
正如预期的那样，使用最复杂的模型可以使样本内数据拟合最小化。
我们不能使样条曲线拟合分类变量。
我们不能将样条曲线拟合到因子，但可以使用一个样条曲线拟合一个连续变量并添加其他预测变量的模型。
1. library(gam) 2. m4 <- gam(...) 3. deviance(m4)
## [1] 4476501
anova(m1, m2, m3, m4)
1. ## Analysis of Variance Table 2. ## 3. ## Model 1: wage ~ maritl 4. ## Model 2: wage ~ jobclass 5. ## Model 3: wage ~ maritl + jobclass 6. ## Model 4: wage ~ maritl + jobclass + s(age, 4) 7. ## Res.Df RSS Df Sum of Sq F Pr(>F) 8. ## 1 2995 4858941 9. ## 2 2998 4998547 -3.0000 -139606 31.082 < 2.2e-16 *** 10. ## 3 2994 4654752 4.0000 343795 57.408 < 2.2e-16 *** 11. ## 4 2990 4476501 4.0002 178252 29.764 < 2.2e-16 *** 12. ## --- 13. ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
F检验表明，我们从模型四到模型一统计显著改善的变量有年龄，wage，maritl，和jobclass。
