【独家发布】An Introduction to Statistical Learning with Applications in R

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

An Introduction to Statistical Learning

with Applications in R

Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani

http://www-bcf.usc.edu/~gareth/ISL/code.html

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏3 回帖

关键词：introduction Applications Application Statistical statistica Robert

1. # Chapter 2 Lab: Introduction to R
   
2. 
   
3. # Basic Commands
   
4. 
   
5. x <- c(1,3,2,5)
   
6. x
   
7. x = c(1,6,2)
   
8. x
   
9. y = c(1,4,3)
   
10. length(x)
    
11. length(y)
    
12. x+y
    
13. ls()
    
14. rm(x,y)
    
15. ls()
    
16. rm(list=ls())
    
17. ?matrix
    
18. x=matrix(data=c(1,2,3,4), nrow=2, ncol=2)
    
19. x
    
20. x=matrix(c(1,2,3,4),2,2)
    
21. matrix(c(1,2,3,4),2,2,byrow=TRUE)
    
22. sqrt(x)
    
23. x^2
    
24. x=rnorm(50)
    
25. y=x+rnorm(50,mean=50,sd=.1)
    
26. cor(x,y)
    
27. set.seed(1303)
    
28. rnorm(50)
    
29. set.seed(3)
    
30. y=rnorm(100)
    
31. mean(y)
    
32. var(y)
    
33. sqrt(var(y))
    
34. sd(y)
    
35. 
    
36. # Graphics
    
37. 
    
38. x=rnorm(100)
    
39. y=rnorm(100)
    
40. plot(x,y)
    
41. plot(x,y,xlab="this is the x-axis",ylab="this is the y-axis",main="Plot of X vs Y")
    
42. pdf("Figure.pdf")
    
43. plot(x,y,col="green")
    
44. dev.off()
    
45. x=seq(1,10)
    
46. x
    
47. x=1:10
    
48. x
    
49. x=seq(-pi,pi,length=50)
    
50. y=x
    
51. f=outer(x,y,function(x,y)cos(y)/(1+x^2))
    
52. contour(x,y,f)
    
53. contour(x,y,f,nlevels=45,add=T)
    
54. fa=(f-t(f))/2
    
55. contour(x,y,fa,nlevels=15)
    
56. image(x,y,fa)
    
57. persp(x,y,fa)
    
58. persp(x,y,fa,theta=30)
    
59. persp(x,y,fa,theta=30,phi=20)
    
60. persp(x,y,fa,theta=30,phi=70)
    
61. persp(x,y,fa,theta=30,phi=40)
    
62. 
    
63. # Indexing Data
    
64. 
    
65. A=matrix(1:16,4,4)
    
66. A
    
67. A[2,3]
    
68. A[c(1,3),c(2,4)]
    
69. A[1:3,2:4]
    
70. A[1:2,]
    
71. A[,1:2]
    
72. A[1,]
    
73. A[-c(1,3),]
    
74. A[-c(1,3),-c(1,3,4)]
    
75. dim(A)
    
76. 
    
77. # Loading Data
    
78. 
    
79. Auto=read.table("Auto.data")
    
80. fix(Auto)
    
81. Auto=read.table("Auto.data",header=T,na.strings="?")
    
82. fix(Auto)
    
83. Auto=read.csv("Auto.csv",header=T,na.strings="?")
    
84. fix(Auto)
    
85. dim(Auto)
    
86. Auto[1:4,]
    
87. Auto=na.omit(Auto)
    
88. dim(Auto)
    
89. names(Auto)
    
90. 
    
91. # Additional Graphical and Numerical Summaries
    
92. 
    
93. plot(cylinders, mpg)
    
94. plot(Auto$cylinders, Auto$mpg)
    
95. attach(Auto)
    
96. plot(cylinders, mpg)
    
97. cylinders=as.factor(cylinders)
    
98. plot(cylinders, mpg)
    
99. plot(cylinders, mpg, col="red")
    
100. plot(cylinders, mpg, col="red", varwidth=T)
     
101. plot(cylinders, mpg, col="red", varwidth=T,horizontal=T)
     
102. plot(cylinders, mpg, col="red", varwidth=T, xlab="cylinders", ylab="MPG")
     
103. hist(mpg)
     
104. hist(mpg,col=2)
     
105. hist(mpg,col=2,breaks=15)
     
106. pairs(Auto)
     
107. pairs(~ mpg + displacement + horsepower + weight + acceleration, Auto)
     
108. plot(horsepower,mpg)
     
109. identify(horsepower,mpg,name)
     
110. summary(Auto)
     
111. summary(mpg)

复制代码

1. # Chapter 3 Lab: Linear Regression
   
2. 
   
3. library(MASS)
   
4. library(ISLR)
   
5. 
   
6. # Simple Linear Regression
   
7. 
   
8. fix(Boston)
   
9. names(Boston)
   
10. lm.fit=lm(medv~lstat)
    
11. lm.fit=lm(medv~lstat,data=Boston)
    
12. attach(Boston)
    
13. lm.fit=lm(medv~lstat)
    
14. lm.fit
    
15. summary(lm.fit)
    
16. names(lm.fit)
    
17. coef(lm.fit)
    
18. confint(lm.fit)
    
19. predict(lm.fit,data.frame(lstat=(c(5,10,15))), interval="confidence")
    
20. predict(lm.fit,data.frame(lstat=(c(5,10,15))), interval="prediction")
    
21. plot(lstat,medv)
    
22. abline(lm.fit)
    
23. abline(lm.fit,lwd=3)
    
24. abline(lm.fit,lwd=3,col="red")
    
25. plot(lstat,medv,col="red")
    
26. plot(lstat,medv,pch=20)
    
27. plot(lstat,medv,pch="+")
    
28. plot(1:20,1:20,pch=1:20)
    
29. par(mfrow=c(2,2))
    
30. plot(lm.fit)
    
31. plot(predict(lm.fit), residuals(lm.fit))
    
32. plot(predict(lm.fit), rstudent(lm.fit))
    
33. plot(hatvalues(lm.fit))
    
34. which.max(hatvalues(lm.fit))
    
35. 
    
36. # Multiple Linear Regression
    
37. 
    
38. lm.fit=lm(medv~lstat+age,data=Boston)
    
39. summary(lm.fit)
    
40. lm.fit=lm(medv~.,data=Boston)
    
41. summary(lm.fit)
    
42. library(car)
    
43. vif(lm.fit)
    
44. lm.fit1=lm(medv~.-age,data=Boston)
    
45. summary(lm.fit1)
    
46. lm.fit1=update(lm.fit, ~.-age)
    
47. 
    
48. # Interaction Terms
    
49. 
    
50. summary(lm(medv~lstat*age,data=Boston))
    
51. 
    
52. # Non-linear Transformations of the Predictors
    
53. 
    
54. lm.fit2=lm(medv~lstat+I(lstat^2))
    
55. summary(lm.fit2)
    
56. lm.fit=lm(medv~lstat)
    
57. anova(lm.fit,lm.fit2)
    
58. par(mfrow=c(2,2))
    
59. plot(lm.fit2)
    
60. lm.fit5=lm(medv~poly(lstat,5))
    
61. summary(lm.fit5)
    
62. summary(lm(medv~log(rm),data=Boston))
    
63. 
    
64. # Qualitative Predictors
    
65. 
    
66. fix(Carseats)
    
67. names(Carseats)
    
68. lm.fit=lm(Sales~.+Income:Advertising+Price:Age,data=Carseats)
    
69. summary(lm.fit)
    
70. attach(Carseats)
    
71. contrasts(ShelveLoc)
    
72. 
    
73. # Writing Functions
    
74. 
    
75. LoadLibraries
    
76. LoadLibraries()
    
77. LoadLibraries=function(){
    
78. library(ISLR)
    
79. library(MASS)
    
80. print("The libraries have been loaded.")
    
81. }
    
82. LoadLibraries
    
83. LoadLibraries()

复制代码

1. # Chapter 4 Lab: Logistic Regression, LDA, QDA, and KNN
   
2. 
   
3. # The Stock Market Data
   
4. 
   
5. library(ISLR)
   
6. names(Smarket)
   
7. dim(Smarket)
   
8. summary(Smarket)
   
9. pairs(Smarket)
   
10. cor(Smarket)
    
11. cor(Smarket[,-9])
    
12. attach(Smarket)
    
13. plot(Volume)
    
14. 
    
15. # Logistic Regression
    
16. 
    
17. glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomial)
    
18. summary(glm.fit)
    
19. coef(glm.fit)
    
20. summary(glm.fit)$coef
    
21. summary(glm.fit)$coef[,4]
    
22. glm.probs=predict(glm.fit,type="response")
    
23. glm.probs[1:10]
    
24. contrasts(Direction)
    
25. glm.pred=rep("Down",1250)
    
26. glm.pred[glm.probs>.5]="Up"
    
27. table(glm.pred,Direction)
    
28. (507+145)/1250
    
29. mean(glm.pred==Direction)
    
30. train=(Year<2005)
    
31. Smarket.2005=Smarket[!train,]
    
32. dim(Smarket.2005)
    
33. Direction.2005=Direction[!train]
    
34. glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomial,subset=train)
    
35. glm.probs=predict(glm.fit,Smarket.2005,type="response")
    
36. glm.pred=rep("Down",252)
    
37. glm.pred[glm.probs>.5]="Up"
    
38. table(glm.pred,Direction.2005)
    
39. mean(glm.pred==Direction.2005)
    
40. mean(glm.pred!=Direction.2005)
    
41. glm.fit=glm(Direction~Lag1+Lag2,data=Smarket,family=binomial,subset=train)
    
42. glm.probs=predict(glm.fit,Smarket.2005,type="response")
    
43. glm.pred=rep("Down",252)
    
44. glm.pred[glm.probs>.5]="Up"
    
45. table(glm.pred,Direction.2005)
    
46. mean(glm.pred==Direction.2005)
    
47. 106/(106+76)
    
48. predict(glm.fit,newdata=data.frame(Lag1=c(1.2,1.5),Lag2=c(1.1,-0.8)),type="response")
    
49. 
    
50. # Linear Discriminant Analysis
    
51. 
    
52. library(MASS)
    
53. lda.fit=lda(Direction~Lag1+Lag2,data=Smarket,subset=train)
    
54. lda.fit
    
55. plot(lda.fit)
    
56. lda.pred=predict(lda.fit, Smarket.2005)
    
57. names(lda.pred)
    
58. lda.class=lda.pred$class
    
59. table(lda.class,Direction.2005)
    
60. mean(lda.class==Direction.2005)
    
61. sum(lda.pred$posterior[,1]>=.5)
    
62. sum(lda.pred$posterior[,1]<.5)
    
63. lda.pred$posterior[1:20,1]
    
64. lda.class[1:20]
    
65. sum(lda.pred$posterior[,1]>.9)
    
66. 
    
67. # Quadratic Discriminant Analysis
    
68. 
    
69. qda.fit=qda(Direction~Lag1+Lag2,data=Smarket,subset=train)
    
70. qda.fit
    
71. qda.class=predict(qda.fit,Smarket.2005)$class
    
72. table(qda.class,Direction.2005)
    
73. mean(qda.class==Direction.2005)
    
74. 
    
75. # K-Nearest Neighbors
    
76. 
    
77. library(class)
    
78. train.X=cbind(Lag1,Lag2)[train,]
    
79. test.X=cbind(Lag1,Lag2)[!train,]
    
80. train.Direction=Direction[train]
    
81. set.seed(1)
    
82. knn.pred=knn(train.X,test.X,train.Direction,k=1)
    
83. table(knn.pred,Direction.2005)
    
84. (83+43)/252
    
85. knn.pred=knn(train.X,test.X,train.Direction,k=3)
    
86. table(knn.pred,Direction.2005)
    
87. mean(knn.pred==Direction.2005)
    
88. 
    
89. # An Application to Caravan Insurance Data
    
90. 
    
91. dim(Caravan)
    
92. attach(Caravan)
    
93. summary(Purchase)
    
94. 348/5822
    
95. standardized.X=scale(Caravan[,-86])
    
96. var(Caravan[,1])
    
97. var(Caravan[,2])
    
98. var(standardized.X[,1])
    
99. var(standardized.X[,2])
    
100. test=1:1000
     
101. train.X=standardized.X[-test,]
     
102. test.X=standardized.X[test,]
     
103. train.Y=Purchase[-test]
     
104. test.Y=Purchase[test]
     
105. set.seed(1)
     
106. knn.pred=knn(train.X,test.X,train.Y,k=1)
     
107. mean(test.Y!=knn.pred)
     
108. mean(test.Y!="No")
     
109. table(knn.pred,test.Y)
     
110. 9/(68+9)
     
111. knn.pred=knn(train.X,test.X,train.Y,k=3)
     
112. table(knn.pred,test.Y)
     
113. 5/26
     
114. knn.pred=knn(train.X,test.X,train.Y,k=5)
     
115. table(knn.pred,test.Y)
     
116. 4/15
     
117. glm.fit=glm(Purchase~.,data=Caravan,family=binomial,subset=-test)
     
118. glm.probs=predict(glm.fit,Caravan[test,],type="response")
     
119. glm.pred=rep("No",1000)
     
120. glm.pred[glm.probs>.5]="Yes"
     
121. table(glm.pred,test.Y)
     
122. glm.pred=rep("No",1000)
     
123. glm.pred[glm.probs>.25]="Yes"
     
124. table(glm.pred,test.Y)
     
125. 11/(22+11)

复制代码

1. # Chaper 5 Lab: Cross-Validation and the Bootstrap
   
2. 
   
3. # The Validation Set Approach
   
4. 
   
5. library(ISLR)
   
6. set.seed(1)
   
7. train=sample(392,196)
   
8. lm.fit=lm(mpg~horsepower,data=Auto,subset=train)
   
9. attach(Auto)
   
10. mean((mpg-predict(lm.fit,Auto))[-train]^2)
    
11. lm.fit2=lm(mpg~poly(horsepower,2),data=Auto,subset=train)
    
12. mean((mpg-predict(lm.fit2,Auto))[-train]^2)
    
13. lm.fit3=lm(mpg~poly(horsepower,3),data=Auto,subset=train)
    
14. mean((mpg-predict(lm.fit3,Auto))[-train]^2)
    
15. set.seed(2)
    
16. train=sample(392,196)
    
17. lm.fit=lm(mpg~horsepower,subset=train)
    
18. mean((mpg-predict(lm.fit,Auto))[-train]^2)
    
19. lm.fit2=lm(mpg~poly(horsepower,2),data=Auto,subset=train)
    
20. mean((mpg-predict(lm.fit2,Auto))[-train]^2)
    
21. lm.fit3=lm(mpg~poly(horsepower,3),data=Auto,subset=train)
    
22. mean((mpg-predict(lm.fit3,Auto))[-train]^2)
    
23. 
    
24. # Leave-One-Out Cross-Validation
    
25. 
    
26. glm.fit=glm(mpg~horsepower,data=Auto)
    
27. coef(glm.fit)
    
28. lm.fit=lm(mpg~horsepower,data=Auto)
    
29. coef(lm.fit)
    
30. library(boot)
    
31. glm.fit=glm(mpg~horsepower,data=Auto)
    
32. cv.err=cv.glm(Auto,glm.fit)
    
33. cv.err$delta
    
34. cv.error=rep(0,5)
    
35. for (i in 1:5){
    
36. glm.fit=glm(mpg~poly(horsepower,i),data=Auto)
    
37. cv.error[i]=cv.glm(Auto,glm.fit)$delta[1]
    
38. }
    
39. cv.error
    
40. 
    
41. # k-Fold Cross-Validation
    
42. 
    
43. set.seed(17)
    
44. cv.error.10=rep(0,10)
    
45. for (i in 1:10){
    
46. glm.fit=glm(mpg~poly(horsepower,i),data=Auto)
    
47. cv.error.10[i]=cv.glm(Auto,glm.fit,K=10)$delta[1]
    
48. }
    
49. cv.error.10
    
50. 
    
51. # The Bootstrap
    
52. 
    
53. alpha.fn=function(data,index){
    
54. X=data$X[index]
    
55. Y=data$Y[index]
    
56. return((var(Y)-cov(X,Y))/(var(X)+var(Y)-2*cov(X,Y)))
    
57. }
    
58. alpha.fn(Portfolio,1:100)
    
59. set.seed(1)
    
60. alpha.fn(Portfolio,sample(100,100,replace=T))
    
61. boot(Portfolio,alpha.fn,R=1000)
    
62. 
    
63. # Estimating the Accuracy of a Linear Regression Model
    
64. 
    
65. boot.fn=function(data,index)
    
66. return(coef(lm(mpg~horsepower,data=data,subset=index)))
    
67. boot.fn(Auto,1:392)
    
68. set.seed(1)
    
69. boot.fn(Auto,sample(392,392,replace=T))
    
70. boot.fn(Auto,sample(392,392,replace=T))
    
71. boot(Auto,boot.fn,1000)
    
72. summary(lm(mpg~horsepower,data=Auto))$coef
    
73. boot.fn=function(data,index)
    
74. coefficients(lm(mpg~horsepower+I(horsepower^2),data=data,subset=index))
    
75. set.seed(1)
    
76. boot(Auto,boot.fn,1000)
    
77. summary(lm(mpg~horsepower+I(horsepower^2),data=Auto))$coef

复制代码

1. # Chapter 6 Lab 1: Subset Selection Methods
   
2. 
   
3. # Best Subset Selection
   
4. 
   
5. library(ISLR)
   
6. fix(Hitters)
   
7. names(Hitters)
   
8. dim(Hitters)
   
9. sum(is.na(Hitters$Salary))
   
10. Hitters=na.omit(Hitters)
    
11. dim(Hitters)
    
12. sum(is.na(Hitters))
    
13. library(leaps)
    
14. regfit.full=regsubsets(Salary~.,Hitters)
    
15. summary(regfit.full)
    
16. regfit.full=regsubsets(Salary~.,data=Hitters,nvmax=19)
    
17. reg.summary=summary(regfit.full)
    
18. names(reg.summary)
    
19. reg.summary$rsq
    
20. par(mfrow=c(2,2))
    
21. plot(reg.summary$rss,xlab="Number of Variables",ylab="RSS",type="l")
    
22. plot(reg.summary$adjr2,xlab="Number of Variables",ylab="Adjusted RSq",type="l")
    
23. which.max(reg.summary$adjr2)
    
24. points(11,reg.summary$adjr2[11], col="red",cex=2,pch=20)
    
25. plot(reg.summary$cp,xlab="Number of Variables",ylab="Cp",type='l')
    
26. which.min(reg.summary$cp)
    
27. points(10,reg.summary$cp[10],col="red",cex=2,pch=20)
    
28. which.min(reg.summary$bic)
    
29. plot(reg.summary$bic,xlab="Number of Variables",ylab="BIC",type='l')
    
30. points(6,reg.summary$bic[6],col="red",cex=2,pch=20)
    
31. plot(regfit.full,scale="r2")
    
32. plot(regfit.full,scale="adjr2")
    
33. plot(regfit.full,scale="Cp")
    
34. plot(regfit.full,scale="bic")
    
35. coef(regfit.full,6)
    
36. 
    
37. # Forward and Backward Stepwise Selection
    
38. 
    
39. regfit.fwd=regsubsets(Salary~.,data=Hitters,nvmax=19,method="forward")
    
40. summary(regfit.fwd)
    
41. regfit.bwd=regsubsets(Salary~.,data=Hitters,nvmax=19,method="backward")
    
42. summary(regfit.bwd)
    
43. coef(regfit.full,7)
    
44. coef(regfit.fwd,7)
    
45. coef(regfit.bwd,7)
    
46. 
    
47. # Choosing Among Models
    
48. 
    
49. set.seed(1)
    
50. train=sample(c(TRUE,FALSE), nrow(Hitters),rep=TRUE)
    
51. test=(!train)
    
52. regfit.best=regsubsets(Salary~.,data=Hitters[train,],nvmax=19)
    
53. test.mat=model.matrix(Salary~.,data=Hitters[test,])
    
54. val.errors=rep(NA,19)
    
55. for(i in 1:19){
    
56.    coefi=coef(regfit.best,id=i)
    
57.    pred=test.mat[,names(coefi)]%*%coefi
    
58.    val.errors[i]=mean((Hitters$Salary[test]-pred)^2)
    
59. }
    
60. val.errors
    
61. which.min(val.errors)
    
62. coef(regfit.best,10)
    
63. predict.regsubsets=function(object,newdata,id,...){
    
64.   form=as.formula(object$call[[2]])
    
65.   mat=model.matrix(form,newdata)
    
66.   coefi=coef(object,id=id)
    
67.   xvars=names(coefi)
    
68.   mat[,xvars]%*%coefi
    
69.   }
    
70. regfit.best=regsubsets(Salary~.,data=Hitters,nvmax=19)
    
71. coef(regfit.best,10)
    
72. k=10
    
73. set.seed(1)
    
74. folds=sample(1:k,nrow(Hitters),replace=TRUE)
    
75. cv.errors=matrix(NA,k,19, dimnames=list(NULL, paste(1:19)))
    
76. for(j in 1:k){
    
77.   best.fit=regsubsets(Salary~.,data=Hitters[folds!=j,],nvmax=19)
    
78.   for(i in 1:19){
    
79.     pred=predict(best.fit,Hitters[folds==j,],id=i)
    
80.     cv.errors[j,i]=mean( (Hitters$Salary[folds==j]-pred)^2)
    
81.     }
    
82.   }
    
83. mean.cv.errors=apply(cv.errors,2,mean)
    
84. mean.cv.errors
    
85. par(mfrow=c(1,1))
    
86. plot(mean.cv.errors,type='b')
    
87. reg.best=regsubsets(Salary~.,data=Hitters, nvmax=19)
    
88. coef(reg.best,11)
    
89. 
    
90. 
    
91. # Chapter 6 Lab 2: Ridge Regression and the Lasso
    
92. 
    
93. x=model.matrix(Salary~.,Hitters)[,-1]
    
94. y=Hitters$Salary
    
95. 
    
96. # Ridge Regression
    
97. 
    
98. library(glmnet)
    
99. grid=10^seq(10,-2,length=100)
    
100. ridge.mod=glmnet(x,y,alpha=0,lambda=grid)
     
101. dim(coef(ridge.mod))
     
102. ridge.mod$lambda[50]
     
103. coef(ridge.mod)[,50]
     
104. sqrt(sum(coef(ridge.mod)[-1,50]^2))
     
105. ridge.mod$lambda[60]
     
106. coef(ridge.mod)[,60]
     
107. sqrt(sum(coef(ridge.mod)[-1,60]^2))
     
108. predict(ridge.mod,s=50,type="coefficients")[1:20,]
     
109. set.seed(1)
     
110. train=sample(1:nrow(x), nrow(x)/2)
     
111. test=(-train)
     
112. y.test=y[test]
     
113. ridge.mod=glmnet(x[train,],y[train],alpha=0,lambda=grid, thresh=1e-12)
     
114. ridge.pred=predict(ridge.mod,s=4,newx=x[test,])
     
115. mean((ridge.pred-y.test)^2)
     
116. mean((mean(y[train])-y.test)^2)
     
117. ridge.pred=predict(ridge.mod,s=1e10,newx=x[test,])
     
118. mean((ridge.pred-y.test)^2)
     
119. ridge.pred=predict(ridge.mod,s=0,newx=x[test,],exact=T)
     
120. mean((ridge.pred-y.test)^2)
     
121. lm(y~x, subset=train)
     
122. predict(ridge.mod,s=0,exact=T,type="coefficients")[1:20,]
     
123. set.seed(1)
     
124. cv.out=cv.glmnet(x[train,],y[train],alpha=0)
     
125. plot(cv.out)
     
126. bestlam=cv.out$lambda.min
     
127. bestlam
     
128. ridge.pred=predict(ridge.mod,s=bestlam,newx=x[test,])
     
129. mean((ridge.pred-y.test)^2)
     
130. out=glmnet(x,y,alpha=0)
     
131. predict(out,type="coefficients",s=bestlam)[1:20,]
     
132. 
     
133. # The Lasso
     
134. 
     
135. lasso.mod=glmnet(x[train,],y[train],alpha=1,lambda=grid)
     
136. plot(lasso.mod)
     
137. set.seed(1)
     
138. cv.out=cv.glmnet(x[train,],y[train],alpha=1)
     
139. plot(cv.out)
     
140. bestlam=cv.out$lambda.min
     
141. lasso.pred=predict(lasso.mod,s=bestlam,newx=x[test,])
     
142. mean((lasso.pred-y.test)^2)
     
143. out=glmnet(x,y,alpha=1,lambda=grid)
     
144. lasso.coef=predict(out,type="coefficients",s=bestlam)[1:20,]
     
145. lasso.coef
     
146. lasso.coef[lasso.coef!=0]
     
147. 
     
148. 
     
149. # Chapter 6 Lab 3: PCR and PLS Regression
     
150. 
     
151. # Principal Components Regression
     
152. 
     
153. library(pls)
     
154. set.seed(2)
     
155. pcr.fit=pcr(Salary~., data=Hitters,scale=TRUE,validation="CV")
     
156. summary(pcr.fit)
     
157. validationplot(pcr.fit,val.type="MSEP")
     
158. set.seed(1)
     
159. pcr.fit=pcr(Salary~., data=Hitters,subset=train,scale=TRUE, validation="CV")
     
160. validationplot(pcr.fit,val.type="MSEP")
     
161. pcr.pred=predict(pcr.fit,x[test,],ncomp=7)
     
162. mean((pcr.pred-y.test)^2)
     
163. pcr.fit=pcr(y~x,scale=TRUE,ncomp=7)
     
164. summary(pcr.fit)
     
165. 
     
166. # Partial Least Squares
     
167. 
     
168. set.seed(1)
     
169. pls.fit=plsr(Salary~., data=Hitters,subset=train,scale=TRUE, validation="CV")
     
170. summary(pls.fit)
     
171. validationplot(pls.fit,val.type="MSEP")
     
172. pls.pred=predict(pls.fit,x[test,],ncomp=2)
     
173. mean((pls.pred-y.test)^2)
     
174. pls.fit=plsr(Salary~., data=Hitters,scale=TRUE,ncomp=2)
     
175. summary(pls.fit)

复制代码

1. # Chapter 7 Lab: Non-linear Modeling
   
2. 
   
3. library(ISLR)
   
4. attach(Wage)
   
5. 
   
6. # Polynomial Regression and Step Functions
   
7. 
   
8. fit=lm(wage~poly(age,4),data=Wage)
   
9. coef(summary(fit))
   
10. fit2=lm(wage~poly(age,4,raw=T),data=Wage)
    
11. coef(summary(fit2))
    
12. fit2a=lm(wage~age+I(age^2)+I(age^3)+I(age^4),data=Wage)
    
13. coef(fit2a)
    
14. fit2b=lm(wage~cbind(age,age^2,age^3,age^4),data=Wage)
    
15. agelims=range(age)
    
16. age.grid=seq(from=agelims[1],to=agelims[2])
    
17. preds=predict(fit,newdata=list(age=age.grid),se=TRUE)
    
18. se.bands=cbind(preds$fit+2*preds$se.fit,preds$fit-2*preds$se.fit)
    
19. par(mfrow=c(1,2),mar=c(4.5,4.5,1,1),oma=c(0,0,4,0))
    
20. plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
    
21. title("Degree-4 Polynomial",outer=T)
    
22. lines(age.grid,preds$fit,lwd=2,col="blue")
    
23. matlines(age.grid,se.bands,lwd=1,col="blue",lty=3)
    
24. preds2=predict(fit2,newdata=list(age=age.grid),se=TRUE)
    
25. max(abs(preds$fit-preds2$fit))
    
26. fit.1=lm(wage~age,data=Wage)
    
27. fit.2=lm(wage~poly(age,2),data=Wage)
    
28. fit.3=lm(wage~poly(age,3),data=Wage)
    
29. fit.4=lm(wage~poly(age,4),data=Wage)
    
30. fit.5=lm(wage~poly(age,5),data=Wage)
    
31. anova(fit.1,fit.2,fit.3,fit.4,fit.5)
    
32. coef(summary(fit.5))
    
33. (-11.983)^2
    
34. fit.1=lm(wage~education+age,data=Wage)
    
35. fit.2=lm(wage~education+poly(age,2),data=Wage)
    
36. fit.3=lm(wage~education+poly(age,3),data=Wage)
    
37. anova(fit.1,fit.2,fit.3)
    
38. fit=glm(I(wage>250)~poly(age,4),data=Wage,family=binomial)
    
39. preds=predict(fit,newdata=list(age=age.grid),se=T)
    
40. pfit=exp(preds$fit)/(1+exp(preds$fit))
    
41. se.bands.logit = cbind(preds$fit+2*preds$se.fit, preds$fit-2*preds$se.fit)
    
42. se.bands = exp(se.bands.logit)/(1+exp(se.bands.logit))
    
43. preds=predict(fit,newdata=list(age=age.grid),type="response",se=T)
    
44. plot(age,I(wage>250),xlim=agelims,type="n",ylim=c(0,.2))
    
45. points(jitter(age), I((wage>250)/5),cex=.5,pch="|",col="darkgrey")
    
46. lines(age.grid,pfit,lwd=2, col="blue")
    
47. matlines(age.grid,se.bands,lwd=1,col="blue",lty=3)
    
48. table(cut(age,4))
    
49. fit=lm(wage~cut(age,4),data=Wage)
    
50. coef(summary(fit))
    
51. 
    
52. # Splines
    
53. 
    
54. library(splines)
    
55. fit=lm(wage~bs(age,knots=c(25,40,60)),data=Wage)
    
56. pred=predict(fit,newdata=list(age=age.grid),se=T)
    
57. plot(age,wage,col="gray")
    
58. lines(age.grid,pred$fit,lwd=2)
    
59. lines(age.grid,pred$fit+2*pred$se,lty="dashed")
    
60. lines(age.grid,pred$fit-2*pred$se,lty="dashed")
    
61. dim(bs(age,knots=c(25,40,60)))
    
62. dim(bs(age,df=6))
    
63. attr(bs(age,df=6),"knots")
    
64. fit2=lm(wage~ns(age,df=4),data=Wage)
    
65. pred2=predict(fit2,newdata=list(age=age.grid),se=T)
    
66. lines(age.grid, pred2$fit,col="red",lwd=2)
    
67. plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
    
68. title("Smoothing Spline")
    
69. fit=smooth.spline(age,wage,df=16)
    
70. fit2=smooth.spline(age,wage,cv=TRUE)
    
71. fit2$df
    
72. lines(fit,col="red",lwd=2)
    
73. lines(fit2,col="blue",lwd=2)
    
74. legend("topright",legend=c("16 DF","6.8 DF"),col=c("red","blue"),lty=1,lwd=2,cex=.8)
    
75. plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
    
76. title("Local Regression")
    
77. fit=loess(wage~age,span=.2,data=Wage)
    
78. fit2=loess(wage~age,span=.5,data=Wage)
    
79. lines(age.grid,predict(fit,data.frame(age=age.grid)),col="red",lwd=2)
    
80. lines(age.grid,predict(fit2,data.frame(age=age.grid)),col="blue",lwd=2)
    
81. legend("topright",legend=c("Span=0.2","Span=0.5"),col=c("red","blue"),lty=1,lwd=2,cex=.8)
    
82. 
    
83. # GAMs
    
84. 
    
85. gam1=lm(wage~ns(year,4)+ns(age,5)+education,data=Wage)
    
86. library(gam)
    
87. gam.m3=gam(wage~s(year,4)+s(age,5)+education,data=Wage)
    
88. par(mfrow=c(1,3))
    
89. plot(gam.m3, se=TRUE,col="blue")
    
90. plot.gam(gam1, se=TRUE, col="red")
    
91. gam.m1=gam(wage~s(age,5)+education,data=Wage)
    
92. gam.m2=gam(wage~year+s(age,5)+education,data=Wage)
    
93. anova(gam.m1,gam.m2,gam.m3,test="F")
    
94. summary(gam.m3)
    
95. preds=predict(gam.m2,newdata=Wage)
    
96. gam.lo=gam(wage~s(year,df=4)+lo(age,span=0.7)+education,data=Wage)
    
97. plot.gam(gam.lo, se=TRUE, col="green")
    
98. gam.lo.i=gam(wage~lo(year,age,span=0.5)+education,data=Wage)
    
99. library(akima)
    
100. plot(gam.lo.i)
     
101. gam.lr=gam(I(wage>250)~year+s(age,df=5)+education,family=binomial,data=Wage)
     
102. par(mfrow=c(1,3))
     
103. plot(gam.lr,se=T,col="green")
     
104. table(education,I(wage>250))
     
105. gam.lr.s=gam(I(wage>250)~year+s(age,df=5)+education,family=binomial,data=Wage,subset=(education!="1. < HS Grad"))
     
106. plot(gam.lr.s,se=T,col="green")

复制代码

1. # Chapter 8 Lab: Decision Trees
   
2. 
   
3. # Fitting Classification Trees
   
4. 
   
5. library(tree)
   
6. library(ISLR)
   
7. attach(Carseats)
   
8. High=ifelse(Sales<=8,"No","Yes")
   
9. Carseats=data.frame(Carseats,High)
   
10. tree.carseats=tree(High~.-Sales,Carseats)
    
11. summary(tree.carseats)
    
12. plot(tree.carseats)
    
13. text(tree.carseats,pretty=0)
    
14. tree.carseats
    
15. set.seed(2)
    
16. train=sample(1:nrow(Carseats), 200)
    
17. Carseats.test=Carseats[-train,]
    
18. High.test=High[-train]
    
19. tree.carseats=tree(High~.-Sales,Carseats,subset=train)
    
20. tree.pred=predict(tree.carseats,Carseats.test,type="class")
    
21. table(tree.pred,High.test)
    
22. (86+57)/200
    
23. set.seed(3)
    
24. cv.carseats=cv.tree(tree.carseats,FUN=prune.misclass)
    
25. names(cv.carseats)
    
26. cv.carseats
    
27. par(mfrow=c(1,2))
    
28. plot(cv.carseats$size,cv.carseats$dev,type="b")
    
29. plot(cv.carseats$k,cv.carseats$dev,type="b")
    
30. prune.carseats=prune.misclass(tree.carseats,best=9)
    
31. plot(prune.carseats)
    
32. text(prune.carseats,pretty=0)
    
33. tree.pred=predict(prune.carseats,Carseats.test,type="class")
    
34. table(tree.pred,High.test)
    
35. (94+60)/200
    
36. prune.carseats=prune.misclass(tree.carseats,best=15)
    
37. plot(prune.carseats)
    
38. text(prune.carseats,pretty=0)
    
39. tree.pred=predict(prune.carseats,Carseats.test,type="class")
    
40. table(tree.pred,High.test)
    
41. (86+62)/200
    
42. 
    
43. # Fitting Regression Trees
    
44. 
    
45. library(MASS)
    
46. set.seed(1)
    
47. train = sample(1:nrow(Boston), nrow(Boston)/2)
    
48. tree.boston=tree(medv~.,Boston,subset=train)
    
49. summary(tree.boston)
    
50. plot(tree.boston)
    
51. text(tree.boston,pretty=0)
    
52. cv.boston=cv.tree(tree.boston)
    
53. plot(cv.boston$size,cv.boston$dev,type='b')
    
54. prune.boston=prune.tree(tree.boston,best=5)
    
55. plot(prune.boston)
    
56. text(prune.boston,pretty=0)
    
57. yhat=predict(tree.boston,newdata=Boston[-train,])
    
58. boston.test=Boston[-train,"medv"]
    
59. plot(yhat,boston.test)
    
60. abline(0,1)
    
61. mean((yhat-boston.test)^2)
    
62. 
    
63. # Bagging and Random Forests
    
64. 
    
65. library(randomForest)
    
66. set.seed(1)
    
67. bag.boston=randomForest(medv~.,data=Boston,subset=train,mtry=13,importance=TRUE)
    
68. bag.boston
    
69. yhat.bag = predict(bag.boston,newdata=Boston[-train,])
    
70. plot(yhat.bag, boston.test)
    
71. abline(0,1)
    
72. mean((yhat.bag-boston.test)^2)
    
73. bag.boston=randomForest(medv~.,data=Boston,subset=train,mtry=13,ntree=25)
    
74. yhat.bag = predict(bag.boston,newdata=Boston[-train,])
    
75. mean((yhat.bag-boston.test)^2)
    
76. set.seed(1)
    
77. rf.boston=randomForest(medv~.,data=Boston,subset=train,mtry=6,importance=TRUE)
    
78. yhat.rf = predict(rf.boston,newdata=Boston[-train,])
    
79. mean((yhat.rf-boston.test)^2)
    
80. importance(rf.boston)
    
81. varImpPlot(rf.boston)
    
82. 
    
83. # Boosting
    
84. 
    
85. library(gbm)
    
86. set.seed(1)
    
87. boost.boston=gbm(medv~.,data=Boston[train,],distribution="gaussian",n.trees=5000,interaction.depth=4)
    
88. summary(boost.boston)
    
89. par(mfrow=c(1,2))
    
90. plot(boost.boston,i="rm")
    
91. plot(boost.boston,i="lstat")
    
92. yhat.boost=predict(boost.boston,newdata=Boston[-train,],n.trees=5000)
    
93. mean((yhat.boost-boston.test)^2)
    
94. boost.boston=gbm(medv~.,data=Boston[train,],distribution="gaussian",n.trees=5000,interaction.depth=4,shrinkage=0.2,verbose=F)
    
95. yhat.boost=predict(boost.boston,newdata=Boston[-train,],n.trees=5000)
    
96. mean((yhat.boost-boston.test)^2)

复制代码

1. # Chapter 9 Lab: Support Vector Machines
   
2. 
   
3. # Support Vector Classifier
   
4. 
   
5. set.seed(1)
   
6. x=matrix(rnorm(20*2), ncol=2)
   
7. y=c(rep(-1,10), rep(1,10))
   
8. x[y==1,]=x[y==1,] + 1
   
9. plot(x, col=(3-y))
   
10. dat=data.frame(x=x, y=as.factor(y))
    
11. library(e1071)
    
12. svmfit=svm(y~., data=dat, kernel="linear", cost=10,scale=FALSE)
    
13. plot(svmfit, dat)
    
14. svmfit$index
    
15. summary(svmfit)
    
16. svmfit=svm(y~., data=dat, kernel="linear", cost=0.1,scale=FALSE)
    
17. plot(svmfit, dat)
    
18. svmfit$index
    
19. set.seed(1)
    
20. tune.out=tune(svm,y~.,data=dat,kernel="linear",ranges=list(cost=c(0.001, 0.01, 0.1, 1,5,10,100)))
    
21. summary(tune.out)
    
22. bestmod=tune.out$best.model
    
23. summary(bestmod)
    
24. xtest=matrix(rnorm(20*2), ncol=2)
    
25. ytest=sample(c(-1,1), 20, rep=TRUE)
    
26. xtest[ytest==1,]=xtest[ytest==1,] + 1
    
27. testdat=data.frame(x=xtest, y=as.factor(ytest))
    
28. ypred=predict(bestmod,testdat)
    
29. table(predict=ypred, truth=testdat$y)
    
30. svmfit=svm(y~., data=dat, kernel="linear", cost=.01,scale=FALSE)
    
31. ypred=predict(svmfit,testdat)
    
32. table(predict=ypred, truth=testdat$y)
    
33. x[y==1,]=x[y==1,]+0.5
    
34. plot(x, col=(y+5)/2, pch=19)
    
35. dat=data.frame(x=x,y=as.factor(y))
    
36. svmfit=svm(y~., data=dat, kernel="linear", cost=1e5)
    
37. summary(svmfit)
    
38. plot(svmfit, dat)
    
39. svmfit=svm(y~., data=dat, kernel="linear", cost=1)
    
40. summary(svmfit)
    
41. plot(svmfit,dat)
    
42. 
    
43. # Support Vector Machine
    
44. 
    
45. set.seed(1)
    
46. x=matrix(rnorm(200*2), ncol=2)
    
47. x[1:100,]=x[1:100,]+2
    
48. x[101:150,]=x[101:150,]-2
    
49. y=c(rep(1,150),rep(2,50))
    
50. dat=data.frame(x=x,y=as.factor(y))
    
51. plot(x, col=y)
    
52. train=sample(200,100)
    
53. svmfit=svm(y~., data=dat[train,], kernel="radial",  gamma=1, cost=1)
    
54. plot(svmfit, dat[train,])
    
55. summary(svmfit)
    
56. svmfit=svm(y~., data=dat[train,], kernel="radial",gamma=1,cost=1e5)
    
57. plot(svmfit,dat[train,])
    
58. set.seed(1)
    
59. tune.out=tune(svm, y~., data=dat[train,], kernel="radial", ranges=list(cost=c(0.1,1,10,100,1000),gamma=c(0.5,1,2,3,4)))
    
60. summary(tune.out)
    
61. table(true=dat[-train,"y"], pred=predict(tune.out$best.model,newx=dat[-train,]))
    
62. 
    
63. # ROC Curves
    
64. 
    
65. library(ROCR)
    
66. rocplot=function(pred, truth, ...){
    
67.    predob = prediction(pred, truth)
    
68.    perf = performance(predob, "tpr", "fpr")
    
69.    plot(perf,...)}
    
70. svmfit.opt=svm(y~., data=dat[train,], kernel="radial",gamma=2, cost=1,decision.values=T)
    
71. fitted=attributes(predict(svmfit.opt,dat[train,],decision.values=TRUE))$decision.values
    
72. par(mfrow=c(1,2))
    
73. rocplot(fitted,dat[train,"y"],main="Training Data")
    
74. svmfit.flex=svm(y~., data=dat[train,], kernel="radial",gamma=50, cost=1, decision.values=T)
    
75. fitted=attributes(predict(svmfit.flex,dat[train,],decision.values=T))$decision.values
    
76. rocplot(fitted,dat[train,"y"],add=T,col="red")
    
77. fitted=attributes(predict(svmfit.opt,dat[-train,],decision.values=T))$decision.values
    
78. rocplot(fitted,dat[-train,"y"],main="Test Data")
    
79. fitted=attributes(predict(svmfit.flex,dat[-train,],decision.values=T))$decision.values
    
80. rocplot(fitted,dat[-train,"y"],add=T,col="red")
    
81. 
    
82. # SVM with Multiple Classes
    
83. 
    
84. set.seed(1)
    
85. x=rbind(x, matrix(rnorm(50*2), ncol=2))
    
86. y=c(y, rep(0,50))
    
87. x[y==0,2]=x[y==0,2]+2
    
88. dat=data.frame(x=x, y=as.factor(y))
    
89. par(mfrow=c(1,1))
    
90. plot(x,col=(y+1))
    
91. svmfit=svm(y~., data=dat, kernel="radial", cost=10, gamma=1)
    
92. plot(svmfit, dat)
    
93. 
    
94. # Application to Gene Expression Data
    
95. 
    
96. library(ISLR)
    
97. names(Khan)
    
98. dim(Khan$xtrain)
    
99. dim(Khan$xtest)
    
100. length(Khan$ytrain)
     
101. length(Khan$ytest)
     
102. table(Khan$ytrain)
     
103. table(Khan$ytest)
     
104. dat=data.frame(x=Khan$xtrain, y=as.factor(Khan$ytrain))
     
105. out=svm(y~., data=dat, kernel="linear",cost=10)
     
106. summary(out)
     
107. table(out$fitted, dat$y)
     
108. dat.te=data.frame(x=Khan$xtest, y=as.factor(Khan$ytest))
     
109. pred.te=predict(out, newdata=dat.te)
     
110. table(pred.te, dat.te$y)

复制代码

1. # Chapter 10 Lab 1: Principal Components Analysis
   
2. 
   
3. states=row.names(USArrests)
   
4. states
   
5. names(USArrests)
   
6. apply(USArrests, 2, mean)
   
7. apply(USArrests, 2, var)
   
8. pr.out=prcomp(USArrests, scale=TRUE)
   
9. names(pr.out)
   
10. pr.out$center
    
11. pr.out$scale
    
12. pr.out$rotation
    
13. dim(pr.out$x)
    
14. biplot(pr.out, scale=0)
    
15. pr.out$rotation=-pr.out$rotation
    
16. pr.out$x=-pr.out$x
    
17. biplot(pr.out, scale=0)
    
18. pr.out$sdev
    
19. pr.var=pr.out$sdev^2
    
20. pr.var
    
21. pve=pr.var/sum(pr.var)
    
22. pve
    
23. plot(pve, xlab="Principal Component", ylab="Proportion of Variance Explained", ylim=c(0,1),type='b')
    
24. plot(cumsum(pve), xlab="Principal Component", ylab="Cumulative Proportion of Variance Explained", ylim=c(0,1),type='b')
    
25. a=c(1,2,8,-3)
    
26. cumsum(a)
    
27. 
    
28. 
    
29. # Chapter 10 Lab 2: Clustering
    
30. 
    
31. # K-Means Clustering
    
32. 
    
33. set.seed(2)
    
34. x=matrix(rnorm(50*2), ncol=2)
    
35. x[1:25,1]=x[1:25,1]+3
    
36. x[1:25,2]=x[1:25,2]-4
    
37. km.out=kmeans(x,2,nstart=20)
    
38. km.out$cluster
    
39. plot(x, col=(km.out$cluster+1), main="K-Means Clustering Results with K=2", xlab="", ylab="", pch=20, cex=2)
    
40. set.seed(4)
    
41. km.out=kmeans(x,3,nstart=20)
    
42. km.out
    
43. plot(x, col=(km.out$cluster+1), main="K-Means Clustering Results with K=3", xlab="", ylab="", pch=20, cex=2)
    
44. set.seed(3)
    
45. km.out=kmeans(x,3,nstart=1)
    
46. km.out$tot.withinss
    
47. km.out=kmeans(x,3,nstart=20)
    
48. km.out$tot.withinss
    
49. 
    
50. # Hierarchical Clustering
    
51. 
    
52. hc.complete=hclust(dist(x), method="complete")
    
53. hc.average=hclust(dist(x), method="average")
    
54. hc.single=hclust(dist(x), method="single")
    
55. par(mfrow=c(1,3))
    
56. plot(hc.complete,main="Complete Linkage", xlab="", sub="", cex=.9)
    
57. plot(hc.average, main="Average Linkage", xlab="", sub="", cex=.9)
    
58. plot(hc.single, main="Single Linkage", xlab="", sub="", cex=.9)
    
59. cutree(hc.complete, 2)
    
60. cutree(hc.average, 2)
    
61. cutree(hc.single, 2)
    
62. cutree(hc.single, 4)
    
63. xsc=scale(x)
    
64. plot(hclust(dist(xsc), method="complete"), main="Hierarchical Clustering with Scaled Features")
    
65. x=matrix(rnorm(30*3), ncol=3)
    
66. dd=as.dist(1-cor(t(x)))
    
67. plot(hclust(dd, method="complete"), main="Complete Linkage with Correlation-Based Distance", xlab="", sub="")
    
68. 
    
69. 
    
70. # Chapter 10 Lab 3: NCI60 Data Example
    
71. 
    
72. # The NCI60 data
    
73. 
    
74. library(ISLR)
    
75. nci.labs=NCI60$labs
    
76. nci.data=NCI60$data
    
77. dim(nci.data)
    
78. nci.labs[1:4]
    
79. table(nci.labs)
    
80. 
    
81. # PCA on the NCI60 Data
    
82. 
    
83. pr.out=prcomp(nci.data, scale=TRUE)
    
84. Cols=function(vec){
    
85.     cols=rainbow(length(unique(vec)))
    
86.     return(cols[as.numeric(as.factor(vec))])
    
87.   }
    
88. par(mfrow=c(1,2))
    
89. plot(pr.out$x[,1:2], col=Cols(nci.labs), pch=19,xlab="Z1",ylab="Z2")
    
90. plot(pr.out$x[,c(1,3)], col=Cols(nci.labs), pch=19,xlab="Z1",ylab="Z3")
    
91. summary(pr.out)
    
92. plot(pr.out)
    
93. pve=100*pr.out$sdev^2/sum(pr.out$sdev^2)
    
94. par(mfrow=c(1,2))
    
95. plot(pve,  type="o", ylab="PVE", xlab="Principal Component", col="blue")
    
96. plot(cumsum(pve), type="o", ylab="Cumulative PVE", xlab="Principal Component", col="brown3")
    
97. 
    
98. # Clustering the Observations of the NCI60 Data
    
99. 
    
100. sd.data=scale(nci.data)
     
101. par(mfrow=c(1,3))
     
102. data.dist=dist(sd.data)
     
103. plot(hclust(data.dist), labels=nci.labs, main="Complete Linkage", xlab="", sub="",ylab="")
     
104. plot(hclust(data.dist, method="average"), labels=nci.labs, main="Average Linkage", xlab="", sub="",ylab="")
     
105. plot(hclust(data.dist, method="single"), labels=nci.labs,  main="Single Linkage", xlab="", sub="",ylab="")
     
106. hc.out=hclust(dist(sd.data))
     
107. hc.clusters=cutree(hc.out,4)
     
108. table(hc.clusters,nci.labs)
     
109. par(mfrow=c(1,1))
     
110. plot(hc.out, labels=nci.labs)
     
111. abline(h=139, col="red")
     
112. hc.out
     
113. set.seed(2)
     
114. km.out=kmeans(sd.data, 4, nstart=20)
     
115. km.clusters=km.out$cluster
     
116. table(km.clusters,hc.clusters)
     
117. hc.out=hclust(dist(pr.out$x[,1:5]))
     
118. plot(hc.out, labels=nci.labs, main="Hier. Clust. on First Five Score Vectors")
     
119. table(cutree(hc.out,4), nci.labs)

复制代码