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

?

1. 執(zhí)行多項式回歸使用age預(yù)測wage。使用交叉驗證為多項式選擇最佳次數(shù)。選擇了什么程度，這與使用進(jìn)行假設(shè)檢驗的結(jié)果相比如何ANOVA？對所得多項式擬合數(shù)據(jù)進(jìn)行繪圖。

加載工資數(shù)據(jù)集。保留所有交叉驗證錯誤的數(shù)組。我們正在執(zhí)行K=10? K倍交叉驗證。

1. rm(list = ls())

2. set.seed(1)

3. library(ISLR)

4. library(boot)

5. 


6. # container of test errors

7. cv.MSE <- NA

8. 


9. # loop over powers of age

10. for (i in 1:15) {

11. glm.fit <- ?glm(wage ~ poly(age, i), data = Wage)

12. # we use cv.glm's cross-validation and keep the vanilla cv test error

13. cv.MSE[i] <- ?cv.glm(Wage, glm.fit, K = 10)$delta[1]

14. }

15. # inspect results object

16. 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"點與線之間的關(guān)系圖來說明結(jié)果。?

1. # illustrate results with a line plot connecting the cv.error dots

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

5. 


6. # horizontal line for 1se to less complexity

7. abline(h = min(cv.MSE) + sd(cv.MSE) , lty = "dotted")

8. 


9. # where is the minimum

10. points( x = which.min(cv.MSE), y = min(cv.MSE), col = "red", pch = "X", cex = 1.5 )

我們再次以較高的年齡權(quán)重對模型進(jìn)行擬合以進(jìn)行方差分析。

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

根據(jù)F檢驗，我們應(yīng)該選擇年齡提高到3的冪的模型，而通過交叉驗證 。

現(xiàn)在，我們繪制多項式擬合的結(jié)果。

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

2. ...

1. 擬合階躍函數(shù)以wage使用進(jìn)行預(yù)測age?。繪制獲得的擬合圖。

1. cv.error <- ?NA

2. ...

3. 


4. # highlight minimum

5. points( 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

1. lm.fit <- ?glm(wage ~ cut(age, 4), data = Wage)

2. ...

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")

?

?

Q2

該Wage數(shù)據(jù)集包含了一些其他的功能，我們還沒有覆蓋，如婚姻狀況（maritl），作業(yè)類（jobclass），等等。探索其中一些其他預(yù)測變量與的關(guān)系wage，并使用非線性擬合技術(shù)將靈活的模型擬合到數(shù)據(jù)中。?

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. # boxplots of relationships

2. par(mfrow=c(1,2))

3. plot(Wage$maritl, Wage$wage, frame.plot = "FALSE")

4. plot(Wage$jobclass, Wage$wage, frame.plot = "FALSE")

?

看來一對已婚夫婦平均比其他群體掙更多的錢?？磥?，信息性工作的工資平均高于工業(yè)性工作。

多項式和階躍函數(shù)

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

正如預(yù)期的那樣，使用最復(fù)雜的模型可以使樣本內(nèi)數(shù)據(jù)擬合最小化。

?

我們不能使樣條曲線適合分類變量。

?

我們不能將樣條曲線擬合到因子，但可以使用一個樣條曲線擬合一個連續(xù)變量并添加其他預(yù)測變量的模型。

1. library(gam)

2. m4 <- ?gam(...)

3. deviance(m4)

## [1] 4476501anova(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檢驗表明，我們從模型四送一統(tǒng)計顯著改善計有年齡花，wage，maritl，和jobclass。

?

?

Boston數(shù)據(jù)回歸

?變量dis（距離五個波士頓就業(yè)中心的加權(quán)平均），?nox?（在每10百萬份的氮氧化物濃度） 。我們將其dis視為預(yù)測因素和nox作為響應(yīng)變量。

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

1. 使用poly()函數(shù)擬合三次多項式回歸來預(yù)測nox使用dis。報告回歸輸出，并繪制結(jié)果數(shù)據(jù)和多項式擬合。

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使用進(jìn)行預(yù)測時，所有多項式項都是有效的dis。該圖顯示了一條平滑的曲線，很好地擬合了數(shù)據(jù)。

1. 繪制多項式適合不同多項式度的范圍（例如，從1到10），并報告相關(guān)的殘差平方和。

我們繪制1到10度的多項式并保存RSS。

1. # container

2. train.rss <- ?NA

3. ...

4. # show model fit in training set

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

正如預(yù)期的那樣，RSS隨多項式次數(shù)單調(diào)遞減。

1. 執(zhí)行交叉驗證或其他方法來選擇多項式的最佳次數(shù)，并解釋您的結(jié)果。

我們執(zhí)行LLOCV并手工編碼：

1. # container

2. cv.error <- matrix(NA, nrow = nrow(Boston), ncol = 10)

3. 


4. ...

5. names(result) <- paste( "^", 1:10, sep= "" )

6. 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平方。

1. 使用bs()函數(shù)擬合回歸樣條曲線以nox使用進(jìn)行預(yù)測dis。使用四個自由度報告適合度的輸出。

[3,6,9][3,6,9]bs()dfknots

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. # plot results

2. ...

3. # all lines at once

4. matlines( dis.grid,

5. ...

6. col = "black", lwd = 2, lty = c("solid", "dashed", "dashed"))

dis>9dis>9

1. 現(xiàn)在針對一定范圍的自由度擬合樣條回歸，并繪制結(jié)果擬合并報告結(jié)果RSS。描述獲得的結(jié)果。

我們使用3到16之間的dfs擬合回歸樣條曲線。

1. box <- ?NA

2. 


3. for (i in 3:16) {

4. ...

5. }

6. 


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

df=14df=14

?

?

ISLR包中的College數(shù)據(jù)集。

1. 將數(shù)據(jù)分為訓(xùn)練集和測試集。使用學(xué)費作為響應(yīng)，使用其他變量作為預(yù)測變量，對訓(xùn)練集執(zhí)行前向逐步選擇，以便確定僅使用預(yù)測變量子集的令人滿意的模型。

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. # train/test split row index numbers

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是該子集的最小大小。但是，根據(jù)1個標(biāo)準(zhǔn)誤差規(guī)則，我們將選擇具有4個預(yù)測變量的模型。

1. ...

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

3. names(coefi)

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

1. 將GAM擬合到訓(xùn)練數(shù)據(jù)上，使用州外學(xué)費作為響應(yīng)，并使用在上一步中選擇的功能作為預(yù)測變量。繪制結(jié)果，并解釋您的發(fā)現(xiàn)。

1. library(gam)

2. ...

3. plot(gam.fit, se=TRUE, col="blue")

?

1. 評估在測試集上獲得的模型，并解釋獲得的結(jié)果。

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

0.770.770.740.74

1. 對于哪些變量（如果有），是否存在與響應(yīng)呈非線性關(guān)系的證據(jù)？

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

非參數(shù)Anova檢驗顯示了響應(yīng)與支出之間存在非線性關(guān)系的有力證據(jù)，以及響應(yīng)與Grad.Rate或PhD之間具有中等強(qiáng)度的非線性關(guān)系（使用p值為0.05）。