Adrian
Adrian

Reputation: 9803

How to predict with cv.glm in R

Suppose I have the dataset:

mydata=data.frame(
  status = as.factor(c(0,0,0,0,0,1,1,1,1,1)),
  a = c(1,3,4,5,6,1,2,3,4,5),
  b = c(4,2,3,6,2,1,3,4,5,6)
)

I construct a glm model, train it on half of the observations, and use that model to predict the response of the other half of the observations.

   train=sample(1:nrow(mydata), nrow(mydata)/2)
   test=mydata[-train,]
   test.response=status[-train]
   fit=glm(status~., data=mydata, family="binomial", subset=train)
   probs=predict(fit, test, type="response")
   pred=rep(0,5)
   pred[probs>0.5]=1
   table(pred, test.response)

And this tells me how many true positives and true negatives I have predicted (I had 2, and 2, respectively).

     test.response
pred 0 1
   0 2 0
   1 1 2

Instead of manually coding the train and test datasets, I turne to cv.glm so that R can just cross-validate for me.

   library(boot)
   fit2=glm(status~., data=mydata, family="binomial")
   cv.fit=cv.glm(mydata, fit2, K=2)

My question is, how can I use the cross validated model to predict my response variable? What I used before was probs=predict(fit, test, type="response"), but in this case I don't know what test is.

Upvotes: 3

Views: 9135

Answers (1)

Hack-R
Hack-R

Reputation: 23210

I believe this is what you want:

require(boot)
require(glmnet)
status <- c(0,0,0,0,0,1,1,1,1,1)
ymat   <- as.matrix(status)
xdata  <- data.frame(a,b)
xmat   <- as.matrix(xdata)
fit.cv <- cv.glmnet(y = ymat, x = xmat, family="binomial")
lmin   <- fit.cv$lambda.min
l1se   <- fit.cv$lambda.1se
net    <- glmnet(y = ymat, x = xmat, family="binomial")
predict(net, s=lmin, type = "nonzero")
colnames(status)[predict(net, s=lmin,type="nonzero")$X1]
plot(cv.glmnet(y = ymat,x = xmat, family= "binomial", type = "class"))
plot(cv.glmnet(y = ymat,x = xmat, family= "binomial", type = "mae"))
plot(cv.glmnet(y = ymat,x = xmat, family= "binomial", type = "deviance"))
plot(cv.glmnet(y = ymat,x = xmat, family= "binomial", type = "mse"))
plot(cv.glmnet(y = ymat,x = xmat, family= "binomial", type = "auc")) #needs more data

The dataset was really tiny so I recommend trying this on a bigger dataset to avoid warnings/errors.

Upvotes: 2

Related Questions