1

I am trying to do best subset selection on the wine dataset, and then I want to get the test error rate using 10 fold CV. The code I used is - 

cost1 <- function(good, pi=0) mean(abs(good-pi) > 0.5)
res.best.logistic <-
    bestglm(Xy = winedata,
            family = binomial,          # binomial family for logistic
            IC = "AIC",                 # Information criteria
            method = "exhaustive")
res.best.logistic$BestModels
best.cv.err<- cv.glm(winedata,res.best.logistic$BestModel,cost1, K=10)

However, this gives the error - 

Error in UseMethod("family") : no applicable method for 'family' applied to an object of class "NULL"

I thought that $BestModel is the lm-object that represents the best fit, and that's what manual also says. If that's the case, then why cant I find the test error on it using 10 fold CV, with the help of cv.glm?

The dataset used is the white wine dataset from https://archive.ics.uci.edu/ml/datasets/Wine+Quality and the package used is the boot package for cv.glm, and the bestglm package.

The data was processed as - 

winedata <- read.delim("winequality-white.csv", sep = ';')
winedata$quality[winedata$quality< 7] <- "0" #recode
winedata$quality[winedata$quality>=7] <- "1" #recode
winedata$quality <- factor(winedata$quality)# Convert the column to a factor
names(winedata)[names(winedata) == "quality"] <- "good"      #rename 'quality' to 'good'
user1993
  • 498
  • 1
  • 10
  • 22
  • Could you edit your post specifying which is the package you are using and if the dataset is in that package? Otherwise your post is not reproducible. – Scipione Sarlo Apr 04 '20 at 10:41
  • @ScipioneSarlo, I have added the dataset and packages used – user1993 Apr 04 '20 at 11:23
  • please, do not share the url of your data but try to respect the reproducible suggestions of SO for a good question: https://stackoverflow.com/help/minimal-reproducible-example – Scipione Sarlo Apr 04 '20 at 11:30
  • @ScipioneSarlo, actually this wine dataset is a very well known dataset for studying different machine learning methods, and many textbooks regularly reference it. Its a public dataset. I have also added the code I used to pre-process the data, for the sake of completeness! – user1993 Apr 04 '20 at 11:38
  • How did you create the `cost1` input into `cv.glm()`? – Len Greski Apr 04 '20 at 12:25
  • @LenGreski, I have added my `cost1` function now; basically if something is >0.5, its considered `good` – user1993 Apr 04 '20 at 12:37

1 Answers1

2

bestglm fit rearranges your data and name your response variable as y, hence if you pass it back into cv.glm, winedata does not have a column y and everything crashes after that

It's always good to check what is the class:

class(res.best.logistic$BestModel)
[1] "glm" "lm"

But if you look at the call of res.best.logistic$BestModel:

res.best.logistic$BestModel$call

glm(formula = y ~ ., family = family, data = Xi, weights = weights)

head(res.best.logistic$BestModel$model)
  y fixed.acidity volatile.acidity citric.acid residual.sugar chlorides
1 0           7.0             0.27        0.36           20.7     0.045
2 0           6.3             0.30        0.34            1.6     0.049
3 0           8.1             0.28        0.40            6.9     0.050
4 0           7.2             0.23        0.32            8.5     0.058
5 0           7.2             0.23        0.32            8.5     0.058
6 0           8.1             0.28        0.40            6.9     0.050
  free.sulfur.dioxide density   pH sulphates
1                  45  1.0010 3.00      0.45
2                  14  0.9940 3.30      0.49
3                  30  0.9951 3.26      0.44
4                  47  0.9956 3.19      0.40
5                  47  0.9956 3.19      0.40
6                  30  0.9951 3.26      0.44

You can substitute things in the call etc, but it's too much of a mess. Fitting is not costly, so make a fit on winedata and pass it to cv.glm:

best_var = apply(res.best.logistic$BestModels[,-ncol(winedata)],1,which)
# take the variable names for best model
best_var = names(best_var[[1]])
new_form = as.formula(paste("good ~", paste(best_var,collapse="+")))
fit = glm(new_form,winedata,family="binomial")

best.cv.err<- cv.glm(winedata,fit,cost1, K=10)
StupidWolf
  • 45,075
  • 17
  • 40
  • 72
  • Thanks a lot!! That makes a lot of sense. Is it possible to use this same method for calculating the 10 fold CV for a ridge regression model? I also posted a question about it - https://stackoverflow.com/questions/61029739/calculating-cross-validation-error-for-ridge-regression-model – user1993 Apr 04 '20 at 14:44
  • the result from glmnet is not of class "glm". you need to use cv.glmnet – StupidWolf Apr 04 '20 at 14:47