{"id":274,"date":"2013-07-07T13:06:37","date_gmt":"2013-07-07T09:36:37","guid":{"rendered":"http:\/\/vua.nadiran.com\/?p=274"},"modified":"2013-07-10T14:52:45","modified_gmt":"2013-07-10T11:22:45","slug":"%d9%86%d9%85%d9%88%d9%86%d9%87-%d9%87%d8%a7%db%8c-%d8%a8%d8%b1%d9%86%d8%a7%d9%85%d9%87-lasso-%d8%af%d8%b1-%d9%be%da%a9%db%8c%da%86-r","status":"publish","type":"post","link":"https:\/\/vua.nadiran.com\/?p=274","title":{"rendered":"\u0646\u0645\u0648\u0646\u0647 \u0647\u0627\u06cc \u0628\u0631\u0646\u0627\u0645\u0647 \u0631\u06af\u0631\u0633\u06cc\u0648\u0646 \u062f\u0631 \u067e\u06a9\u06cc\u0686 R"},"content":{"rendered":"

\n

## Practical session: Linear regression<\/p>\n

## Jean-Philippe.Vert@mines.org<\/p>\n

##<\/p>\n

## In this practical session to test several linear regression methods on the prostate cancer dataset<\/p>\n

##<\/p>\n

####################################<\/p>\n

## Prepare data<\/p>\n

####################################<\/p>\n

# Download prostate data<\/p>\n

con = url (“http:\/\/www-stat.stanford.edu\/~tibs\/ElemStatLearn\/datasets\/prostate.data”)<\/p>\n

prost=read.csv(con,row.names=1,sep=”\\t”)<\/p>\n

# Alternatively, load the file and read from local file as follows<\/p>\n

# prost=read.csv(‘prostate.data.txt’,row.names=1,sep=”\\t”)<\/p>\n

# Scale data and prepare train\/test split<\/p>\n

prost.std <- data.frame(cbind(scale(prost[,1:8]),prost$lpsa))<\/p>\n

names(prost.std)[9] <- ‘lpsa’<\/p>\n

data.train <- prost.std[prost$train,]<\/p>\n

data.test <- prost.std[!prost$train,]<\/p>\n

y.test <- data.test$lpsa<\/p>\n

n.train <- nrow(data.train)<\/p>\n

\n

####################################<\/p>\n

## Ordinary least squares<\/p>\n

####################################<\/p>\n

m.ols <- lm(lpsa ~ . , data=data.train)<\/p>\n

summary(m.ols)<\/p>\n

y.pred.ols <- predict(m.ols,data.test)<\/p>\n

summary((y.pred.ols – y.test)^2)<\/p>\n

\n

# the importance of a feature depends on other features. Compare the following. Explain (visualize).<\/p>\n

summary(lm(lpsa~svi,data=datatrain))<\/p>\n

summary(lm(lpsa~svi+lcavol,data=datatrain))<\/p>\n

\n

####################################<\/p>\n

## Best subset selection<\/p>\n

####################################<\/p>\n

library(leaps)<\/p>\n

l <- leaps(data.train[,1:8],data.train[,9],method=’r2′)<\/p>\n

plot(l$size,l$r2)<\/p>\n

l <- leaps(data.train[,1:8],data.train[,9],method=’Cp’)<\/p>\n

plot(l$size,l$Cp)<\/p>\n

# Select best model according to Cp<\/p>\n

bestfeat <- l$which[which.min(l$Cp),]<\/p>\n

# Train and test the model on the best subset<\/p>\n

m.bestsubset <- lm(lpsa ~ .,data=data.train[,bestfeat])<\/p>\n

summary(m.bestsubset)<\/p>\n

y.pred.bestsubset <- predict(m.bestsubset,data.test[,bestfeat])<\/p>\n

summary((y.pred.bestsubset – y.test)^2)<\/p>\n

####################################<\/p>\n

## Greedy subset selection<\/p>\n

####################################<\/p>\n

library(MASS)<\/p>\n

# Forward selection<\/p>\n

m.init <- lm(lpsa ~ 1 , data=data.train)<\/p>\n

m.forward <- stepAIC(m.init, scope=list(upper=~lcavol+lweight+age+lbph+svi+lcp+gleason+pgg45,lower=~1) , direction=”forward”)<\/p>\n

y.pred.forward <- predict(m.forward,data.test)<\/p>\n

summary((y.pred.forward – y.test)^2)<\/p>\n

# Backward selection<\/p>\n

m.init <- lm(lpsa ~ . , data=data.train)<\/p>\n

m.backward <- stepAIC(m.init , direction=”backward”)<\/p>\n

y.pred.backward <- predict(m.backward,data.test)<\/p>\n

summary((y.pred.backward – y.test)^2)<\/p>\n

# Hybrid selection<\/p>\n

m.init <- lm(lpsa ~ 1 , data=data.train)<\/p>\n

m.hybrid <- stepAIC(m.init, scope=list(upper=~lcavol+lweight+age+lbph+svi+lcp+gleason+pgg45,lower=~1) , direction=”both”)<\/p>\n

y.pred.hybrid <- predict(m.hybrid,data.test)<\/p>\n

summary((y.pred.hybrid – y.test)^2)<\/p>\n

####################################<\/p>\n

## Ridge regression<\/p>\n

####################################<\/p>\n

library(MASS)<\/p>\n

m.ridge <- lm.ridge(lpsa ~ .,data=data.train, lambda = seq(0,20,0.1))<\/p>\n

plot(m.ridge)<\/p>\n

# select parameter by minimum GCV<\/p>\n

plot(m.ridge$GCV)<\/p>\n

# Predict is not implemented so we need to do it ourselves<\/p>\n

y.pred.ridge = scale(data.test[,1:8],center = F, scale = m.ridge$scales)%*% m.ridge$coef[,which.min(m.ridge$GCV)] + m.ridge$ym<\/p>\n

summary((y.pred.ridge – y.test)^2)<\/p>\n

####################################<\/p>\n

## Lasso<\/p>\n

####################################<\/p>\n

library(lars)<\/p>\n

m.lasso <- lars(as.matrix(data.train[,1:8]),data.train$lpsa)<\/p>\n

plot(m.lasso)<\/p>\n

# Cross-validation<\/p>\n

r <- cv.lars(as.matrix(data.train[,1:8]),data.train$lpsa)<\/p>\n

bestfraction <- r$fraction[which.min(r$cv)]<\/p>\n

##### Note 5\/8\/11: in the new lars package 0.9-8, the field r$fraction seems to have been replaced by r$index. The previous line should therefore be replaced by:<\/p>\n

# bestfraction <- r$index[which.min(r$cv)]<\/p>\n

# Observe coefficients<\/p>\n

coef.lasso <- predict(m.lasso,as.matrix(data.test[,1:8]),s=bestfraction,type=”coefficient”,mode=”fraction”)<\/p>\n

coef.lasso<\/p>\n

# Prediction<\/p>\n

y.pred.lasso <- predict(m.lasso,as.matrix(data.test[,1:8]),s=bestfraction,type=”fit”,mode=”fraction”)$fit<\/p>\n

summary((y.pred.lasso – y.test)^2)<\/p>\n

\n

####################################<\/p>\n

## PCR<\/p>\n

####################################<\/p>\n

library(pls)<\/p>\n

m.pcr <- pcr(lpsa ~ .,data=data.train , validation=”CV”)<\/p>\n

# select number of components (by CV)<\/p>\n

ncomp <- which.min(m.pcr$validation$adj)<\/p>\n

# predict<\/p>\n

y.pred.pcr <- predict(m.pcr,data.test , ncomp=ncomp)<\/p>\n

summary((y.pred.pcr – y.test)^2)<\/p>\n

\n

####################################<\/p>\n

## PLS<\/p>\n

####################################<\/p>\n

library(pls)<\/p>\n

m.pls <- plsr(lpsa ~ .,data=data.train , validation=”CV”)<\/p>\n

# select number of components (by CV)<\/p>\n

ncomp <- which.min(m.pls$validation$adj)<\/p>\n

# predict<\/p>\n

y.pred.pcr <- predict(m.pls,data.test , ncomp=ncomp)<\/p>\n

summary((y.pred.pcr – y.test)^2)<\/p>\n

\n

\r\n\r\n\r\n\r\n<\/pre>\n
\n
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<\/p>\n
%<\/p>\n
% File: Lasso.m<\/p>\n
% Author: Jinzhu Jia<\/p>\n
%<\/p>\n
% One simple (Gradiant Descent) implementation of the Lasso<\/p>\n
% The objective function is:<\/p>\n
% \u06f1\/\u06f2*sum((Y – X * beta -beta0).^2) + lambda * sum(abs(beta))<\/p>\n
% Comparison with CVX code is done<\/p>\n
% Reference: To be uploaded<\/p>\n
%<\/p>\n
% Version 4.23.2010<\/p>\n
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<\/p>\n
function [beta0,beta] = Lasso(X,Y,lambda)<\/p>\n
n = size(X,1);<\/p>\n
p = size(X,2);<\/p>\n
if(size(Y,1) ~= n)<\/p>\n
fprintf( ‘Error: dim of X and dim of Y are not match!!\\n’)<\/p>\n
quit cancel<\/p>\n
end<\/p>\n
beta0 = 0;<\/p>\n
beta= zeros(p,1);<\/p>\n
crit = 1;<\/p>\n
step = 0;<\/p>\n
while(crit > 1E-5)<\/p>\n
step = step + 1;<\/p>\n
obj_ini = 1\/2*sum((Y – X * beta -beta0).^2) + lambda * sum(abs(beta));<\/p>\n
beta0 = mean(Y – X*beta);<\/p>\n
for(j = 1:p)<\/p>\n
a = sum(X(:,j).^2);<\/p>\n
b = sum((Y – X*beta).* X(:,j)) + beta(j) * a;<\/p>\n
if lambda >= abs(b)<\/p>\n
beta(j) = 0;<\/p>\n
else<\/p>\n
if b>0<\/p>\n
beta(j) = (b – lambda) \/ a;<\/p>\n
else<\/p>\n
beta(j) = (b+lambda) \/a;<\/p>\n
end<\/p>\n
end<\/p>\n
end<\/p>\n
obj_now = 1\/2*sum((Y – X * beta -beta0).^2) + lambda * sum(abs(beta));<\/p>\n
crit = abs(obj_now – obj_ini) \/ abs(obj_ini);<\/p>\n
end<\/p>\n
\n
\n
\u0645\u062b\u0627\u0644 \u062e\u0631\u0648\u062c\u06cc :<\/p>\n
\n
\"example-regress-matlab\"<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"
## Practical session: Linear regression ## Jean-Philippe.Vert@mines.org ## ## In this practical session to test several linear regression methods on the prostate cancer dataset ## #################################### ## Prepare data #################################### # Download prostate data con = url (“http:\/\/www-stat.stanford.edu\/~tibs\/ElemStatLearn\/datasets\/prostate.data”) prost=read.csv(con,row.names=1,sep=”\\t”) # Alternatively, load the file and read from local file as follows # prost=read.csv(‘prostate.data.txt’,row.names=1,sep=”\\t”) # Scale […]<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[19],"tags":[],"class_list":["post-274","post","type-post","status-publish","format-standard","hentry","category-19","category-19-id","post-seq-1","post-parity-odd","meta-position-corners","fix"],"_links":{"self":[{"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=\/wp\/v2\/posts\/274"}],"collection":[{"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=274"}],"version-history":[{"count":4,"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=\/wp\/v2\/posts\/274\/revisions"}],"predecessor-version":[{"id":276,"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=\/wp\/v2\/posts\/274\/revisions\/276"}],"wp:attachment":[{"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=274"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=274"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/vua.nadiran.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=274"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}