library(data.table) drugsens <- fread("data/drug_sens.txt",header = TRUE,drop=1,data.table=F) rownames(drugsens) <- fread("data/drug_sens.txt",header = TRUE,select=1,data.table=F)[,1] drugsens <- as.matrix(drugsens) y <- drugsens[,1] x <- drugsens[,-1] # at first simple linear regression # try small #genes dat <- data.frame(y,x[,1:50]) lmfit <- lm(y~., dat) summary(lmfit) select_var <- step(lmfit) # next elastic lasso library(glmnet) fitmodel <- glmnet(x = x, y = y, family = "gaussian", lambda = 2.5, alpha = 0.01) fitmodel$beta beta <- as.numeric(fitmodel$beta) hist(beta) #parameter tuning alpha <- c(0.01,0.1,0.5,0.8) param <- cbind(alpha,mse=NA) for (i in 1:length(alpha)) { cvfit <- cv.glmnet(x = x, y = y, family = "gaussian", alpha = alpha[i]) param[i,] <- c(alpha = alpha[i],mse = min(cvfit$cvm)) cat(i,"\n") } best_alpha <- param[which.min(param[,2]),1] cvfit_bestalpha <- cv.glmnet(x = x, y = y, family = "gaussian", alpha = best_alpha) best_lambda <- cvfit_bestalpha$lambda.min fitmodel <- glmnet(x = x, y = y, family = "gaussian", lambda = best_lambda, alpha = best_alpha) fitmodel$beta beta <- as.numeric(fitmodel$beta) hist(beta) nonzero_idx <- beta > 0 select_genes <- rownames(fitmodel$beta)[nonzero_idx]