乐趣区

关于r语言:R语言ISLR工资数据进行多项式回归和样条回归分析

jiezi

 原文链接：http://tecdat.cn/?p=8531

  1. 执行多项式回归应用 age 预测wage。应用穿插验证为多项式抉择最佳次数。抉择了什么水平，这与应用进行假设检验的后果相比如何ANOVA？对所得多项式拟合数据进行绘图。

加载工资数据集。保留所有穿插验证谬误的数组。咱们正在执行 K =10  K 倍穿插验证。

rm(list = ls())set.seed(1)library(ISLR)library(boot)# container of test errorscv.MSE <- NA# loop over powers of agefor (i in 1:15) {glm.fit <-  glm(wage ~ poly(age, i), data = Wage)  # we use cv.glm's cross-validation and keep the vanilla cv test error  cv.MSE[i] <-  cv.glm(Wage, glm.fit, K = 10)$delta[1]}# inspect results objectcv.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" 点与线之间的关系图来阐明后果。

# illustrate results with a line plot connecting the cv.error dotsplot(x = 1:15, y = cv.MSE, xlab = "power of age", ylab = "CV error",       type = "b", pch = 19, lwd = 2, bty = "n",       ylim = c( min(cv.MSE) - sd(cv.MSE), max(cv.MSE) + sd(cv.MSE) ) )# horizontal line for 1se to less complexityabline(h = min(cv.MSE) + sd(cv.MSE) , lty = "dotted")# where is the minimumpoints(x = which.min(cv.MSE), y = min(cv.MSE), col = "red", pch = "X", cex = 1.5 )

咱们再次以较高的年龄权重对模型进行拟合以进行方差分析。

`

  1. Analysis of Deviance Table

  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 的幂的模型，而通过穿插验证。

当初，咱们绘制多项式拟合的后果。

plot(wage ~ age, data = Wage, col = "darkgrey",  bty = "n")...

  1. 拟合阶跃函数以 wage 应用进行预测age。绘制取得的拟合图。
cv.error <-  NA...# highlight minimumpoints(x = which.min(cv.error), y = min(cv.error, na.rm = TRUE), col = "red", pch = "X", cex = 1.5 )

k=8 k=81sd1sdk=4k=4

44

lm.fit <-  glm(wage ~ cut(age, 4), data = Wage)...matlines(age.grid, cbind( lm.pred$fit + 2* lm.pred$se.fit,                          lm.pred$fit - 2* lm.pred$se.fit),         col = "red", lty ="dashed")

Q2

该 Wage 数据集蕴含了一些其余的性能，咱们还没有笼罩，如婚姻状况（maritl），作业类（jobclass），等等。摸索其中一些其余预测变量与的关系wage，并应用非线性拟合技术将灵便的模型拟合到数据中。

...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

`

# boxplots of relationshipspar(mfrow=c(1,2))plot(Wage$maritl, Wage$wage, frame.plot = "FALSE")plot(Wage$jobclass, Wage$wage, frame.plot = "FALSE")

看来一对已婚夫妇均匀比其余群体挣更多的钱。看来，信息性工作的工资均匀高于工业性工作。

多项式和阶跃函数

m1 <-  lm(wage ~ maritl, data = Wage)deviance(m1) # fit (RSS in linear; -2*logLik)

## [1] 4858941

m2 <-  lm(wage ~ jobclass, data = Wage)deviance(m2)

## [1] 4998547

m3 <-  lm(wage ~ maritl + jobclass, data = Wage)deviance(m3)

## [1] 4654752

正如预期的那样，应用最简单的模型能够使样本内数据拟合最小化。

咱们不能使样条曲线适宜分类变量。

咱们不能将样条曲线拟合到因子，但能够应用一个样条曲线拟合一个连续变量并增加其余预测变量的模型。

library(gam)m4 <-  gam(...)deviance(m4)

## [1] 4476501

anova(m1, m2, m3, m4)

`

  1. Analysis of Variance Table

  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。

Boston 数据回归

变量 dis（间隔五个波士顿待业核心的加权均匀），nox（在每 10 百万份的氮氧化物浓度）。咱们将其dis 视为预测因素和 nox 作为 响应变量。

rm(list = ls())set.seed(1)library(MASS)attach(Boston)

`

  1. The following objects are masked from Boston (pos = 14):

  3. age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio,

  4. rad, rm, tax, zn

`

  1. 应用 poly() 函数拟合三次多项式回归来预测 nox 应用dis。报告回归输入，并绘制后果数据和多项式拟合。
m1 <-  lm(nox ~ poly(dis, 3))summary(m1)

`

  2. Call:

  3. lm(formula = nox ~ poly(dis, 3))

  5. Residuals:

  6. Min 1Q Median 3Q Max

  7. -0.121130 -0.040619 -0.009738 0.023385 0.194904

  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

  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

`

dislim <-  range(dis)...lines(x = dis.grid, y = lm.pred$fit, col = "red", lwd = 2)matlines(x = dis.grid, y = cbind(lm.pred$fit + 2* lm.pred$se.fit,                                 lm.pred$fit - 2* lm.pred$se.fit)          , col = "red", lwd = 1.5, lty = "dashed")

摘要显示，在 nox 应用进行预测时，所有多项式项都是无效的dis。该图显示了一条平滑的曲线，很好地拟合了数据。

  1. 绘制多项式适宜不同多项式度的范畴（例如，从 1 到 10），并报告相干的残差平方和。

咱们绘制 1 到 10 度的多项式并保留 RSS。

# containertrain.rss <-  NA...# show model fit in training settrain.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 随多项式次数枯燥递加。

  1. 执行穿插验证或其余办法来抉择多项式的最佳次数，并解释您的后果。

咱们执行 LLOCV 并手工编码：

# containercv.error <- matrix(NA, nrow = nrow(Boston), ncol = 10)...names(result) <- paste("^", 1:10, sep= "")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

`

plot(result ~ seq(1,10), type = "b", pch = 19, bty = "n", xlab = "powers of dist to empl. center",     ylab = "cv error")abline(h = min(cv.error) + sd(cv.error), col = "red", lty = "dashed")

基于穿插验证，咱们将抉择 dis 平方。

  1. 应用 bs() 函数拟合回归样条曲线以 nox 应用进行预测dis。应用四个自由度报告适宜度的输入。

[3,6,9][3,6,9]bs()`df`knots

library(splines)m4 <-  lm(nox ~ bs(dis, knots = c(3, 6, 9)))summary(m4)

`

  2. Call:

  3. lm(formula = nox ~ bs(dis, knots = c(3, 6, 9)))

  5. Residuals:

  6. Min 1Q Median 3Q Max

  7. -0.132134 -0.039466 -0.009042 0.025344 0.187258

  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

  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

`

# plot results...# all lines at oncematlines(dis.grid,...          col = "black", lwd = 2, lty = c("solid", "dashed", "dashed"))

dis>9dis>9

  1. 当初针对肯定范畴的自由度拟合样条回归，并绘制后果拟合并报告后果 RSS。形容取得的后果。

咱们应用 3 到 16 之间的 dfs 拟合回归样条曲线。

box <-  NAfor (i in 3:16) {...}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

`

df=14df=14

ISLR包中的 College 数据集。

  1. 将数据分为训练集和测试集。应用学费作为响应，应用其余变量作为预测变量，对训练集执行前向逐渐抉择，以便确定仅应用预测变量子集的令人满意的模型。
rm(list = ls())set.seed(1)library(leaps)attach(College)

`

  1. The following objects are masked from College (pos = 14):

  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

`

# train/test split row index numberstrain <-  sample(length(Outstate), length(Outstate)/2)test <-  -train...abline(h=max.adjr2 - std.adjr2, col="red", lty=2)

所有 cp，BIC 和 adjr2 得分均显示大小 6 是该子集的最小大小。然而，依据 1 个标准误差规定，咱们将抉择具备 4 个预测变量的模型。

...coefi <-  coef(m5, id = 4)names(coefi)

## [1] "(Intercept)" "PrivateYes" "Room.Board" "perc.alumni" "Expend"

  1. 将 GAM 拟合到训练数据上，应用州外学费作为响应，并应用在上一步中抉择的性能作为预测变量。绘制后果，并解释您的发现。
library(gam)...plot(gam.fit, se=TRUE, col="blue")

  1. 评估在测试集上取得的模型，并解释取得的后果。
gam.pred <-  predict(gam.fit, College.test)gam.err <-  mean((College.test$Outstate - gam.pred)^2)gam.err

## [1] 3745460

gam.tss <-  mean((College.test$Outstate - mean(College.test$Outstate))^2)test.rss <-  1 - gam.err / gam.tsstest.rss

## [1] 0.7696916

0.770.770.740.74

  1. 对于哪些变量（如果有），是否存在与响应呈非线性关系的证据？
summary(gam.fit)

`

  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

  9. (Dispersion Parameter for gaussian family taken to be 3300711)

  11. Null Deviance: 6221998532 on 387 degrees of freedom

  12. Residual Deviance: 1231165118 on 373 degrees of freedom

  13. AIC: 6941.542

  15. Number of Local Scoring Iterations: 2

  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

  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）。

退出移动版