0% found this document useful (0 votes)
16 views18 pages

All Labs R File

All Labs R File

Uploaded by

ratan203
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as TXT, PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
16 views18 pages

All Labs R File

All Labs R File

Uploaded by

ratan203
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as TXT, PDF, TXT or read online on Scribd

# Chapter 2 Lab: Introduction to R

# Basic Commands
x <- c(1,3,2,5)
x
x = c(1,6,2)
x
y = c(1,4,3)
length(x)
length(y)
x+y
ls()
rm(x,y)
ls()
rm(list=ls())
?matrix
x=matrix(data=c(1,2,3,4), nrow=2, ncol=2)
x
x=matrix(c(1,2,3,4),2,2)
matrix(c(1,2,3,4),2,2,byrow=TRUE)
sqrt(x)
x^2
x=rnorm(50)
y=x+rnorm(50,mean=50,sd=.1)
cor(x,y)
[Link](1303)
rnorm(50)
[Link](3)
y=rnorm(100)
mean(y)
var(y)
sqrt(var(y))
sd(y)
# Graphics
x=rnorm(100)
y=rnorm(100)
plot(x,y)
plot(x,y,xlab="this is the x-axis",ylab="this is the y-axis",main="Plot of X vs
Y")
pdf("[Link]")
plot(x,y,col="green")
[Link]()
x=seq(1,10)
x
x=1:10
x
x=seq(-pi,pi,length=50)
y=x
f=outer(x,y,function(x,y)cos(y)/(1+x^2))
contour(x,y,f)
contour(x,y,f,nlevels=45,add=T)
fa=(f-t(f))/2
contour(x,y,fa,nlevels=15)
image(x,y,fa)
persp(x,y,fa)
persp(x,y,fa,theta=30)
persp(x,y,fa,theta=30,phi=20)

persp(x,y,fa,theta=30,phi=70)
persp(x,y,fa,theta=30,phi=40)
# Indexing Data
A=matrix(1:16,4,4)
A
A[2,3]
A[c(1,3),c(2,4)]
A[1:3,2:4]
A[1:2,]
A[,1:2]
A[1,]
A[-c(1,3),]
A[-c(1,3),-c(1,3,4)]
dim(A)
# Loading Data
Auto=[Link]("[Link]")
fix(Auto)
Auto=[Link]("[Link]",header=T,[Link]="?")
fix(Auto)
Auto=[Link]("[Link]",header=T,[Link]="?")
fix(Auto)
dim(Auto)
Auto[1:4,]
Auto=[Link](Auto)
dim(Auto)
names(Auto)
# Additional Graphical and Numerical Summaries
plot(cylinders, mpg)
plot(Auto$cylinders, Auto$mpg)
attach(Auto)
plot(cylinders, mpg)
cylinders=[Link](cylinders)
plot(cylinders, mpg)
plot(cylinders, mpg, col="red")
plot(cylinders, mpg, col="red", varwidth=T)
plot(cylinders, mpg, col="red", varwidth=T,horizontal=T)
plot(cylinders, mpg, col="red", varwidth=T, xlab="cylinders", ylab="MPG")
hist(mpg)
hist(mpg,col=2)
hist(mpg,col=2,breaks=15)
pairs(Auto)
pairs(~ mpg + displacement + horsepower + weight + acceleration, Auto)
plot(horsepower,mpg)
identify(horsepower,mpg,name)
summary(Auto)
summary(mpg)

# Chapter 3 Lab: Linear Regression


library(MASS)
library(ISLR)

# Simple Linear Regression


fix(Boston)
names(Boston)
[Link]=lm(medv~lstat)
[Link]=lm(medv~lstat,data=Boston)
attach(Boston)
[Link]=lm(medv~lstat)
[Link]
summary([Link])
names([Link])
coef([Link])
confint([Link])
predict([Link],[Link](lstat=(c(5,10,15))), interval="confidence")
predict([Link],[Link](lstat=(c(5,10,15))), interval="prediction")
plot(lstat,medv)
abline([Link])
abline([Link],lwd=3)
abline([Link],lwd=3,col="red")
plot(lstat,medv,col="red")
plot(lstat,medv,pch=20)
plot(lstat,medv,pch="+")
plot(1:20,1:20,pch=1:20)
par(mfrow=c(2,2))
plot([Link])
plot(predict([Link]), residuals([Link]))
plot(predict([Link]), rstudent([Link]))
plot(hatvalues([Link]))
[Link](hatvalues([Link]))
# Multiple Linear Regression
[Link]=lm(medv~lstat+age,data=Boston)
summary([Link])
[Link]=lm(medv~.,data=Boston)
summary([Link])
library(car)
vif([Link])
lm.fit1=lm(medv~.-age,data=Boston)
summary(lm.fit1)
lm.fit1=update([Link], ~.-age)
# Interaction Terms
summary(lm(medv~lstat*age,data=Boston))
# Non-linear Transformations of the Predictors
lm.fit2=lm(medv~lstat+I(lstat^2))
summary(lm.fit2)
[Link]=lm(medv~lstat)
anova([Link],lm.fit2)
par(mfrow=c(2,2))
plot(lm.fit2)
lm.fit5=lm(medv~poly(lstat,5))
summary(lm.fit5)
summary(lm(medv~log(rm),data=Boston))
# Qualitative Predictors

fix(Carseats)
names(Carseats)
[Link]=lm(Sales~.+Income:Advertising+Price:Age,data=Carseats)
summary([Link])
attach(Carseats)
contrasts(ShelveLoc)
# Writing Functions
LoadLibraries
LoadLibraries()
LoadLibraries=function(){
library(ISLR)
library(MASS)
print("The libraries have been loaded.")
}
LoadLibraries
LoadLibraries()
# Chapter 4 Lab: Logistic Regression, LDA, QDA, and KNN
# The Stock Market Data
library(ISLR)
names(Smarket)
dim(Smarket)
summary(Smarket)
pairs(Smarket)
cor(Smarket)
cor(Smarket[,-9])
attach(Smarket)
plot(Volume)
# Logistic Regression
[Link]=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomi
al)
summary([Link])
coef([Link])
summary([Link])$coef
summary([Link])$coef[,4]
[Link]=predict([Link],type="response")
[Link][1:10]
contrasts(Direction)
[Link]=rep("Down",1250)
[Link][[Link]>.5]="Up"
table([Link],Direction)
(507+145)/1250
mean([Link]==Direction)
train=(Year<2005)
Smarket.2005=Smarket[!train,]
dim(Smarket.2005)
Direction.2005=Direction[!train]
[Link]=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomi
al,subset=train)
[Link]=predict([Link],Smarket.2005,type="response")
[Link]=rep("Down",252)
[Link][[Link]>.5]="Up"
table([Link],Direction.2005)

mean([Link]==Direction.2005)
mean([Link]!=Direction.2005)
[Link]=glm(Direction~Lag1+Lag2,data=Smarket,family=binomial,subset=train)
[Link]=predict([Link],Smarket.2005,type="response")
[Link]=rep("Down",252)
[Link][[Link]>.5]="Up"
table([Link],Direction.2005)
mean([Link]==Direction.2005)
106/(106+76)
predict([Link],newdata=[Link](Lag1=c(1.2,1.5),Lag2=c(1.1,-0.8)),type="respo
nse")
# Linear Discriminant Analysis
library(MASS)
[Link]=lda(Direction~Lag1+Lag2,data=Smarket,subset=train)
[Link]
plot([Link])
[Link]=predict([Link], Smarket.2005)
names([Link])
[Link]=[Link]$class
table([Link],Direction.2005)
mean([Link]==Direction.2005)
sum([Link]$posterior[,1]>=.5)
sum([Link]$posterior[,1]<.5)
[Link]$posterior[1:20,1]
[Link][1:20]
sum([Link]$posterior[,1]>.9)
# Quadratic Discriminant Analysis
[Link]=qda(Direction~Lag1+Lag2,data=Smarket,subset=train)
[Link]
[Link]=predict([Link],Smarket.2005)$class
table([Link],Direction.2005)
mean([Link]==Direction.2005)
# K-Nearest Neighbors
library(class)
train.X=cbind(Lag1,Lag2)[train,]
test.X=cbind(Lag1,Lag2)[!train,]
[Link]=Direction[train]
[Link](1)
[Link]=knn(train.X,test.X,[Link],k=1)
table([Link],Direction.2005)
(83+43)/252
[Link]=knn(train.X,test.X,[Link],k=3)
table([Link],Direction.2005)
mean([Link]==Direction.2005)
# An Application to Caravan Insurance Data
dim(Caravan)
attach(Caravan)
summary(Purchase)
348/5822
standardized.X=scale(Caravan[,-86])
var(Caravan[,1])
var(Caravan[,2])

var(standardized.X[,1])
var(standardized.X[,2])
test=1:1000
train.X=standardized.X[-test,]
test.X=standardized.X[test,]
train.Y=Purchase[-test]
test.Y=Purchase[test]
[Link](1)
[Link]=knn(train.X,test.X,train.Y,k=1)
mean(test.Y!=[Link])
mean(test.Y!="No")
table([Link],test.Y)
9/(68+9)
[Link]=knn(train.X,test.X,train.Y,k=3)
table([Link],test.Y)
5/26
[Link]=knn(train.X,test.X,train.Y,k=5)
table([Link],test.Y)
4/15
[Link]=glm(Purchase~.,data=Caravan,family=binomial,subset=-test)
[Link]=predict([Link],Caravan[test,],type="response")
[Link]=rep("No",1000)
[Link][[Link]>.5]="Yes"
table([Link],test.Y)
[Link]=rep("No",1000)
[Link][[Link]>.25]="Yes"
table([Link],test.Y)
11/(22+11)

# Chaper 5 Lab: Cross-Validation and the Bootstrap


# The Validation Set Approach
library(ISLR)
[Link](1)
train=sample(392,196)
[Link]=lm(mpg~horsepower,data=Auto,subset=train)
attach(Auto)
mean((mpg-predict([Link],Auto))[-train]^2)
lm.fit2=lm(mpg~poly(horsepower,2),data=Auto,subset=train)
mean((mpg-predict(lm.fit2,Auto))[-train]^2)
lm.fit3=lm(mpg~poly(horsepower,3),data=Auto,subset=train)
mean((mpg-predict(lm.fit3,Auto))[-train]^2)
[Link](2)
train=sample(392,196)
[Link]=lm(mpg~horsepower,subset=train)
mean((mpg-predict([Link],Auto))[-train]^2)
lm.fit2=lm(mpg~poly(horsepower,2),data=Auto,subset=train)
mean((mpg-predict(lm.fit2,Auto))[-train]^2)
lm.fit3=lm(mpg~poly(horsepower,3),data=Auto,subset=train)
mean((mpg-predict(lm.fit3,Auto))[-train]^2)
# Leave-One-Out Cross-Validation
[Link]=glm(mpg~horsepower,data=Auto)
coef([Link])
[Link]=lm(mpg~horsepower,data=Auto)
coef([Link])

library(boot)
[Link]=glm(mpg~horsepower,data=Auto)
[Link]=[Link](Auto,[Link])
[Link]$delta
[Link]=rep(0,5)
for (i in 1:5){
[Link]=glm(mpg~poly(horsepower,i),data=Auto)
[Link][i]=[Link](Auto,[Link])$delta[1]
}
[Link]
# k-Fold Cross-Validation
[Link](17)
[Link].10=rep(0,10)
for (i in 1:10){
[Link]=glm(mpg~poly(horsepower,i),data=Auto)
[Link].10[i]=[Link](Auto,[Link],K=10)$delta[1]
}
[Link].10
# The Bootstrap
[Link]=function(data,index){
X=data$X[index]
Y=data$Y[index]
return((var(Y)-cov(X,Y))/(var(X)+var(Y)-2*cov(X,Y)))
}
[Link](Portfolio,1:100)
[Link](1)
[Link](Portfolio,sample(100,100,replace=T))
boot(Portfolio,[Link],R=1000)
# Estimating the Accuracy of a Linear Regression Model
[Link]=function(data,index)
return(coef(lm(mpg~horsepower,data=data,subset=index)))
[Link](Auto,1:392)
[Link](1)
[Link](Auto,sample(392,392,replace=T))
[Link](Auto,sample(392,392,replace=T))
boot(Auto,[Link],1000)
summary(lm(mpg~horsepower,data=Auto))$coef
[Link]=function(data,index)
coefficients(lm(mpg~horsepower+I(horsepower^2),data=data,subset=index))
[Link](1)
boot(Auto,[Link],1000)
summary(lm(mpg~horsepower+I(horsepower^2),data=Auto))$coef

# Chapter 6 Lab 1: Subset Selection Methods


# Best Subset Selection
library(ISLR)
fix(Hitters)
names(Hitters)
dim(Hitters)
sum([Link](Hitters$Salary))

Hitters=[Link](Hitters)
dim(Hitters)
sum([Link](Hitters))
library(leaps)
[Link]=regsubsets(Salary~.,Hitters)
summary([Link])
[Link]=regsubsets(Salary~.,data=Hitters,nvmax=19)
[Link]=summary([Link])
names([Link])
[Link]$rsq
par(mfrow=c(2,2))
plot([Link]$rss,xlab="Number of Variables",ylab="RSS",type="l")
plot([Link]$adjr2,xlab="Number of Variables",ylab="Adjusted RSq",type="l")
[Link]([Link]$adjr2)
points(11,[Link]$adjr2[11], col="red",cex=2,pch=20)
plot([Link]$cp,xlab="Number of Variables",ylab="Cp",type='l')
[Link]([Link]$cp)
points(10,[Link]$cp[10],col="red",cex=2,pch=20)
[Link]([Link]$bic)
plot([Link]$bic,xlab="Number of Variables",ylab="BIC",type='l')
points(6,[Link]$bic[6],col="red",cex=2,pch=20)
plot([Link],scale="r2")
plot([Link],scale="adjr2")
plot([Link],scale="Cp")
plot([Link],scale="bic")
coef([Link],6)
# Forward and Backward Stepwise Selection
[Link]=regsubsets(Salary~.,data=Hitters,nvmax=19,method="forward")
summary([Link])
[Link]=regsubsets(Salary~.,data=Hitters,nvmax=19,method="backward")
summary([Link])
coef([Link],7)
coef([Link],7)
coef([Link],7)
# Choosing Among Models
[Link](1)
train=sample(c(TRUE,FALSE), nrow(Hitters),rep=TRUE)
test=(!train)
[Link]=regsubsets(Salary~.,data=Hitters[train,],nvmax=19)
[Link]=[Link](Salary~.,data=Hitters[test,])
[Link]=rep(NA,19)
for(i in 1:19){
coefi=coef([Link],id=i)
pred=[Link][,names(coefi)]%*%coefi
[Link][i]=mean((Hitters$Salary[test]-pred)^2)
}
[Link]
[Link]([Link])
coef([Link],10)
[Link]=function(object,newdata,id,...){
form=[Link](object$call[[2]])
mat=[Link](form,newdata)
coefi=coef(object,id=id)
xvars=names(coefi)
mat[,xvars]%*%coefi
}

[Link]=regsubsets(Salary~.,data=Hitters,nvmax=19)
coef([Link],10)
k=10
[Link](1)
folds=sample(1:k,nrow(Hitters),replace=TRUE)
[Link]=matrix(NA,k,19, dimnames=list(NULL, paste(1:19)))
for(j in 1:k){
[Link]=regsubsets(Salary~.,data=Hitters[folds!=j,],nvmax=19)
for(i in 1:19){
pred=predict([Link],Hitters[folds==j,],id=i)
[Link][j,i]=mean( (Hitters$Salary[folds==j]-pred)^2)
}
}
[Link]=apply([Link],2,mean)
[Link]
par(mfrow=c(1,1))
plot([Link],type='b')
[Link]=regsubsets(Salary~.,data=Hitters, nvmax=19)
coef([Link],11)
# Chapter 6 Lab 2: Ridge Regression and the Lasso
x=[Link](Salary~.,Hitters)[,-1]
y=Hitters$Salary
# Ridge Regression
library(glmnet)
grid=10^seq(10,-2,length=100)
[Link]=glmnet(x,y,alpha=0,lambda=grid)
dim(coef([Link]))
[Link]$lambda[50]
coef([Link])[,50]
sqrt(sum(coef([Link])[-1,50]^2))
[Link]$lambda[60]
coef([Link])[,60]
sqrt(sum(coef([Link])[-1,60]^2))
predict([Link],s=50,type="coefficients")[1:20,]
[Link](1)
train=sample(1:nrow(x), nrow(x)/2)
test=(-train)
[Link]=y[test]
[Link]=glmnet(x[train,],y[train],alpha=0,lambda=grid, thresh=1e-12)
[Link]=predict([Link],s=4,newx=x[test,])
mean(([Link])^2)
mean((mean(y[train])-[Link])^2)
[Link]=predict([Link],s=1e10,newx=x[test,])
mean(([Link])^2)
[Link]=predict([Link],s=0,newx=x[test,],exact=T)
mean(([Link])^2)
lm(y~x, subset=train)
predict([Link],s=0,exact=T,type="coefficients")[1:20,]
[Link](1)
[Link]=[Link](x[train,],y[train],alpha=0)
plot([Link])
bestlam=[Link]$[Link]
bestlam
[Link]=predict([Link],s=bestlam,newx=x[test,])
mean(([Link])^2)

out=glmnet(x,y,alpha=0)
predict(out,type="coefficients",s=bestlam)[1:20,]
# The Lasso
[Link]=glmnet(x[train,],y[train],alpha=1,lambda=grid)
plot([Link])
[Link](1)
[Link]=[Link](x[train,],y[train],alpha=1)
plot([Link])
bestlam=[Link]$[Link]
[Link]=predict([Link],s=bestlam,newx=x[test,])
mean(([Link])^2)
out=glmnet(x,y,alpha=1,lambda=grid)
[Link]=predict(out,type="coefficients",s=bestlam)[1:20,]
[Link]
[Link][[Link]!=0]
# Chapter 6 Lab 3: PCR and PLS Regression
# Principal Components Regression
library(pls)
[Link](2)
[Link]=pcr(Salary~., data=Hitters,scale=TRUE,validation="CV")
summary([Link])
validationplot([Link],[Link]="MSEP")
[Link](1)
[Link]=pcr(Salary~., data=Hitters,subset=train,scale=TRUE, validation="CV")
validationplot([Link],[Link]="MSEP")
[Link]=predict([Link],x[test,],ncomp=7)
mean(([Link])^2)
[Link]=pcr(y~x,scale=TRUE,ncomp=7)
summary([Link])
# Partial Least Squares
[Link](1)
[Link]=plsr(Salary~., data=Hitters,subset=train,scale=TRUE, validation="CV")
summary([Link])
validationplot([Link],[Link]="MSEP")
[Link]=predict([Link],x[test,],ncomp=2)
mean(([Link])^2)
[Link]=plsr(Salary~., data=Hitters,scale=TRUE,ncomp=2)
summary([Link])

# Chapter 7 Lab: Non-linear Modeling


library(ISLR)
attach(Wage)
# Polynomial Regression and Step Functions
fit=lm(wage~poly(age,4),data=Wage)
coef(summary(fit))
fit2=lm(wage~poly(age,4,raw=T),data=Wage)
coef(summary(fit2))

fit2a=lm(wage~age+I(age^2)+I(age^3)+I(age^4),data=Wage)
coef(fit2a)
fit2b=lm(wage~cbind(age,age^2,age^3,age^4),data=Wage)
agelims=range(age)
[Link]=seq(from=agelims[1],to=agelims[2])
preds=predict(fit,newdata=list(age=[Link]),se=TRUE)
[Link]=cbind(preds$fit+2*preds$[Link],preds$fit-2*preds$[Link])
par(mfrow=c(1,2),mar=c(4.5,4.5,1,1),oma=c(0,0,4,0))
plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
title("Degree-4 Polynomial",outer=T)
lines([Link],preds$fit,lwd=2,col="blue")
matlines([Link],[Link],lwd=1,col="blue",lty=3)
preds2=predict(fit2,newdata=list(age=[Link]),se=TRUE)
max(abs(preds$fit-preds2$fit))
fit.1=lm(wage~age,data=Wage)
fit.2=lm(wage~poly(age,2),data=Wage)
fit.3=lm(wage~poly(age,3),data=Wage)
fit.4=lm(wage~poly(age,4),data=Wage)
fit.5=lm(wage~poly(age,5),data=Wage)
anova(fit.1,fit.2,fit.3,fit.4,fit.5)
coef(summary(fit.5))
(-11.983)^2
fit.1=lm(wage~education+age,data=Wage)
fit.2=lm(wage~education+poly(age,2),data=Wage)
fit.3=lm(wage~education+poly(age,3),data=Wage)
anova(fit.1,fit.2,fit.3)
fit=glm(I(wage>250)~poly(age,4),data=Wage,family=binomial)
preds=predict(fit,newdata=list(age=[Link]),se=T)
pfit=exp(preds$fit)/(1+exp(preds$fit))
[Link] = cbind(preds$fit+2*preds$[Link], preds$fit-2*preds$[Link])
[Link] = exp([Link])/(1+exp([Link]))
preds=predict(fit,newdata=list(age=[Link]),type="response",se=T)
plot(age,I(wage>250),xlim=agelims,type="n",ylim=c(0,.2))
points(jitter(age), I((wage>250)/5),cex=.5,pch="|",col="darkgrey")
lines([Link],pfit,lwd=2, col="blue")
matlines([Link],[Link],lwd=1,col="blue",lty=3)
table(cut(age,4))
fit=lm(wage~cut(age,4),data=Wage)
coef(summary(fit))
# Splines
library(splines)
fit=lm(wage~bs(age,knots=c(25,40,60)),data=Wage)
pred=predict(fit,newdata=list(age=[Link]),se=T)
plot(age,wage,col="gray")
lines([Link],pred$fit,lwd=2)
lines([Link],pred$fit+2*pred$se,lty="dashed")
lines([Link],pred$fit-2*pred$se,lty="dashed")
dim(bs(age,knots=c(25,40,60)))
dim(bs(age,df=6))
attr(bs(age,df=6),"knots")
fit2=lm(wage~ns(age,df=4),data=Wage)
pred2=predict(fit2,newdata=list(age=[Link]),se=T)
lines([Link], pred2$fit,col="red",lwd=2)
plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
title("Smoothing Spline")
fit=[Link](age,wage,df=16)
fit2=[Link](age,wage,cv=TRUE)
fit2$df

lines(fit,col="red",lwd=2)
lines(fit2,col="blue",lwd=2)
legend("topright",legend=c("16 DF","6.8 DF"),col=c("red","blue"),lty=1,lwd=2,cex
=.8)
plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
title("Local Regression")
fit=loess(wage~age,span=.2,data=Wage)
fit2=loess(wage~age,span=.5,data=Wage)
lines([Link],predict(fit,[Link](age=[Link])),col="red",lwd=2)
lines([Link],predict(fit2,[Link](age=[Link])),col="blue",lwd=2)
legend("topright",legend=c("Span=0.2","Span=0.5"),col=c("red","blue"),lty=1,lwd=
2,cex=.8)
# GAMs
gam1=lm(wage~ns(year,4)+ns(age,5)+education,data=Wage)
library(gam)
gam.m3=gam(wage~s(year,4)+s(age,5)+education,data=Wage)
par(mfrow=c(1,3))
plot(gam.m3, se=TRUE,col="blue")
[Link](gam1, se=TRUE, col="red")
gam.m1=gam(wage~s(age,5)+education,data=Wage)
gam.m2=gam(wage~year+s(age,5)+education,data=Wage)
anova(gam.m1,gam.m2,gam.m3,test="F")
summary(gam.m3)
preds=predict(gam.m2,newdata=Wage)
[Link]=gam(wage~s(year,df=4)+lo(age,span=0.7)+education,data=Wage)
[Link]([Link], se=TRUE, col="green")
[Link].i=gam(wage~lo(year,age,span=0.5)+education,data=Wage)
library(akima)
plot([Link].i)
[Link]=gam(I(wage>250)~year+s(age,df=5)+education,family=binomial,data=Wage)
par(mfrow=c(1,3))
plot([Link],se=T,col="green")
table(education,I(wage>250))
[Link].s=gam(I(wage>250)~year+s(age,df=5)+education,family=binomial,data=Wage,su
bset=(education!="1. < HS Grad"))
plot([Link].s,se=T,col="green")

# Chapter 8 Lab: Decision Trees


# Fitting Classification Trees
library(tree)
library(ISLR)
attach(Carseats)
High=ifelse(Sales<=8,"No","Yes")
Carseats=[Link](Carseats,High)
[Link]=tree(High~.-Sales,Carseats)
summary([Link])
plot([Link])
text([Link],pretty=0)
[Link]
[Link](2)
train=sample(1:nrow(Carseats), 200)
[Link]=Carseats[-train,]
[Link]=High[-train]
[Link]=tree(High~.-Sales,Carseats,subset=train)

[Link]=predict([Link],[Link],type="class")
table([Link],[Link])
(86+57)/200
[Link](3)
[Link]=[Link]([Link],FUN=[Link])
names([Link])
[Link]
par(mfrow=c(1,2))
plot([Link]$size,[Link]$dev,type="b")
plot([Link]$k,[Link]$dev,type="b")
[Link]=[Link]([Link],best=9)
plot([Link])
text([Link],pretty=0)
[Link]=predict([Link],[Link],type="class")
table([Link],[Link])
(94+60)/200
[Link]=[Link]([Link],best=15)
plot([Link])
text([Link],pretty=0)
[Link]=predict([Link],[Link],type="class")
table([Link],[Link])
(86+62)/200
# Fitting Regression Trees
library(MASS)
[Link](1)
train = sample(1:nrow(Boston), nrow(Boston)/2)
[Link]=tree(medv~.,Boston,subset=train)
summary([Link])
plot([Link])
text([Link],pretty=0)
[Link]=[Link]([Link])
plot([Link]$size,[Link]$dev,type='b')
[Link]=[Link]([Link],best=5)
plot([Link])
text([Link],pretty=0)
yhat=predict([Link],newdata=Boston[-train,])
[Link]=Boston[-train,"medv"]
plot(yhat,[Link])
abline(0,1)
mean(([Link])^2)
# Bagging and Random Forests
library(randomForest)
[Link](1)
[Link]=randomForest(medv~.,data=Boston,subset=train,mtry=13,importance=TRUE)
[Link]
[Link] = predict([Link],newdata=Boston[-train,])
plot([Link], [Link])
abline(0,1)
mean(([Link])^2)
[Link]=randomForest(medv~.,data=Boston,subset=train,mtry=13,ntree=25)
[Link] = predict([Link],newdata=Boston[-train,])
mean(([Link])^2)
[Link](1)
[Link]=randomForest(medv~.,data=Boston,subset=train,mtry=6,importance=TRUE)
[Link] = predict([Link],newdata=Boston[-train,])
mean(([Link])^2)

importance([Link])
varImpPlot([Link])
# Boosting
library(gbm)
[Link](1)
[Link]=gbm(medv~.,data=Boston[train,],distribution="gaussian",[Link]=5000
,[Link]=4)
summary([Link])
par(mfrow=c(1,2))
plot([Link],i="rm")
plot([Link],i="lstat")
[Link]=predict([Link],newdata=Boston[-train,],[Link]=5000)
mean(([Link])^2)
[Link]=gbm(medv~.,data=Boston[train,],distribution="gaussian",[Link]=5000
,[Link]=4,shrinkage=0.2,verbose=F)
[Link]=predict([Link],newdata=Boston[-train,],[Link]=5000)
mean(([Link])^2)

# Chapter 9 Lab: Support Vector Machines


# Support Vector Classifier
[Link](1)
x=matrix(rnorm(20*2), ncol=2)
y=c(rep(-1,10), rep(1,10))
x[y==1,]=x[y==1,] + 1
plot(x, col=(3-y))
dat=[Link](x=x, y=[Link](y))
library(e1071)
svmfit=svm(y~., data=dat, kernel="linear", cost=10,scale=FALSE)
plot(svmfit, dat)
svmfit$index
summary(svmfit)
svmfit=svm(y~., data=dat, kernel="linear", cost=0.1,scale=FALSE)
plot(svmfit, dat)
svmfit$index
[Link](1)
[Link]=tune(svm,y~.,data=dat,kernel="linear",ranges=list(cost=c(0.001, 0.01, 0
.1, 1,5,10,100)))
summary([Link])
bestmod=[Link]$[Link]
summary(bestmod)
xtest=matrix(rnorm(20*2), ncol=2)
ytest=sample(c(-1,1), 20, rep=TRUE)
xtest[ytest==1,]=xtest[ytest==1,] + 1
testdat=[Link](x=xtest, y=[Link](ytest))
ypred=predict(bestmod,testdat)
table(predict=ypred, truth=testdat$y)
svmfit=svm(y~., data=dat, kernel="linear", cost=.01,scale=FALSE)
ypred=predict(svmfit,testdat)
table(predict=ypred, truth=testdat$y)
x[y==1,]=x[y==1,]+0.5
plot(x, col=(y+5)/2, pch=19)
dat=[Link](x=x,y=[Link](y))
svmfit=svm(y~., data=dat, kernel="linear", cost=1e5)
summary(svmfit)

plot(svmfit, dat)
svmfit=svm(y~., data=dat, kernel="linear", cost=1)
summary(svmfit)
plot(svmfit,dat)
# Support Vector Machine
[Link](1)
x=matrix(rnorm(200*2), ncol=2)
x[1:100,]=x[1:100,]+2
x[101:150,]=x[101:150,]-2
y=c(rep(1,150),rep(2,50))
dat=[Link](x=x,y=[Link](y))
plot(x, col=y)
train=sample(200,100)
svmfit=svm(y~., data=dat[train,], kernel="radial", gamma=1, cost=1)
plot(svmfit, dat[train,])
summary(svmfit)
svmfit=svm(y~., data=dat[train,], kernel="radial",gamma=1,cost=1e5)
plot(svmfit,dat[train,])
[Link](1)
[Link]=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)))
summary([Link])
table(true=dat[-train,"y"], pred=predict([Link]$[Link],newx=dat[-train,]))
# ROC Curves
library(ROCR)
rocplot=function(pred, truth, ...){
predob = prediction(pred, truth)
perf = performance(predob, "tpr", "fpr")
plot(perf,...)}
[Link]=svm(y~., data=dat[train,], kernel="radial",gamma=2, cost=1,decision.v
alues=T)
fitted=attributes(predict([Link],dat[train,],[Link]=TRUE))$decision
.values
par(mfrow=c(1,2))
rocplot(fitted,dat[train,"y"],main="Training Data")
[Link]=svm(y~., data=dat[train,], kernel="radial",gamma=50, cost=1, decisio
[Link]=T)
fitted=attributes(predict([Link],dat[train,],[Link]=T))$decision.v
alues
rocplot(fitted,dat[train,"y"],add=T,col="red")
fitted=attributes(predict([Link],dat[-train,],[Link]=T))$decision.v
alues
rocplot(fitted,dat[-train,"y"],main="Test Data")
fitted=attributes(predict([Link],dat[-train,],[Link]=T))$decision.
values
rocplot(fitted,dat[-train,"y"],add=T,col="red")
# SVM with Multiple Classes
[Link](1)
x=rbind(x, matrix(rnorm(50*2), ncol=2))
y=c(y, rep(0,50))
x[y==0,2]=x[y==0,2]+2
dat=[Link](x=x, y=[Link](y))
par(mfrow=c(1,1))
plot(x,col=(y+1))

svmfit=svm(y~., data=dat, kernel="radial", cost=10, gamma=1)


plot(svmfit, dat)
# Application to Gene Expression Data
library(ISLR)
names(Khan)
dim(Khan$xtrain)
dim(Khan$xtest)
length(Khan$ytrain)
length(Khan$ytest)
table(Khan$ytrain)
table(Khan$ytest)
dat=[Link](x=Khan$xtrain, y=[Link](Khan$ytrain))
out=svm(y~., data=dat, kernel="linear",cost=10)
summary(out)
table(out$fitted, dat$y)
[Link]=[Link](x=Khan$xtest, y=[Link](Khan$ytest))
[Link]=predict(out, newdata=[Link])
table([Link], [Link]$y)

# Chapter 10 Lab 1: Principal Components Analysis


states=[Link](USArrests)
states
names(USArrests)
apply(USArrests, 2, mean)
apply(USArrests, 2, var)
[Link]=prcomp(USArrests, scale=TRUE)
names([Link])
[Link]$center
[Link]$scale
[Link]$rotation
dim([Link]$x)
biplot([Link], scale=0)
[Link]$rotation=-[Link]$rotation
[Link]$x=-[Link]$x
biplot([Link], scale=0)
[Link]$sdev
[Link]=[Link]$sdev^2
[Link]
pve=[Link]/sum([Link])
pve
plot(pve, xlab="Principal Component", ylab="Proportion of Variance Explained", y
lim=c(0,1),type='b')
plot(cumsum(pve), xlab="Principal Component", ylab="Cumulative Proportion of Var
iance Explained", ylim=c(0,1),type='b')
a=c(1,2,8,-3)
cumsum(a)
# Chapter 10 Lab 2: Clustering
# K-Means Clustering
[Link](2)
x=matrix(rnorm(50*2), ncol=2)
x[1:25,1]=x[1:25,1]+3

x[1:25,2]=x[1:25,2]-4
[Link]=kmeans(x,2,nstart=20)
[Link]$cluster
plot(x, col=([Link]$cluster+1), main="K-Means Clustering Results with K=2", xlab
="", ylab="", pch=20, cex=2)
[Link](4)
[Link]=kmeans(x,3,nstart=20)
[Link]
plot(x, col=([Link]$cluster+1), main="K-Means Clustering Results with K=3", xlab
="", ylab="", pch=20, cex=2)
[Link](3)
[Link]=kmeans(x,3,nstart=1)
[Link]$[Link]
[Link]=kmeans(x,3,nstart=20)
[Link]$[Link]
# Hierarchical Clustering
[Link]=hclust(dist(x), method="complete")
[Link]=hclust(dist(x), method="average")
[Link]=hclust(dist(x), method="single")
par(mfrow=c(1,3))
plot([Link],main="Complete Linkage", xlab="", sub="", cex=.9)
plot([Link], main="Average Linkage", xlab="", sub="", cex=.9)
plot([Link], main="Single Linkage", xlab="", sub="", cex=.9)
cutree([Link], 2)
cutree([Link], 2)
cutree([Link], 2)
cutree([Link], 4)
xsc=scale(x)
plot(hclust(dist(xsc), method="complete"), main="Hierarchical Clustering with Sc
aled Features")
x=matrix(rnorm(30*3), ncol=3)
dd=[Link](1-cor(t(x)))
plot(hclust(dd, method="complete"), main="Complete Linkage with Correlation-Base
d Distance", xlab="", sub="")
# Chapter 10 Lab 3: NCI60 Data Example
# The NCI60 data
library(ISLR)
[Link]=NCI60$labs
[Link]=NCI60$data
dim([Link])
[Link][1:4]
table([Link])
# PCA on the NCI60 Data
[Link]=prcomp([Link], scale=TRUE)
Cols=function(vec){
cols=rainbow(length(unique(vec)))
return(cols[[Link]([Link](vec))])
}
par(mfrow=c(1,2))
plot([Link]$x[,1:2], col=Cols([Link]), pch=19,xlab="Z1",ylab="Z2")
plot([Link]$x[,c(1,3)], col=Cols([Link]), pch=19,xlab="Z1",ylab="Z3")
summary([Link])

plot([Link])
pve=100*[Link]$sdev^2/sum([Link]$sdev^2)
par(mfrow=c(1,2))
plot(pve, type="o", ylab="PVE", xlab="Principal Component", col="blue")
plot(cumsum(pve), type="o", ylab="Cumulative PVE", xlab="Principal Component", c
ol="brown3")
# Clustering the Observations of the NCI60 Data
[Link]=scale([Link])
par(mfrow=c(1,3))
[Link]=dist([Link])
plot(hclust([Link]), labels=[Link], main="Complete Linkage", xlab="", sub="
",ylab="")
plot(hclust([Link], method="average"), labels=[Link], main="Average Linkage
", xlab="", sub="",ylab="")
plot(hclust([Link], method="single"), labels=[Link], main="Single Linkage"
, xlab="", sub="",ylab="")
[Link]=hclust(dist([Link]))
[Link]=cutree([Link],4)
table([Link],[Link])
par(mfrow=c(1,1))
plot([Link], labels=[Link])
abline(h=139, col="red")
[Link]
[Link](2)
[Link]=kmeans([Link], 4, nstart=20)
[Link]=[Link]$cluster
table([Link],[Link])
[Link]=hclust(dist([Link]$x[,1:5]))
plot([Link], labels=[Link], main="Hier. Clust. on First Five Score Vectors")
table(cutree([Link],4), [Link])

You might also like