Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I was able to reproduce table 3.1 from ESL. However, when I tried to reproduce table 3.2, my estimated coefficients were way off (shown below):

         [,1]
 [1,]  0.4292
 [2,]  0.5765
 [3,]  0.6140
 [4,] -0.0190
 [5,]  0.1448
 [6,]  0.7372
 [7,] -0.2063
 [8,] -0.0295
 [9,]  0.0095

The results from table 3.2 in ESL are as follow: enter image description here

I have attached my code here:

res <- read.table("prostate.data", sep = "")
fix(res)
XTraining = subset(res, train)
XTesting = subset(res, train == FALSE)


nrow = dim( XTraining )[1]
p = dim( XTraining )[2] - 1 # the last column is the response 

D = XTraining[,1:p-1] # get the predictor data

# This gives the raw data for Table 3.1 from the book:
# 
print(cor(D),digts=3)

library(xtable)
xtable( cor(D), caption="Duplication of the values from Table 3.1 from the book", digits=3 )

#
# Duplicate Table 3.2 from the book 
#

# Append a column of ones:
# 
Dp = cbind( matrix(1,nrow,1), as.matrix( D ) )

lpsa = XTraining[,p]

library(MASS)

betaHat = ginv( t(Dp) %*% Dp ) %*% t(Dp) %*% as.matrix(lpsa)

# this is basically the first column in Table 3.2:
#
print('first column: beta estimates')
print(betaHat,digits=2)

Does anyone have similar issues ? I would like to know what went wrong with my output.

share|improve this question

1 Answer 1

up vote 5 down vote accepted

The problem is you have not scaled the data as described in the information file:

http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/prostate.info.txt

The features must first be scaled to have mean zero and variance 96 (=n) before the analyses in Tables 3.1 and beyond. That is, if x is the 96 by 8 matrix of features, we compute xp <- scale(x,TRUE,TRUE)

This gives the correct answer:

X <- read.table("prostate.data", sep = "")
train <- X[,10]
y <- X[,9]
X <- X[,-(9:10)]

X.scaled <- scale(X,TRUE,TRUE)

betaHat <- lm(y[train]~X.scaled[train,])$coef
print(betaHat)

         (Intercept)  X.scaled[train, ]lcavol X.scaled[train, ]lweight     X.scaled[train, ]age    X.scaled[train, ]lbph 
          2.46493292               0.67952814               0.26305307              -0.14146483               0.21014656 
X.scaled[train, ]svi     X.scaled[train, ]lcp X.scaled[train, ]gleason   X.scaled[train, ]pgg45 
          0.30520060              -0.28849277              -0.02130504               0.26695576 
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.