I have created this post to spark a discussion around an idea I had regarding the choice of the seed number used to replicate the results of a statistical model. Here is some background on how I came up with this idea and how it may be applied.
I have a data set, called data, which has 506 observations, 13 predictor variables, and 1 continuous response variable. I split the data into a training and test set by randomly sampling 380 indexes without replacement, storing the observations in the rows equal to the 380 indexes from the data dataframe into a dataframe called train, and storing the remaining 126 observations into a dataframe called test.
Once I created the training and test sets, I built a linear model using the glm function, predicted the outcomes of the test set using the predict function, and computed the MSE between the test set's predicted values and its actual values. In this instance, I obtained an MSE of 16.65786. I left the code here for the night as it was time to go home.
The R code for this is:
library(MASS)
data <- Boston
index <- sample(x = 1:nrow(data), size = round(nrow(data)*0.75), replace = FALSE)
train <- data[index, ]
test <- data[-index, ]
lm1.fit <- glm(medv ~ ., data = train)
pr.lm1 <- predict(lm1.fit, test)
(MSE.lm1 <- sum((pr.lm1 - test$medv)^2) / nrow(test))
The following day, I realized that I should have set a seed so that the random sample of indexes (rows) used to create the training and test sets are the same each time the code is run. So I chose a seed of 1 and reran the code above. This time I received an MSE of 28.72933, a huge difference from the MSE I obtained the night before. This sparked the question of what if the seed I chose was one that lead to an outlier case or that the unseeded result that I obtained was an outlier case. Maybe both results were outlier cases! So I did some investigating.
I used a for loop 10,000 times, each time getting a new random sample of indices, created a new training and test set, built a linear model, predicted outcomes for the test set, and calculated the MSE. I stored each of the 10,000 MSEs in a dataframe and computed the mean and median MSE. I also looked at the distribution of the MSEs which appeared to be approximately normal. I reran the forloop again just to make sure that the mean and median of the MSEs did not change much, they didn't.
The first 10,000 MSEs had a mean of 24.13 and a median of 23.51.
The second 10,000 MSEs had a mean of 24.01 and a median of 23.48.
The R code for this is:
error <- data.frame(err = rep(0, 10000))
for(i in 1:10000)
{
index <- sample(x = 1:nrow(data), size = round(nrow(data)*0.75),
replace = FALSE)
train <- data[index, ]
test <- data[-index, ]
#Linear Model 1 - No Data Cleaning
lm1.fit <- glm(medv ~ ., data = train)
summary(lm1.fit)
pr.lm1 <- predict(lm1.fit, test)
error[i, 1] <- sum((pr.lm1 - test$medv)^2) / nrow(test)
}
summary(error[, 1])
The last thing I did was loop through this same model building process until I found a seed that lead to an MSE in between the mean and the median of the 20,000 MSEs.
The R code for this is:
for(i in 1:10000)
{
set.seed(i)
index <- sample(x = 1:nrow(data), size = round(nrow(data)*0.75),
replace = FALSE)
train <- data[index, ]
test <- data[-index, ]
#Linear Model 1 - No Data Cleaning
lm1.fit <- glm(medv ~ ., data = train)
summary(lm1.fit)
pr.lm1 <- predict(lm1.fit, test)
error <- sum((pr.lm1 - test$medv)^2) / nrow(test)
if(error > 23 & error < 24){print(i)}
}
I chose to use a seed of 47 which lead to an MSE of 23.63085. I believe that this method for obtaining a seed more accurately depicts a model's ability to accurately predict versus just randomly picking a seed or finding a seed which gives the best results (which I view as "cheating").
My question to all you statisticians out there is does this seem like a reasonable ideology / approach? Is there a more appropriate method out there? Is there a downside to picking a seed this way?
Thank you everyone in advance for your input!
[Scortchi] Thanks - I was trying to be pithy & conflated two points: (1) Specifically the OP's re-inventing Monte-Carlo cross-validation; & (2) the sensitivity of the out-of-sample MSE estimate to the seed used shows that they should be using some form of re-sampling validation rather than a train/test split. – Scortchi
[Brandon] @Scortchi So what you're saying is I should simply use the following code instead:
library(MASS)
data <- Boston
lm1.fit <- glm(medv ~ ., data = data)
library(boot)
set.seed(1)
cvResult <- cv.glm(data = data, glmfit = lm1.fit, K = 10)
cvResult$delta
Which leads to a cv prediction error of 23.38074, very close to the prediction error obtained using my long method above. Am I correct?
