rm(list=ls())
library(ElemStatLearn)
data(spam)
# names(spam) <- c("make", "address", "all", "3d", "our",
# "over", "remove", "internet", "order", "mail",
# "receive", "will", "people", "report", "addresses",
# "free", "business", "email", "you", "credit",
# "your", "font", "000", "money", "hp",
# "hpl", "george", "650", "lab", "labs",
# "telnet", "857", "data", "415", "85",
# "technology", "1999", "parts", "pm",
# "direct", "cs", "meeting", "original", "project",
# "re", "edu", "table", "conference", ";:",
# "(:", "[:", "!:", "$:", "#:",
# "CRave", "CRlong", "CRtotal", "spam")
names(spam) <- c("make", "address", "all", "x3d", "our",
"over", "remove", "internet", "order", "mail",
"receive", "will", "people", "report", "addresses",
"free", "business", "email", "you", "credit",
"your", "font", "x000", "money", "hp",
"hpl", "george", "x650", "lab", "labs",
"telnet", "x857", "data", "x415", "x85",
"technology", "x1999", "parts", "pm",
"direct", "cs", "meeting", "original", "project",
"re", "edu", "table", "conference", "p1",
"p2", "p3", "p4", "p5", "p6",
"CRave", "CRlong", "CRtotal", "spam")
summary(spam)
spam.test_indx = read.delim("http://www-stat.stanford.edu/~tibs/ElemStatLearn/datasets/spam.traintest",
sep="\n", header=FALSE)
#########################################################################
#########################################################################
Y = as.data.frame(matrix(rep(0,4601),nrow=4601,ncol=1))
Y[spam$spam == "spam", ] = 1
names(Y) = "spam"
Y[,1]=factor(Y[,1])
data = cbind(spam[,-58],Y)
data.train = data[spam.test_indx == 0,]
data.test = data[spam.test_indx == 1,]
summary(data.train)
summary(data.test)
rm(Y, data, spam, spam.test_indx)
###############################################################################
###############################################################################
library(rpart)
library(pROC)
names(data.train)
spam.tree = rpart(spam~.
, data = data.train
, method = "class"
, xval = 5
, cp = 0.00001
, minsplit = 1
, parms=list(split='information')
, na.action = na.exclude
)
printcp(spam.tree)
plotcp(spam.tree, upper = "size")
spam.prune = prune.rpart(spam.tree, 0.0025)
print(spam.prune)
plot(spam.prune)
plot(spam.prune, compress=T, uniform=T, branch=0.4, margin=0.01)
text(spam.prune)
spam.prune = prune.rpart(spam.tree, 0.003)
plot(spam.prune, compress=T, uniform=T, branch=0.4, margin=0.01)
text(spam.prune)
summary(spam.prune)
plotcp(spam.prune)
y.hat = predict(spam.prune, data.test, type="prob")[,2]
head(y.hat)
roc(data.test$spam,y.hat,plot=T)
############################################################################
spam.tree = rpart(spam~.
, data = data.train
, method = "class"
, xval = 5
, cp = 0.00001
, minsplit = 1
, parms=list(split='gini')
, na.action = na.exclude
)
plotcp(spam.tree)
spam.prune = prune.rpart(spam.tree, 0.002)
plotcp(spam.prune)
plot(spam.prune, compress=T, uniform=T, branch=0.4, margin=0.01)
text(spam.prune)
summary(spam.prune)
print(spam.prune)
y.hat = predict(spam.prune, data.test, type="prob")[,2]
roc(data.test$spam,y.hat,plot=T)
#Remark: for cross-entropy, smaller trees get better ROC than Gini
Monday, April 29, 2013
GAM: predict email spam
rm(list=ls())
library(ElemStatLearn)
data(spam)
# names(spam) <- c("make", "address", "all", "3d", "our",
# "over", "remove", "internet", "order", "mail",
# "receive", "will", "people", "report", "addresses",
# "free", "business", "email", "you", "credit",
# "your", "font", "000", "money", "hp",
# "hpl", "george", "650", "lab", "labs",
# "telnet", "857", "data", "415", "85",
# "technology", "1999", "parts", "pm",
# "direct", "cs", "meeting", "original", "project",
# "re", "edu", "table", "conference", ";:",
# "(:", "[:", "!:", "$:", "#:",
# "CRave", "CRlong", "CRtotal", "spam")
names(spam) <- c("make", "address", "all", "x3d", "our",
"over", "remove", "internet", "order", "mail",
"receive", "will", "people", "report", "addresses",
"free", "business", "email", "you", "credit",
"your", "font", "x000", "money", "hp",
"hpl", "george", "x650", "lab", "labs",
"telnet", "x857", "data", "x415", "x85",
"technology", "x1999", "parts", "pm",
"direct", "cs", "meeting", "original", "project",
"re", "edu", "table", "conference", "p1",
"p2", "p3", "p4", "p5", "p6",
"CRave", "CRlong", "CRtotal", "spam")
summary(spam)
spam.test_indx = read.delim("http://www-stat.stanford.edu/~tibs/ElemStatLearn/datasets/spam.traintest",
sep="\n", header=FALSE)
#########################################################################
#########################################################################
Y = as.data.frame(matrix(rep(0,4601),nrow=4601,ncol=1))
Y[spam$spam == "spam", ] = 1
names(Y) = "spam"
Y[,1]=factor(Y[,1])
data = cbind(spam[,-58],Y)
data.train = data[spam.test_indx == 0,]
data.test = data[spam.test_indx == 1,]
summary(data.train)
summary(data.test)
###########################################################################
library(gam)
# var.name = names(data)
# for (i in c(1:(length(var.name)-1))){
# cat("+s(", var.name[i], ", 4)\n", append=TRUE, sep = "", collapse="")
# }
spamgam = gam(spam~
+s(make, 4)
+s(address, 4)
+s(all, 4)
+s(x3d, 4)
+s(our, 4)
+s(over, 4)
+s(remove, 4)
+s(internet, 4)
+s(order, 4)
+s(mail, 4)
+s(receive, 4)
+s(will, 4)
+s(people, 4)
+s(report, 4)
+s(addresses, 4)
+s(free, 4)
+s(business, 4)
+s(email, 4)
+s(you, 4)
+s(credit, 4)
+s(your, 4)
+s(font, 4)
+s(x000, 4)
+s(money, 4)
+s(hp, 4)
+s(hpl, 4)
+s(george, 4)
+s(x650, 4)
+s(lab, 4)
+s(labs, 4)
+s(telnet, 4)
+s(x857, 4)
+s(data, 4)
+s(x415, 4)
+s(x85, 4)
+s(technology, 4)
+s(x1999, 4)
+s(parts, 4)
+s(pm, 4)
+s(direct, 4)
+s(cs, 4)
+s(meeting, 4)
+s(original, 4)
+s(project, 4)
+s(re, 4)
+s(edu, 4)
+s(table, 4)
+s(conference, 4)
+s(p1, 4)
+s(p2, 4)
+s(p3, 4)
+s(p4, 4)
+s(p5, 4)
+s(p6, 4)
+s(CRave, 4)
+s(CRlong, 4)
+s(CRtotal, 4)
,data = data.train
,family = binomial(link = "logit")
)
plot(spamgam, residual=TRUE, se=TRUE,pch=".")
summary(spamgam)
library(pROC)
y.hat = predict(spamgam, data.test, type="response")
roc(data.test$spam,y.hat,plot=T)
############################################################
for (i in c(1:(length(var.name)-1))){
cat("+s(log(", var.name[i], "+0.1))\n", append=TRUE, sep = "", collapse="")
}
spamgam = gam(spam~
+s(log(make+0.1))
+s(log(address+0.1))
+s(log(all+0.1))
+s(log(x3d+0.1))
+s(log(our+0.1))
+s(log(over+0.1))
+s(log(remove+0.1))
+s(log(internet+0.1))
+s(log(order+0.1))
+s(log(mail+0.1))
+s(log(receive+0.1))
+s(log(will+0.1))
+s(log(people+0.1))
+s(log(report+0.1))
+s(log(addresses+0.1))
+s(log(free+0.1))
+s(log(business+0.1))
+s(log(email+0.1))
+s(log(you+0.1))
+s(log(credit+0.1))
+s(log(your+0.1))
+s(log(font+0.1))
+s(log(x000+0.1))
+s(log(money+0.1))
+s(log(hp+0.1))
+s(log(hpl+0.1))
+s(log(george+0.1))
+s(log(x650+0.1))
+s(log(lab+0.1))
+s(log(labs+0.1))
+s(log(telnet+0.1))
+s(log(x857+0.1))
+s(log(data+0.1))
+s(log(x415+0.1))
+s(log(x85+0.1))
+s(log(technology+0.1))
+s(log(x1999+0.1))
+s(log(parts+0.1))
+s(log(pm+0.1))
+s(log(direct+0.1))
+s(log(cs+0.1))
+s(log(meeting+0.1))
+s(log(original+0.1))
+s(log(project+0.1))
+s(log(re+0.1))
+s(log(edu+0.1))
+s(log(table+0.1))
+s(log(conference+0.1))
+s(log(p1+0.1))
+s(log(p2+0.1))
+s(log(p3+0.1))
+s(log(p4+0.1))
+s(log(p5+0.1))
+s(log(p6+0.1))
+s(log(CRave+0.1))
+s(log(CRlong+0.1))
+s(log(CRtotal+0.1))
,data = data.train
,family = binomial
)
plot(spamgam, residual=TRUE, se=TRUE,pch=".")
summary(spamgam)
y.hat = predict(spamgam, data.test, type="response")
# head(y.hat,3)
# head(data.test$spam,3)
library(pROC)
roc(data.test$spam,y.hat,plot=T)
# penalized GAM is availabe in mgcv, but extremely slow.
library(ElemStatLearn)
data(spam)
# names(spam) <- c("make", "address", "all", "3d", "our",
# "over", "remove", "internet", "order", "mail",
# "receive", "will", "people", "report", "addresses",
# "free", "business", "email", "you", "credit",
# "your", "font", "000", "money", "hp",
# "hpl", "george", "650", "lab", "labs",
# "telnet", "857", "data", "415", "85",
# "technology", "1999", "parts", "pm",
# "direct", "cs", "meeting", "original", "project",
# "re", "edu", "table", "conference", ";:",
# "(:", "[:", "!:", "$:", "#:",
# "CRave", "CRlong", "CRtotal", "spam")
names(spam) <- c("make", "address", "all", "x3d", "our",
"over", "remove", "internet", "order", "mail",
"receive", "will", "people", "report", "addresses",
"free", "business", "email", "you", "credit",
"your", "font", "x000", "money", "hp",
"hpl", "george", "x650", "lab", "labs",
"telnet", "x857", "data", "x415", "x85",
"technology", "x1999", "parts", "pm",
"direct", "cs", "meeting", "original", "project",
"re", "edu", "table", "conference", "p1",
"p2", "p3", "p4", "p5", "p6",
"CRave", "CRlong", "CRtotal", "spam")
summary(spam)
spam.test_indx = read.delim("http://www-stat.stanford.edu/~tibs/ElemStatLearn/datasets/spam.traintest",
sep="\n", header=FALSE)
#########################################################################
#########################################################################
Y = as.data.frame(matrix(rep(0,4601),nrow=4601,ncol=1))
Y[spam$spam == "spam", ] = 1
names(Y) = "spam"
Y[,1]=factor(Y[,1])
data = cbind(spam[,-58],Y)
data.train = data[spam.test_indx == 0,]
data.test = data[spam.test_indx == 1,]
summary(data.train)
summary(data.test)
###########################################################################
library(gam)
# var.name = names(data)
# for (i in c(1:(length(var.name)-1))){
# cat("+s(", var.name[i], ", 4)\n", append=TRUE, sep = "", collapse="")
# }
spamgam = gam(spam~
+s(make, 4)
+s(address, 4)
+s(all, 4)
+s(x3d, 4)
+s(our, 4)
+s(over, 4)
+s(remove, 4)
+s(internet, 4)
+s(order, 4)
+s(mail, 4)
+s(receive, 4)
+s(will, 4)
+s(people, 4)
+s(report, 4)
+s(addresses, 4)
+s(free, 4)
+s(business, 4)
+s(email, 4)
+s(you, 4)
+s(credit, 4)
+s(your, 4)
+s(font, 4)
+s(x000, 4)
+s(money, 4)
+s(hp, 4)
+s(hpl, 4)
+s(george, 4)
+s(x650, 4)
+s(lab, 4)
+s(labs, 4)
+s(telnet, 4)
+s(x857, 4)
+s(data, 4)
+s(x415, 4)
+s(x85, 4)
+s(technology, 4)
+s(x1999, 4)
+s(parts, 4)
+s(pm, 4)
+s(direct, 4)
+s(cs, 4)
+s(meeting, 4)
+s(original, 4)
+s(project, 4)
+s(re, 4)
+s(edu, 4)
+s(table, 4)
+s(conference, 4)
+s(p1, 4)
+s(p2, 4)
+s(p3, 4)
+s(p4, 4)
+s(p5, 4)
+s(p6, 4)
+s(CRave, 4)
+s(CRlong, 4)
+s(CRtotal, 4)
,data = data.train
,family = binomial(link = "logit")
)
plot(spamgam, residual=TRUE, se=TRUE,pch=".")
summary(spamgam)
library(pROC)
y.hat = predict(spamgam, data.test, type="response")
roc(data.test$spam,y.hat,plot=T)
############################################################
for (i in c(1:(length(var.name)-1))){
cat("+s(log(", var.name[i], "+0.1))\n", append=TRUE, sep = "", collapse="")
}
spamgam = gam(spam~
+s(log(make+0.1))
+s(log(address+0.1))
+s(log(all+0.1))
+s(log(x3d+0.1))
+s(log(our+0.1))
+s(log(over+0.1))
+s(log(remove+0.1))
+s(log(internet+0.1))
+s(log(order+0.1))
+s(log(mail+0.1))
+s(log(receive+0.1))
+s(log(will+0.1))
+s(log(people+0.1))
+s(log(report+0.1))
+s(log(addresses+0.1))
+s(log(free+0.1))
+s(log(business+0.1))
+s(log(email+0.1))
+s(log(you+0.1))
+s(log(credit+0.1))
+s(log(your+0.1))
+s(log(font+0.1))
+s(log(x000+0.1))
+s(log(money+0.1))
+s(log(hp+0.1))
+s(log(hpl+0.1))
+s(log(george+0.1))
+s(log(x650+0.1))
+s(log(lab+0.1))
+s(log(labs+0.1))
+s(log(telnet+0.1))
+s(log(x857+0.1))
+s(log(data+0.1))
+s(log(x415+0.1))
+s(log(x85+0.1))
+s(log(technology+0.1))
+s(log(x1999+0.1))
+s(log(parts+0.1))
+s(log(pm+0.1))
+s(log(direct+0.1))
+s(log(cs+0.1))
+s(log(meeting+0.1))
+s(log(original+0.1))
+s(log(project+0.1))
+s(log(re+0.1))
+s(log(edu+0.1))
+s(log(table+0.1))
+s(log(conference+0.1))
+s(log(p1+0.1))
+s(log(p2+0.1))
+s(log(p3+0.1))
+s(log(p4+0.1))
+s(log(p5+0.1))
+s(log(p6+0.1))
+s(log(CRave+0.1))
+s(log(CRlong+0.1))
+s(log(CRtotal+0.1))
,data = data.train
,family = binomial
)
plot(spamgam, residual=TRUE, se=TRUE,pch=".")
summary(spamgam)
y.hat = predict(spamgam, data.test, type="response")
# head(y.hat,3)
# head(data.test$spam,3)
library(pROC)
roc(data.test$spam,y.hat,plot=T)
# penalized GAM is availabe in mgcv, but extremely slow.
Sunday, March 31, 2013
Spline Regression
Chapter 5
rm(list=ls())
library(ElemStatLearn)
library(splines)
# Basic Examples
set.seed(0)
x = seq(0, 4*pi, length.out=50)
y = cos(x) + 0.3*rnorm(length(x))
par(mfrow=c(2,3))
# linear, no interior knots
plot(x,y,type="p", main = "Deg = 1, Df = 1")
m1 = lm(y ~ bs(x,degree=1,df=1))
lines(x, fitted(m1))
# linear, 1 interior knot
plot(x,y,type="p", main = "Deg = 1, Df = 2")
m2 = lm(y ~ bs(x,degree=1,df=2))
lines(x, fitted(m2))
# linear, 2 interior knots
plot(x, y,type="p", main = "Deg = 1, Df = 3")
m3 = lm(y ~ bs(x,degree=1,df=3))
lines(x, fitted(m3))
# quadratic, no interior knots
plot(x, y,type="p", main = "Deg = 2, Df = 2")
m4 = lm( y ~ bs(x,degree=2,df=2))
lines(x, fitted(m4))
# quadratic, 1 interior knot
plot(x, y,type="p", main = "Deg = 2, Df = 3")
m5 = lm(y ~ bs(x,degree=2,df=3))
lines(x, fitted(m5))
# quadratic, 2 interior knots
plot(x, y,type="p", main = "Deg = 2, Df = 4")
m6 = lm(y ~ bs(x,degree=2,df=4))
lines(x, fitted(m6))
par(mfrow=c(1,1))
##########################################################################
##########################################################################
################ South African Heart Disease #########################
##########################################################################
##########################################################################
form="chd~ns(sbp,df=4)+ns(tobacco,df=4)+ns(ldl,df=4)+famhist+ns(obesity,df=4)+ns(alcohol,df=4)+ns(age,df=4)"
form = formula(form)
sbp.ns = ns(SAheart$sbp,df=4)
attr(sbp.ns,"knots")
summary(SAheart$sbp)
# Natural Cubic Splines
model.ncs = glm(form, data=SAheart, family=binomial("logit"))
summary(model.ncs)
names(model.ncs)
# model.ncs$coefficients
# # model.ncs$effects
# model.ncs$xlevels
# model.ncs$method
# # model.ncs$data
# model.ncs$terms # knots
# # model.ncs$model
# recreate Table 5.1
drop1( model.ncs, scope=form, test="Chisq" )
# create the model in Table 5.1
library(MASS)
stepAIC(model.ncs, scope=form, k=2)
##########################################################################
##########################################################################
####################### Phoneme Recognition ##########################
############ Example: Filtering and Feature Extraction #################
##########################################################################
##########################################################################
rm(list=ls())
data(phoneme)
set.seed(0)
aa_indx = which(phoneme$g == "aa")
ao_indx = which(phoneme$g == "ao")
AA_data = phoneme[aa_indx[1:15], 1:256]
AO_data = phoneme[ao_indx[1:15], 1:256]
min_l = min(c(min(AA_data), min(AO_data)))
max_l = max(c(max(AA_data), max(AO_data)))
# Figure 5.5 on page 149
ii=1
plot( as.double(AA_data[ii,]),
ylim=c(min_l,max_l), type="l", col="red",
xlab="Frequency", ylab="Log-periodogram" )
for( ii in 2:dim(AA_data)[1] ){
lines( as.double(AA_data[ii,]), col="red" )
}
for( ii in 1:dim(AO_data)[1] ){
lines( as.double(AO_data[ii,]), col="blue" )
}
rm(list=ls())
aa_indx = which(phoneme$g == "aa")
ao_indx = which(phoneme$g == "ao")
###########################################################
AA_data = phoneme[aa_indx,]
train_inds = grep("^train", AA_data$speaker)
test_inds = grep("^test", AA_data$speaker)
AA_data_train = AA_data[train_inds, 1:256]
AA_data_test = AA_data[test_inds, 1:256]
#####################################
###### call AA class 1 #############
#####################################
AA_data_train$Y = rep(1, dim(AA_data_train)[1])
AA_data_test$Y = rep(1, dim(AA_data_test)[1])
###########################################################
AO_data = phoneme[ao_indx,]
train_inds = grep("^train", AO_data$speaker)
test_inds = grep("^test", AO_data$speaker)
AO_data_train = AO_data[train_inds, 1:256]
AO_data_test = AO_data[test_inds, 1:256]
#####################################
###### call AO class 1 #############
#####################################
AO_data_train$Y = rep(0, dim(AO_data_train)[1])
AO_data_test$Y = rep(0, dim(AO_data_test)[1])
####################################################
################## Combine Data ##################
####################################################
Data_train = rbind(AA_data_train, AO_data_train)
Data_test = rbind(AA_data_test, AO_data_test)
#########################################################################
############# Logistic Regression Without Intercept ##################
#########################################################################
form = paste("Y ~ -1 + ", paste(colnames(Data_train)[1:256], collapse="+"))
m = glm(form, family=binomial("logit"), data=Data_train)
mc = as.double(coefficients(m))
m.smooth = smooth.spline(mc,df=15)
plot(mc, ylim=c(-0.8,+0.6),type="l",
xlab="Frequency",ylab="Logistic Regression Coefficients")
lines(m.smooth, col="blue" )
abline(h=0)
Y_hat_train = predict(m, Data_train[,1:256], type="response")
predicted_class_label_train = as.double(Y_hat_train > 0.5)
Y_hat_test = predict(m, Data_test[,1:256], type="response")
predicted_class_label_test = as.double(Y_hat_test > 0.5)
err_rate_train1 = mean(Data_train[,257] != predicted_class_label_train)
err_rate_test1 = mean(Data_test[,257] != predicted_class_label_test)
print(sprintf('Logistic Regression (Raw): err_rate_train= %10.6f; err_rate_test= %10.5f',
err_rate_train1, err_rate_test1))
beta = matrix(m.smooth$y,ncol=1)
logit_train = as.matrix(Data_train[,1:256])%*%beta
Y_smooth_train = exp(logit_train)/(1+exp(logit_train))
predicted_class_smooth_train = as.double(Y_smooth_train > 0.5)
logit_test = as.matrix(Data_test[,1:256])%*%beta
Y_smooth_test = exp(logit_test)/(1+exp(logit_test))
predicted_class_smooth_test = as.double(Y_smooth_test > 0.5)
err_rate_train2 = mean(Data_train[,257] != predicted_class_smooth_train)
err_rate_test2 = mean(Data_test[,257] != predicted_class_smooth_test)
print(sprintf('Logistic Regression (Raw): err_rate_train= %10.6f; err_rate_test= %10.5f',
err_rate_train1, err_rate_test1))
print(sprintf('Logistic Regression (Regularized): err_rate_train= %10.6f; err_rate_test= %10.5f',
err_rate_train2, err_rate_test2))
#######################################################################
############# BMD Example ############################################
#######################################################################
males = bone$gender == "male"
females = bone$gender == "female"
# df is monotone in lambda for smoothing splines
boneMaleSmooth = smooth.spline( bone[bone$gender == "male","age"],
bone[bone$gender == "male","spnbmd"], df=12)
boneFemaleSmooth = smooth.spline( bone[bone$gender == "female","age"],
bone[bone$gender == "female","spnbmd"], df=12)
plot(boneMaleSmooth, ylim=c(-0.05,0.20), col="blue", type="l", xlab="Age", ylab="spnbmd")
points(bone[males,c(2,4)], col="blue", pch=20)
lines(boneFemaleSmooth, ylim=c(-0.05,0.20), col="red")
points(bone[females,c(2,4)], col="red", pch=20)
rm(list=ls())
library(ElemStatLearn)
library(splines)
# Basic Examples
set.seed(0)
x = seq(0, 4*pi, length.out=50)
y = cos(x) + 0.3*rnorm(length(x))
par(mfrow=c(2,3))
# linear, no interior knots
plot(x,y,type="p", main = "Deg = 1, Df = 1")
m1 = lm(y ~ bs(x,degree=1,df=1))
lines(x, fitted(m1))
# linear, 1 interior knot
plot(x,y,type="p", main = "Deg = 1, Df = 2")
m2 = lm(y ~ bs(x,degree=1,df=2))
lines(x, fitted(m2))
# linear, 2 interior knots
plot(x, y,type="p", main = "Deg = 1, Df = 3")
m3 = lm(y ~ bs(x,degree=1,df=3))
lines(x, fitted(m3))
# quadratic, no interior knots
plot(x, y,type="p", main = "Deg = 2, Df = 2")
m4 = lm( y ~ bs(x,degree=2,df=2))
lines(x, fitted(m4))
# quadratic, 1 interior knot
plot(x, y,type="p", main = "Deg = 2, Df = 3")
m5 = lm(y ~ bs(x,degree=2,df=3))
lines(x, fitted(m5))
# quadratic, 2 interior knots
plot(x, y,type="p", main = "Deg = 2, Df = 4")
m6 = lm(y ~ bs(x,degree=2,df=4))
lines(x, fitted(m6))
par(mfrow=c(1,1))
##########################################################################
##########################################################################
################ South African Heart Disease #########################
##########################################################################
##########################################################################
form="chd~ns(sbp,df=4)+ns(tobacco,df=4)+ns(ldl,df=4)+famhist+ns(obesity,df=4)+ns(alcohol,df=4)+ns(age,df=4)"
form = formula(form)
sbp.ns = ns(SAheart$sbp,df=4)
attr(sbp.ns,"knots")
summary(SAheart$sbp)
# Natural Cubic Splines
model.ncs = glm(form, data=SAheart, family=binomial("logit"))
summary(model.ncs)
names(model.ncs)
# model.ncs$coefficients
# # model.ncs$effects
# model.ncs$xlevels
# model.ncs$method
# # model.ncs$data
# model.ncs$terms # knots
# # model.ncs$model
# recreate Table 5.1
drop1( model.ncs, scope=form, test="Chisq" )
# create the model in Table 5.1
library(MASS)
stepAIC(model.ncs, scope=form, k=2)
##########################################################################
##########################################################################
####################### Phoneme Recognition ##########################
############ Example: Filtering and Feature Extraction #################
##########################################################################
##########################################################################
rm(list=ls())
data(phoneme)
set.seed(0)
aa_indx = which(phoneme$g == "aa")
ao_indx = which(phoneme$g == "ao")
AA_data = phoneme[aa_indx[1:15], 1:256]
AO_data = phoneme[ao_indx[1:15], 1:256]
min_l = min(c(min(AA_data), min(AO_data)))
max_l = max(c(max(AA_data), max(AO_data)))
# Figure 5.5 on page 149
ii=1
plot( as.double(AA_data[ii,]),
ylim=c(min_l,max_l), type="l", col="red",
xlab="Frequency", ylab="Log-periodogram" )
for( ii in 2:dim(AA_data)[1] ){
lines( as.double(AA_data[ii,]), col="red" )
}
for( ii in 1:dim(AO_data)[1] ){
lines( as.double(AO_data[ii,]), col="blue" )
}
rm(list=ls())
aa_indx = which(phoneme$g == "aa")
ao_indx = which(phoneme$g == "ao")
###########################################################
AA_data = phoneme[aa_indx,]
train_inds = grep("^train", AA_data$speaker)
test_inds = grep("^test", AA_data$speaker)
AA_data_train = AA_data[train_inds, 1:256]
AA_data_test = AA_data[test_inds, 1:256]
#####################################
###### call AA class 1 #############
#####################################
AA_data_train$Y = rep(1, dim(AA_data_train)[1])
AA_data_test$Y = rep(1, dim(AA_data_test)[1])
###########################################################
AO_data = phoneme[ao_indx,]
train_inds = grep("^train", AO_data$speaker)
test_inds = grep("^test", AO_data$speaker)
AO_data_train = AO_data[train_inds, 1:256]
AO_data_test = AO_data[test_inds, 1:256]
#####################################
###### call AO class 1 #############
#####################################
AO_data_train$Y = rep(0, dim(AO_data_train)[1])
AO_data_test$Y = rep(0, dim(AO_data_test)[1])
####################################################
################## Combine Data ##################
####################################################
Data_train = rbind(AA_data_train, AO_data_train)
Data_test = rbind(AA_data_test, AO_data_test)
#########################################################################
############# Logistic Regression Without Intercept ##################
#########################################################################
form = paste("Y ~ -1 + ", paste(colnames(Data_train)[1:256], collapse="+"))
m = glm(form, family=binomial("logit"), data=Data_train)
mc = as.double(coefficients(m))
m.smooth = smooth.spline(mc,df=15)
plot(mc, ylim=c(-0.8,+0.6),type="l",
xlab="Frequency",ylab="Logistic Regression Coefficients")
lines(m.smooth, col="blue" )
abline(h=0)
Y_hat_train = predict(m, Data_train[,1:256], type="response")
predicted_class_label_train = as.double(Y_hat_train > 0.5)
Y_hat_test = predict(m, Data_test[,1:256], type="response")
predicted_class_label_test = as.double(Y_hat_test > 0.5)
err_rate_train1 = mean(Data_train[,257] != predicted_class_label_train)
err_rate_test1 = mean(Data_test[,257] != predicted_class_label_test)
print(sprintf('Logistic Regression (Raw): err_rate_train= %10.6f; err_rate_test= %10.5f',
err_rate_train1, err_rate_test1))
beta = matrix(m.smooth$y,ncol=1)
logit_train = as.matrix(Data_train[,1:256])%*%beta
Y_smooth_train = exp(logit_train)/(1+exp(logit_train))
predicted_class_smooth_train = as.double(Y_smooth_train > 0.5)
logit_test = as.matrix(Data_test[,1:256])%*%beta
Y_smooth_test = exp(logit_test)/(1+exp(logit_test))
predicted_class_smooth_test = as.double(Y_smooth_test > 0.5)
err_rate_train2 = mean(Data_train[,257] != predicted_class_smooth_train)
err_rate_test2 = mean(Data_test[,257] != predicted_class_smooth_test)
print(sprintf('Logistic Regression (Raw): err_rate_train= %10.6f; err_rate_test= %10.5f',
err_rate_train1, err_rate_test1))
print(sprintf('Logistic Regression (Regularized): err_rate_train= %10.6f; err_rate_test= %10.5f',
err_rate_train2, err_rate_test2))
#######################################################################
############# BMD Example ############################################
#######################################################################
males = bone$gender == "male"
females = bone$gender == "female"
# df is monotone in lambda for smoothing splines
boneMaleSmooth = smooth.spline( bone[bone$gender == "male","age"],
bone[bone$gender == "male","spnbmd"], df=12)
boneFemaleSmooth = smooth.spline( bone[bone$gender == "female","age"],
bone[bone$gender == "female","spnbmd"], df=12)
plot(boneMaleSmooth, ylim=c(-0.05,0.20), col="blue", type="l", xlab="Age", ylab="spnbmd")
points(bone[males,c(2,4)], col="blue", pch=20)
lines(boneFemaleSmooth, ylim=c(-0.05,0.20), col="red")
points(bone[females,c(2,4)], col="red", pch=20)
Tuesday, March 26, 2013
Simple LDA and QDA in R
rm(list=ls())
library(MASS)
library(ElemStatLearn)
data(mixture.example)
data = mixture.example
rm(mixture.example)
x = data$x
y = data$y
summary(factor(y))
rm(data)
plot(x,col = 2*y + 1, xlab = "X1", ylab = "X2")
# LDA
model.lda = lda(x, y)
summary(model.lda)
model.lda
model.lda$prior
model.lda$means
model.lda$scaling
abline(model.lda$scaling)
points(model.lda$means, pch = 20)
# QDA
model.qda = qda(x, y)
model.qda
scaling.matrix1 = model.qda$scaling[,,1]
scaling.matrix2 = model.qda$scaling[,,2]
names(model.qda)
model2.qda = qda(x, y, CV = TRUE)
model2.qda
names(model2.qda)
model2.qda$class
library(MASS)
library(ElemStatLearn)
data(mixture.example)
data = mixture.example
rm(mixture.example)
x = data$x
y = data$y
summary(factor(y))
rm(data)
plot(x,col = 2*y + 1, xlab = "X1", ylab = "X2")
# LDA
model.lda = lda(x, y)
summary(model.lda)
model.lda
model.lda$prior
model.lda$means
model.lda$scaling
abline(model.lda$scaling)
points(model.lda$means, pch = 20)
# QDA
model.qda = qda(x, y)
model.qda
scaling.matrix1 = model.qda$scaling[,,1]
scaling.matrix2 = model.qda$scaling[,,2]
names(model.qda)
model2.qda = qda(x, y, CV = TRUE)
model2.qda
names(model2.qda)
model2.qda$class
Monday, March 25, 2013
Logistic and LASSO Regression
# The following code is for the book The Elements of Statistical Learning,
chapter 4
# Example: South African Heart Disease (Page: 122)
# load data
rm(list=ls())
library(ElemStatLearn)
data(SAheart)
data = SAheart[,c(1:3,5,7:10)]
# to convert factor variables into dummy variables
temp = matrix(0,nrow(data),1)
for (i in c(1:nrow(data))){
if (data[i,4] == 'Present'){temp[i] = 1}
else {temp[i] = 0}
}
temp = as.data.frame(temp)
names(temp) = "famhist"
data.new = cbind(data[,1:3], temp, data[,5:8])
# in order to apply glmnet function
# change dataframe to matrix
Y = as.matrix(data.new[,8])
X = as.matrix(data.new[,-8])
logit = glm(chd ~ ., family=binomial("logit"),data=data.new,na.action=na.exclude)
# table 4.2 on page 122
summary(logit)
# Figure 4.12
pairs(data.new, main="South African heart disease data",pch = 21)
# reduced model: stepwise (same as the linear regression)
logit.reduced = glm(chd ~ tobacco + ldl + famhist + age,
family=binomial("logit"),data=data.new,na.action=na.exclude)
# table 4.3 on page 124
summary(logit.reduced)
library(glmnet)
# not work if X has a factor variable
# This is the only package I know to apply lasso for binormial response
# work
lasso = glmnet(scale(X), Y, family = "binomial", alpha = 1, standardize = FALSE, intercept=TRUE)
# not work
lasso = glmnet(X, Y, family = "binomial", alpha = 1, standardize = TRUE, intercept=TRUE)
# FIGURE 4.13 on page 126
plot(lasso, xvar = "norm", label=TRUE)
abline(h=0)
# Example: South African Heart Disease (Page: 122)
# load data
rm(list=ls())
library(ElemStatLearn)
data(SAheart)
data = SAheart[,c(1:3,5,7:10)]
# to convert factor variables into dummy variables
temp = matrix(0,nrow(data),1)
for (i in c(1:nrow(data))){
if (data[i,4] == 'Present'){temp[i] = 1}
else {temp[i] = 0}
}
temp = as.data.frame(temp)
names(temp) = "famhist"
data.new = cbind(data[,1:3], temp, data[,5:8])
# in order to apply glmnet function
# change dataframe to matrix
Y = as.matrix(data.new[,8])
X = as.matrix(data.new[,-8])
logit = glm(chd ~ ., family=binomial("logit"),data=data.new,na.action=na.exclude)
# table 4.2 on page 122
summary(logit)
# Figure 4.12
pairs(data.new, main="South African heart disease data",pch = 21)
# reduced model: stepwise (same as the linear regression)
logit.reduced = glm(chd ~ tobacco + ldl + famhist + age,
family=binomial("logit"),data=data.new,na.action=na.exclude)
# table 4.3 on page 124
summary(logit.reduced)
library(glmnet)
# not work if X has a factor variable
# This is the only package I know to apply lasso for binormial response
# work
lasso = glmnet(scale(X), Y, family = "binomial", alpha = 1, standardize = FALSE, intercept=TRUE)
# not work
lasso = glmnet(X, Y, family = "binomial", alpha = 1, standardize = TRUE, intercept=TRUE)
# FIGURE 4.13 on page 126
plot(lasso, xvar = "norm", label=TRUE)
abline(h=0)
Regression IV: Principal Components Regression
rm(list=ls())
library(ElemStatLearn)
data(prostate)
data = prostate[,1:9]
train = subset(prostate, train==T)[,1:9]
test = subset(prostate, train!=T)[,1:9]
# scale function in R substract mean,
# and divide by std dev (with df = N-1)
train.s = scale(train)
# fit ridge regression manually
Y = as.numeric(train.s[,9])
X = as.matrix(train.s[,-9])
x.svd = svd(X)
#scatterplots of samples PCs
par(mar=c(1,1,1,1))
layout(matrix(1:64,8,8))
mycols = rainbow(length(Y))
orY = order(Y)
for(i in 1:8){
for(j in 1:8){
plot(x.svd$u[,i],x.svd$u[,j],type="p",pch=16,col=mycols[orY])
}
}
# amount of variance explained
var = 0
var.cum = 0;
for(i in 1:8){
var[i] = x.svd$d[i]/sum(x.svd$d)
var.cum[i] = sum(var)
}
par(mfrow=c(1,2))
par(mar=c(5,4,4,2))
barplot(var,ylab="Amount of Var Explained",xlab="PCs")
barplot(var.cum,ylab="Cummulative Var Explained",xlab="PCs")
# PC direction weights
par(mfrow=c(3,3))
par(mar=c(5,4,3,2))
for(i in 1:8){
barplot(x.svd$v[,i],names.arg=names(train)[1:8])
}
#PC regression
beta.pcr = diag(1/x.svd$d)%*%t(x.svd$u)%*%Y
Z = X%*%x.svd$v
#training error
err.pcr = 0
for(i in 1:8){
err.pcr[i] = sum(as.numeric(Y - Z[,1:i,drop=FALSE]%*%beta.pcr[1:i,1])^2)/length(Y)
}
err.pcr
library(ElemStatLearn)
data(prostate)
data = prostate[,1:9]
train = subset(prostate, train==T)[,1:9]
test = subset(prostate, train!=T)[,1:9]
# scale function in R substract mean,
# and divide by std dev (with df = N-1)
train.s = scale(train)
# fit ridge regression manually
Y = as.numeric(train.s[,9])
X = as.matrix(train.s[,-9])
x.svd = svd(X)
#scatterplots of samples PCs
par(mar=c(1,1,1,1))
layout(matrix(1:64,8,8))
mycols = rainbow(length(Y))
orY = order(Y)
for(i in 1:8){
for(j in 1:8){
plot(x.svd$u[,i],x.svd$u[,j],type="p",pch=16,col=mycols[orY])
}
}
# amount of variance explained
var = 0
var.cum = 0;
for(i in 1:8){
var[i] = x.svd$d[i]/sum(x.svd$d)
var.cum[i] = sum(var)
}
par(mfrow=c(1,2))
par(mar=c(5,4,4,2))
barplot(var,ylab="Amount of Var Explained",xlab="PCs")
barplot(var.cum,ylab="Cummulative Var Explained",xlab="PCs")
# PC direction weights
par(mfrow=c(3,3))
par(mar=c(5,4,3,2))
for(i in 1:8){
barplot(x.svd$v[,i],names.arg=names(train)[1:8])
}
#PC regression
beta.pcr = diag(1/x.svd$d)%*%t(x.svd$u)%*%Y
Z = X%*%x.svd$v
#training error
err.pcr = 0
for(i in 1:8){
err.pcr[i] = sum(as.numeric(Y - Z[,1:i,drop=FALSE]%*%beta.pcr[1:i,1])^2)/length(Y)
}
err.pcr
Regression III: LASSO
rm(list=ls())
library(ElemStatLearn)
data(prostate)
data = prostate[,1:9]
train = subset(prostate, train==T)[,1:9]
test = subset(prostate, train!=T)[,1:9]
# scale function in R substract mean,
# and divide by std dev (with df = N-1)
train.s = scale(train)
center = attributes(train.s)$'scaled:center'
scale = attributes(train.s)$'scaled:scale'
Y = as.numeric(train.s[,9])
X = as.matrix(train.s[,-9])
# scale testing data based on center and scale of training data
test.s = t((t(test) - center)/scale)
Y1 = as.numeric(test.s[,9])
X1 = as.matrix(test.s[,-9])
# unscaled data
Y0 = as.numeric(train[,9])
X0 = as.matrix(train[,-9])
Y01 = as.numeric(test[,9])
X01 = as.matrix(test[,-9])
#################################################################################
########### fit LASSO regression by glmnet package #############################
#################################################################################
library(glmnet)
# Way 1
lasso.glm = glmnet(X,Y,
family = "gaussian",
alpha = 1,
standardize = FALSE,
intercept = FALSE,
standardize.response = FALSE
)
# # Way 2: NOT WORK
lasso.glm = glmnet(X0,Y0,
family = "gaussian",
alpha = 1,
standardize = TRUE,
intercept = TRUE,
standardize.response = TRUE
)
# # Way 3: NOT WORK
lasso.glm = glmnet(X0,Y0,
family = "gaussian",
alpha = 1
)
# Based on my experiments, I have to scale data first
names(lasso.glm)
plot(lasso.glm, xvar = "norm", label=TRUE)
abline(h=0)
lasso.glm$a0
lasso.glm$beta
lasso.glm$lambda
y = predict(lasso.glm, X, type="link")
y1 = predict(lasso.glm, X1, type="link")
lasso.mse = matrix(0, length(lasso.glm$lambda), 1)
lasso1.mse = matrix(0, length(lasso.glm$lambda), 1)
for (i in 1:length(lasso.glm$lambda)){
lasso.mse[i] = mean((Y - y[,i])^2)
lasso1.mse[i] = mean((Y1 - y1[,i])^2)
}
############################################################################
################### Figure 3.10 on Page 70 #############################
############################################################################
plot(lasso.glm$lambda, lasso.mse)
lines(lasso.glm$lambda, lasso1.mse)
# cross-validation
lasso.cv = cv.glmnet(X,Y,
family = "gaussian",
alpha = 1,
nfolds = 10,
standardize = FALSE,
intercept = FALSE,
standardize.response = FALSE,
type.measure = "mse"
)
names(lasso.cv)
plot(lasso.cv)
plot(lasso.cv, -1)
lasso.cv$lambda.1se
lasso.cv$cvm
lasso.cv$cvsd
lasso.cv$glmnet.fit
y.cv = predict(lasso.cv, X, s="lambda.1se")
###########################################################################
############# fit LASSO regression by lars package ###############################
###########################################################################
library(lars)
# cross-validation is available in this package
# scaled version
model.ridge = lars(X,Y,type="lasso",
trace=TRUE,
normalize=FALSE,
intercept=FALSE
)
# unscaled version
model.ridge = lars(X0,Y0,type="lasso",
trace=TRUE,
normalize=TRUE,
intercept=TRUE
)
summary(model.ridge)
model.ridge
plot(model.ridge,
xvar="arc.length",
breaks=TRUE,
plottype="coefficients"
)
predict(model.ridge,
X.test,
type="fit",
mode="step")
coef(model.ridge, s=5.6,mode="step")
# always scale X and Y manually before applying functions
library(ElemStatLearn)
data(prostate)
data = prostate[,1:9]
train = subset(prostate, train==T)[,1:9]
test = subset(prostate, train!=T)[,1:9]
# scale function in R substract mean,
# and divide by std dev (with df = N-1)
train.s = scale(train)
center = attributes(train.s)$'scaled:center'
scale = attributes(train.s)$'scaled:scale'
Y = as.numeric(train.s[,9])
X = as.matrix(train.s[,-9])
# scale testing data based on center and scale of training data
test.s = t((t(test) - center)/scale)
Y1 = as.numeric(test.s[,9])
X1 = as.matrix(test.s[,-9])
# unscaled data
Y0 = as.numeric(train[,9])
X0 = as.matrix(train[,-9])
Y01 = as.numeric(test[,9])
X01 = as.matrix(test[,-9])
#################################################################################
########### fit LASSO regression by glmnet package #############################
#################################################################################
library(glmnet)
# Way 1
lasso.glm = glmnet(X,Y,
family = "gaussian",
alpha = 1,
standardize = FALSE,
intercept = FALSE,
standardize.response = FALSE
)
# # Way 2: NOT WORK
lasso.glm = glmnet(X0,Y0,
family = "gaussian",
alpha = 1,
standardize = TRUE,
intercept = TRUE,
standardize.response = TRUE
)
# # Way 3: NOT WORK
lasso.glm = glmnet(X0,Y0,
family = "gaussian",
alpha = 1
)
# Based on my experiments, I have to scale data first
names(lasso.glm)
plot(lasso.glm, xvar = "norm", label=TRUE)
abline(h=0)
lasso.glm$a0
lasso.glm$beta
lasso.glm$lambda
y = predict(lasso.glm, X, type="link")
y1 = predict(lasso.glm, X1, type="link")
lasso.mse = matrix(0, length(lasso.glm$lambda), 1)
lasso1.mse = matrix(0, length(lasso.glm$lambda), 1)
for (i in 1:length(lasso.glm$lambda)){
lasso.mse[i] = mean((Y - y[,i])^2)
lasso1.mse[i] = mean((Y1 - y1[,i])^2)
}
############################################################################
################### Figure 3.10 on Page 70 #############################
############################################################################
plot(lasso.glm$lambda, lasso.mse)
lines(lasso.glm$lambda, lasso1.mse)
# cross-validation
lasso.cv = cv.glmnet(X,Y,
family = "gaussian",
alpha = 1,
nfolds = 10,
standardize = FALSE,
intercept = FALSE,
standardize.response = FALSE,
type.measure = "mse"
)
names(lasso.cv)
plot(lasso.cv)
plot(lasso.cv, -1)
lasso.cv$lambda.1se
lasso.cv$cvm
lasso.cv$cvsd
lasso.cv$glmnet.fit
y.cv = predict(lasso.cv, X, s="lambda.1se")
###########################################################################
############# fit LASSO regression by lars package ###############################
###########################################################################
library(lars)
# cross-validation is available in this package
# scaled version
model.ridge = lars(X,Y,type="lasso",
trace=TRUE,
normalize=FALSE,
intercept=FALSE
)
# unscaled version
model.ridge = lars(X0,Y0,type="lasso",
trace=TRUE,
normalize=TRUE,
intercept=TRUE
)
summary(model.ridge)
model.ridge
plot(model.ridge,
xvar="arc.length",
breaks=TRUE,
plottype="coefficients"
)
predict(model.ridge,
X.test,
type="fit",
mode="step")
coef(model.ridge, s=5.6,mode="step")
# always scale X and Y manually before applying functions
Subscribe to:
Posts (Atom)