LogisticRegressionForTrialClassification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: Matrix\n",
      "Loading required package: foreach\n",
      "Loaded glmnet 2.0-16\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(boot)\n",
    "library(glmnet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "set.seed(2019)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "train <- read.csv(\"trainC.csv\")\n",
    "test <- read.csv(\"testC.csv\")\n",
    "train <- subset(train, select = -c(sessionDate, trialNum, timeSinceKetamine, animalName))\n",
    "test <- subset(test, select = -c(sessionDate, trialNum, timeSinceKetamine, animalName))\n",
    "\n",
    "#TRAIN 1: ALL COVARIATES PLUS INTERACTION TERMS\n",
    "train1 <- read.csv(\"trainC.csv\")\n",
    "test1 <- read.csv(\"testC.csv\")\n",
    "train1 <- subset(train1, select = -c(sessionDate, trialNum, timeSinceKetamine, animalName))\n",
    "test1 <- subset(test1, select = -c(sessionDate, trialNum, timeSinceKetamine, animalName))\n",
    "\n",
    "#TRAIN 2: ALL COVARIATES NO INTERACTION TERMS\n",
    "train2 <- subset(train1, select = c(totalCellNum,gender,genotype,weight_g,ketamine_day,\n",
    "                                    correlationScore,lickAccuracy,lickNumber,avgFR,\n",
    "                                    avgSingleCellVariance,varianceFR,avgTrialSpeed,\n",
    "                                    varianceSpeed,medianCellDepth,ketBool))\n",
    "test2 <- subset(test1, select = c(totalCellNum,gender,genotype,weight_g,ketamine_day,\n",
    "                                    correlationScore,lickAccuracy,lickNumber,avgFR,\n",
    "                                    avgSingleCellVariance,varianceFR,avgTrialSpeed,\n",
    "                                    varianceSpeed,medianCellDepth,ketBool))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, let's do a 50% split on the training data to determine the best lambda\n",
    "n = length(train[,1])\n",
    "n50 = round(n/2)\n",
    "train50A = train[1:n50,]\n",
    "train50B = train[(n50+1):n,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Logistic Regression Model with Interaction Terms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate test error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================================================================\n",
      "Logistic Regression Model with Interaction Terms\n",
      "\n",
      "Zero-One Loss (10-fold Cross-Validation Average): 0.09182746 \n",
      "Accuracy (10-fold Cross-Validation Average): 0.9081725 \n",
      "=====================================================================\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "n = length(train1[,1])\n",
    "fsize = round(n/k)\n",
    "rmse = rep(0,k)\n",
    "zoloss = rep(0,k)\n",
    "for (i in 1:(k-1)){\n",
    "    # Get train and validation sets\n",
    "    df_train <- train1[-(((i-1)*fsize+1):(i*fsize)),]\n",
    "    df_val <- train1[((i-1)*fsize+1):(i*fsize),]\n",
    "    # Fit model on training and make predictions on validation\n",
    "    model_cv <- glm(ketBool ~ ., data=df_train, family='binomial')\n",
    "    lr_pred_lo <- predict(model_cv,df_val) # lo : log odds\n",
    "    num_val = length(df_val$ketBool)\n",
    "    lr_pred = rep(0,num_val)\n",
    "    actual = rep(0,num_val)\n",
    "    for (j in 1:num_val){\n",
    "        if (lr_pred_lo[j]>0){\n",
    "            lr_pred[j]=1\n",
    "        }\n",
    "    actual[j] = df_val$ketBool[j]\n",
    "    }\n",
    "    # Compute 0-1 loss for each observation\n",
    "    lr_loss = abs(lr_pred-actual) # loss is 0 if NB_pred=actual, 1 otherwise\n",
    "    # Compute mean 0-1 loss on the val set\n",
    "    zoloss[i] = mean(lr_loss)\n",
    "}\n",
    "df_train <- train1[-(((k-1)*fsize+1):n),]\n",
    "df_val <- train1[((k-1)*fsize+1):n,]\n",
    "# Fit model on training and make predictions on validation\n",
    "model_cv <- glm(ketBool ~ ., data=df_train, family='binomial')\n",
    "lr_pred_lo <- predict(model_cv,df_val) # lo : log odds\n",
    "num_val = length(df_val$ketBool)\n",
    "lr_pred = rep(0,num_val)\n",
    "actual = rep(0,num_val)\n",
    "for (j in 1:num_val){\n",
    "    if (lr_pred_lo[j]>0){\n",
    "        lr_pred[j]=1\n",
    "    }\n",
    "    actual[j] = df_val$ketBool[j]\n",
    "}\n",
    "lr_loss = abs(lr_pred-actual)\n",
    "zoloss[k] = mean(lr_loss)\n",
    "test_error_est = mean(zoloss)\n",
    "\n",
    "cat(\"=====================================================================\\n\")\n",
    "cat(\"Logistic Regression Model with Interaction Terms\\n\\n\")\n",
    "cat(\"Zero-One Loss (10-fold Cross-Validation Average):\",test_error_est,\"\\n\")\n",
    "cat(\"Accuracy (10-fold Cross-Validation Average):\",1-test_error_est,\"\\n\")\n",
    "cat(\"=====================================================================\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reduced dataset to match Lasso and Ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================================================================\n",
      "Logistic Regression Model with Interaction Terms\n",
      "\n",
      "Zero-One Loss (10-fold Cross-Validation Average): 0.09505025 \n",
      "Accuracy (10-fold Cross-Validation Average): 0.9049497 \n",
      "=====================================================================\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "n = length(train50B[,1])\n",
    "fsize = round(n/k)\n",
    "rmse = rep(0,k)\n",
    "zoloss = rep(0,k)\n",
    "for (i in 1:(k-1)){\n",
    "    # Get train and validation sets\n",
    "    df_train <- train50B[-(((i-1)*fsize+1):(i*fsize)),]\n",
    "    df_val <- train50B[((i-1)*fsize+1):(i*fsize),]\n",
    "    # Fit model on training and make predictions on validation\n",
    "    model_cv <- glm(ketBool ~ ., data=df_train, family='binomial')\n",
    "    lr_pred_lo <- predict(model_cv,df_val) # lo : log odds\n",
    "    num_val = length(df_val$ketBool)\n",
    "    lr_pred = rep(0,num_val)\n",
    "    actual = rep(0,num_val)\n",
    "    for (j in 1:num_val){\n",
    "        if (lr_pred_lo[j]>0){\n",
    "            lr_pred[j]=1\n",
    "        }\n",
    "    actual[j] = df_val$ketBool[j]\n",
    "    }\n",
    "    # Compute 0-1 loss for each observation\n",
    "    lr_loss = abs(lr_pred-actual) # loss is 0 if NB_pred=actual, 1 otherwise\n",
    "    # Compute mean 0-1 loss on the val set\n",
    "    zoloss[i] = mean(lr_loss)\n",
    "}\n",
    "df_train <- train50B[-(((k-1)*fsize+1):n),]\n",
    "df_val <- train50B[((k-1)*fsize+1):n,]\n",
    "# Fit model on training and make predictions on validation\n",
    "model_cv <- glm(ketBool ~ ., data=df_train, family='binomial')\n",
    "lr_pred_lo <- predict(model_cv,df_val) # lo : log odds\n",
    "num_val = length(df_val$ketBool)\n",
    "lr_pred = rep(0,num_val)\n",
    "actual = rep(0,num_val)\n",
    "for (j in 1:num_val){\n",
    "    if (lr_pred_lo[j]>0){\n",
    "        lr_pred[j]=1\n",
    "    }\n",
    "    actual[j] = df_val$ketBool[j]\n",
    "}\n",
    "lr_loss = abs(lr_pred-actual)\n",
    "zoloss[k] = mean(lr_loss)\n",
    "test_error_est = mean(zoloss)\n",
    "\n",
    "cat(\"=====================================================================\\n\")\n",
    "cat(\"Logistic Regression Model with Interaction Terms\\n\\n\")\n",
    "cat(\"Zero-One Loss (10-fold Cross-Validation Average):\",test_error_est,\"\\n\")\n",
    "cat(\"Accuracy (10-fold Cross-Validation Average):\",1-test_error_est,\"\\n\")\n",
    "cat(\"=====================================================================\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Logistic Regression without Interaction Terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================================================================\n",
      "Logistic Regression Model without Interaction Terms\n",
      "\n",
      "Zero-One Loss (10-fold Cross-Validation Average): 0.1413709 \n",
      "Accuracy (10-fold Cross-Validation Average): 0.8586291 \n",
      "=====================================================================\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "n = length(train2[,1])\n",
    "fsize = round(n/k)\n",
    "rmse = rep(0,k)\n",
    "zoloss = rep(0,k)\n",
    "for (i in 1:(k-1)){\n",
    "    # Get train and validation sets\n",
    "    df_train <- train2[-(((i-1)*fsize+1):(i*fsize)),]\n",
    "    df_val <- train2[((i-1)*fsize+1):(i*fsize),]\n",
    "    # Fit model on training and make predictions on validation\n",
    "    model_cv <- glm(ketBool ~ ., data=df_train, family='binomial')\n",
    "    lr_pred_lo <- predict(model_cv,df_val) # lo : log odds\n",
    "    num_val = length(df_val$ketBool)\n",
    "    lr_pred = rep(0,num_val)\n",
    "    actual = rep(0,num_val)\n",
    "    for (j in 1:num_val){\n",
    "        if (lr_pred_lo[j]>0){\n",
    "            lr_pred[j]=1\n",
    "        }\n",
    "    actual[j] = df_val$ketBool[j]\n",
    "    }\n",
    "    # Compute 0-1 loss for each observation\n",
    "    lr_loss = abs(lr_pred-actual) # loss is 0 if NB_pred=actual, 1 otherwise\n",
    "    # Compute mean 0-1 loss on the val set\n",
    "    zoloss[i] = mean(lr_loss)\n",
    "}\n",
    "df_train <- train2[-(((k-1)*fsize+1):n),]\n",
    "df_val <- train2[((k-1)*fsize+1):n,]\n",
    "# Fit model on training and make predictions on validation\n",
    "model_cv <- glm(ketBool ~ ., data=df_train, family='binomial')\n",
    "lr_pred_lo <- predict(model_cv,df_val) # lo : log odds\n",
    "num_val = length(df_val$ketBool)\n",
    "lr_pred = rep(0,num_val)\n",
    "actual = rep(0,num_val)\n",
    "for (j in 1:num_val){\n",
    "    if (lr_pred_lo[j]>0){\n",
    "        lr_pred[j]=1\n",
    "    }\n",
    "    actual[j] = df_val$ketBool[j]\n",
    "}\n",
    "lr_loss = abs(lr_pred-actual)\n",
    "zoloss[k] = mean(lr_loss)\n",
    "test_error_est = mean(zoloss)\n",
    "\n",
    "cat(\"=====================================================================\\n\")\n",
    "cat(\"Logistic Regression Model without Interaction Terms\\n\\n\")\n",
    "cat(\"Zero-One Loss (10-fold Cross-Validation Average):\",test_error_est,\"\\n\")\n",
    "cat(\"Accuracy (10-fold Cross-Validation Average):\",1-test_error_est,\"\\n\")\n",
    "cat(\"=====================================================================\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GLMNET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, let's do a 50% split on the training data to determine the best lambda\n",
    "n = length(train[,1])\n",
    "n50 = round(n/2)\n",
    "train50A = train[1:n50,]\n",
    "train50B = train[(n50+1):n,]\n",
    "\n",
    "xA = as.matrix(train50A[,-length(train50A)])\n",
    "yA = as.matrix(train50A$ketBool)\n",
    "xB = as.matrix(train50B[,-length(train50B)])\n",
    "yB = as.matrix(train50B$ketBool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0004042946"
     ]
    }
   ],
   "source": [
    "# Select regularization parameter over trainA (50% of training data)\n",
    "model_lasso <- cv.glmnet(xA, yA, family='binomial',alpha=1)\n",
    "lambda_min = model_lasso$lambda.min\n",
    "lambda_1se = model_lasso$lambda.1se\n",
    "cat(lambda_min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================================================================\n",
      "GLMNET Lasso Logistic Regression Model with lambda.min\n",
      "\n",
      "Zero-One Loss (10-fold Cross-Validation Average): 0.09205025 \n",
      "Accuracy (10-fold Cross-Validation Average): 0.9079497 \n",
      "=====================================================================\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "n = length(train50B[,1])\n",
    "fsize = round(n/k)\n",
    "rmse = rep(0,k)\n",
    "zoloss = rep(0,k)\n",
    "for (i in 1:(k-1)){\n",
    "    # Get train and validation sets\n",
    "    xB_train = xB[-(((i-1)*fsize+1):(i*fsize)),]\n",
    "    yB_train = yB[-(((i-1)*fsize+1):(i*fsize)),]\n",
    "    xB_val = xB[((i-1)*fsize+1):(i*fsize),]\n",
    "    yB_val = yB[((i-1)*fsize+1):(i*fsize),]\n",
    "    # Fit model on training and make predictions on validation\n",
    "    model_cv <- glmnet(xB_train, yB_train, family='binomial',alpha=1,lambda=lambda_min)\n",
    "    pred_lo = predict(model_cv, newx = xB_val)\n",
    "    num_val = length(yB_val)\n",
    "    lr_pred = rep(0,num_val)\n",
    "    actual = rep(0,num_val)\n",
    "    for (j in 1:num_val){\n",
    "        if (pred_lo[j]>0){\n",
    "            lr_pred[j]=1\n",
    "        }\n",
    "        actual[j] = yB_val[j]\n",
    "    }\n",
    "    # Compute 0-1 loss for each observation\n",
    "    lr_loss = abs(lr_pred-actual) # loss is 0 if NB_pred=actual, 1 otherwise\n",
    "    # Compute mean 0-1 loss on the val set\n",
    "    zoloss[i] = mean(lr_loss)\n",
    "}\n",
    "xB_train = xB[-(((k-1)*fsize+1):(length(yB))),]\n",
    "yB_train = yB[-(((k-1)*fsize+1):(length(yB))),]\n",
    "xB_val = xB[((k-1)*fsize+1):(length(yB)),]\n",
    "yB_val = yB[((k-1)*fsize+1):(length(yB)),]\n",
    "# Fit model on training and make predictions on validation\n",
    "model_cv <- glmnet(xB_train, yB_train, family='binomial',alpha=1,lambda=lambda_min)\n",
    "pred_lo = predict(model_cv, newx = xB_val)\n",
    "num_val = length(yB_val)\n",
    "lr_pred = rep(0,num_val)\n",
    "actual = rep(0,num_val)\n",
    "for (j in 1:num_val){\n",
    "    if (pred_lo[j]>0){\n",
    "        lr_pred[j]=1\n",
    "    }\n",
    "    actual[j] = yB_val[j]\n",
    "}\n",
    "# Compute 0-1 loss for each observation\n",
    "lr_loss = abs(lr_pred-actual) # loss is 0 if NB_pred=actual, 1 otherwise\n",
    "# Compute mean 0-1 loss on the val set\n",
    "zoloss[k] = mean(lr_loss)\n",
    "test_error_est = mean(zoloss)\n",
    "\n",
    "cat(\"=====================================================================\\n\")\n",
    "cat(\"GLMNET Lasso Logistic Regression Model with lambda.min\\n\\n\")\n",
    "cat(\"Zero-One Loss (10-fold Cross-Validation Average):\",test_error_est,\"\\n\")\n",
    "cat(\"Accuracy (10-fold Cross-Validation Average):\",1-test_error_est,\"\\n\")\n",
    "cat(\"=====================================================================\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:  glmnet(x = xB_train, y = yB_train, family = \"binomial\", alpha = 1,      lambda = lambda_min) \n",
       "\n",
       "     Df   %Dev    Lambda\n",
       "[1,] 55 0.6862 0.0004043"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model_cv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.02988029"
     ]
    }
   ],
   "source": [
    "# Select regularization parameter over trainA (50% of training data)\n",
    "model_lasso <- cv.glmnet(xA, yA, family='binomial',alpha=0)\n",
    "lambda_min = model_lasso$lambda.min\n",
    "lambda_1se = model_lasso$lambda.1se\n",
    "cat(lambda_min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================================================================\n",
      "GLMNET Ridge Logistic Regression Model with lambda.min\n",
      "\n",
      "Zero-One Loss (10-fold Cross-Validation Average): 0.1260879 \n",
      "Accuracy (10-fold Cross-Validation Average): 0.8739121 \n",
      "=====================================================================\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "n = length(train50B[,1])\n",
    "fsize = round(n/k)\n",
    "rmse = rep(0,k)\n",
    "zoloss = rep(0,k)\n",
    "for (i in 1:(k-1)){\n",
    "    # Get train and validation sets\n",
    "    xB_train = xB[-(((i-1)*fsize+1):(i*fsize)),]\n",
    "    yB_train = yB[-(((i-1)*fsize+1):(i*fsize)),]\n",
    "    xB_val = xB[((i-1)*fsize+1):(i*fsize),]\n",
    "    yB_val = yB[((i-1)*fsize+1):(i*fsize),]\n",
    "    # Fit model on training and make predictions on validation\n",
    "    model_cv <- glmnet(xB_train, yB_train, family='binomial',alpha=0,lambda=lambda_min)\n",
    "    pred_lo = predict(model_cv, newx = xB_val)\n",
    "    num_val = length(yB_val)\n",
    "    lr_pred = rep(0,num_val)\n",
    "    actual = rep(0,num_val)\n",
    "    for (j in 1:num_val){\n",
    "        if (pred_lo[j]>0){\n",
    "            lr_pred[j]=1\n",
    "        }\n",
    "        actual[j] = yB_val[j]\n",
    "    }\n",
    "    # Compute 0-1 loss for each observation\n",
    "    lr_loss = abs(lr_pred-actual) # loss is 0 if NB_pred=actual, 1 otherwise\n",
    "    # Compute mean 0-1 loss on the val set\n",
    "    zoloss[i] = mean(lr_loss)\n",
    "}\n",
    "xB_train = xB[-(((k-1)*fsize+1):(length(yB))),]\n",
    "yB_train = yB[-(((k-1)*fsize+1):(length(yB))),]\n",
    "xB_val = xB[((k-1)*fsize+1):(length(yB)),]\n",
    "yB_val = yB[((k-1)*fsize+1):(length(yB)),]\n",
    "# Fit model on training and make predictions on validation\n",
    "model_cv <- glmnet(xB_train, yB_train, family='binomial',alpha=0,lambda=lambda_min)\n",
    "pred_lo = predict(model_cv, newx = xB_val)\n",
    "num_val = length(yB_val)\n",
    "lr_pred = rep(0,num_val)\n",
    "actual = rep(0,num_val)\n",
    "for (j in 1:num_val){\n",
    "    if (pred_lo[j]>0){\n",
    "        lr_pred[j]=1\n",
    "    }\n",
    "    actual[j] = yB_val[j]\n",
    "}\n",
    "# Compute 0-1 loss for each observation\n",
    "lr_loss = abs(lr_pred-actual) # loss is 0 if NB_pred=actual, 1 otherwise\n",
    "# Compute mean 0-1 loss on the val set\n",
    "zoloss[k] = mean(lr_loss)\n",
    "test_error_est = mean(zoloss)\n",
    "\n",
    "cat(\"=====================================================================\\n\")\n",
    "cat(\"GLMNET Ridge Logistic Regression Model with lambda.min\\n\\n\")\n",
    "cat(\"Zero-One Loss (10-fold Cross-Validation Average):\",test_error_est,\"\\n\")\n",
    "cat(\"Accuracy (10-fold Cross-Validation Average):\",1-test_error_est,\"\\n\")\n",
    "cat(\"=====================================================================\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "$lambda\n",
       " [1] 272.25808473 248.07142688 226.03344505 205.95325678 187.65693710\n",
       " [6] 170.98601202 155.79608599 141.95559111 129.34464765 117.85402566\n",
       "[11] 107.38419884  97.84448259  89.15224844  81.23220841  74.01576290\n",
       "[16]  67.44040653  61.44918669  55.99021031  51.01619435  46.48405625\n",
       "[21]  42.35454081  38.59188015  35.16348388  32.03965689  29.19334208\n",
       "[26]  26.59988603  24.23682547  22.08369270  20.12183830  18.33426965\n",
       "[31]  16.70550367  15.22143276  13.86920263  12.63710090  11.51445568\n",
       "[36]  10.49154318   9.55950340   8.71026346   7.93646766   7.23141375\n",
       "[41]   6.58899488   6.00364674   5.47029931   4.98433300   4.54153859\n",
       "[46]   4.13808082   3.77046513   3.43550740   3.13030639   2.85221860\n",
       "[51]   2.59883536   2.36796199   2.15759877   1.96592364   1.79127640\n",
       "[56]   1.63214435   1.48714914   1.35503491   1.23465734   1.12497379\n",
       "[61]   1.02503422   0.93397300   0.85100142   0.77540080   0.70651633\n",
       "[66]   0.64375136   0.58656227   0.53445369   0.48697430   0.44371284\n",
       "[71]   0.40429462   0.36837820   0.33565249   0.30583405   0.27866460\n",
       "[76]   0.25390880   0.23135224   0.21079955   0.19207270   0.17500949\n",
       "[81]   0.15946213   0.14529595   0.13238826   0.12062725   0.10991106\n",
       "[86]   0.10014686   0.09125009   0.08314369   0.07575743   0.06902735\n",
       "[91]   0.06289515   0.05730772   0.05221666   0.04757787   0.04335119\n",
       "[96]   0.03949999   0.03599092   0.03279358   0.02988029\n",
       "\n",
       "$cvm\n",
       " [1] 1.3882764 1.3806430 1.3795975 1.3787386 1.3777988 1.3767709 1.3756469\n",
       " [8] 1.3744183 1.3730759 1.3716096 1.3700090 1.3682624 1.3663577 1.3642818\n",
       "[15] 1.3620207 1.3595598 1.3568834 1.3539752 1.3508181 1.3473940 1.3436879\n",
       "[22] 1.3397679 1.3354508 1.3307902 1.3257697 1.3203703 1.3145736 1.3083619\n",
       "[29] 1.3017189 1.2946299 1.2870820 1.2790649 1.2705711 1.2615964 1.2521402\n",
       "[36] 1.2422059 1.2318012 1.2209384 1.2096341 1.1979099 1.1857918 1.1733096\n",
       "[43] 1.1604976 1.1473933 1.1340373 1.1204723 1.1067429 1.0928944 1.0789725\n",
       "[50] 1.0650224 1.0510882 1.0372123 1.0234352 1.0097943 0.9963245 0.9830571\n",
       "[57] 0.9700203 0.9572386 0.9447331 0.9325213 0.9206174 0.9090322 0.8977736\n",
       "[64] 0.8868464 0.8762529 0.8659929 0.8560640 0.8464618 0.8371803 0.8282083\n",
       "[71] 0.8195339 0.8111630 0.8030765 0.7952626 0.7877091 0.7804033 0.7733321\n",
       "[78] 0.7664818 0.7598389 0.7533895 0.7471200 0.7410197 0.7350696 0.7292573\n",
       "[85] 0.7235693 0.7179961 0.7125239 0.7071407 0.7018280 0.6965918 0.6914066\n",
       "[92] 0.6862732 0.6811821 0.6761314 0.6711201 0.6661358 0.6611965 0.6562907\n",
       "[99] 0.6514272\n",
       "\n",
       "$cvsd\n",
       " [1] 0.0008850502 0.0009310694 0.0008341026 0.0008359960 0.0008385218\n",
       " [6] 0.0008417756 0.0008459140 0.0008511214 0.0008576129 0.0008656379\n",
       "[11] 0.0008754828 0.0008874731 0.0009019749 0.0009193944 0.0009401766\n",
       "[16] 0.0009648019 0.0009937809 0.0010276479 0.0010669529 0.0011122530\n",
       "[21] 0.0011641489 0.0012206087 0.0012870024 0.0013611949 0.0014439253\n",
       "[26] 0.0015356391 0.0016367494 0.0017476351 0.0018686379 0.0020000598\n",
       "[31] 0.0021421590 0.0022951461 0.0024591793 0.0026343595 0.0028207254\n",
       "[36] 0.0030182484 0.0032268285 0.0034462908 0.0036763828 0.0039167743\n",
       "[41] 0.0041670430 0.0044267137 0.0046952174 0.0049719148 0.0052560999\n",
       "[46] 0.0055470046 0.0058438054 0.0061456315 0.0064515738 0.0067606944\n",
       "[51] 0.0070720372 0.0073846387 0.0076975390 0.0080097934 0.0083204826\n",
       "[56] 0.0086287240 0.0089336807 0.0092345710 0.0095306764 0.0098213485\n",
       "[61] 0.0101060159 0.0103841889 0.0106554649 0.0109195321 0.0111761733\n",
       "[66] 0.0114252691 0.0116668015 0.0119008581 0.0121276364 0.0123479384\n",
       "[71] 0.0125625279 0.0127704243 0.0129729128 0.0131708447 0.0133651072\n",
       "[76] 0.0135567133 0.0137467952 0.0139365897 0.0141274160 0.0143206478\n",
       "[81] 0.0145177081 0.0147203087 0.0149286153 0.0151446705 0.0153698034\n",
       "[86] 0.0156047731 0.0158507532 0.0161077823 0.0163776294 0.0166594380\n",
       "[91] 0.0169556979 0.0172639272 0.0175854042 0.0179195045 0.0182650488\n",
       "[96] 0.0186208457 0.0189860081 0.0193599095 0.0197371500\n",
       "\n",
       "$cvup\n",
       " [1] 1.3891614 1.3815741 1.3804316 1.3795746 1.3786373 1.3776127 1.3764928\n",
       " [8] 1.3752694 1.3739335 1.3724753 1.3708845 1.3691499 1.3672597 1.3652012\n",
       "[15] 1.3629609 1.3605246 1.3578772 1.3550029 1.3518850 1.3485062 1.3448521\n",
       "[22] 1.3409885 1.3367378 1.3321514 1.3272136 1.3219060 1.3162103 1.3101095\n",
       "[29] 1.3035876 1.2966299 1.2892241 1.2813600 1.2730303 1.2642307 1.2549609\n",
       "[36] 1.2452241 1.2350280 1.2243846 1.2133105 1.2018266 1.1899588 1.1777363\n",
       "[43] 1.1651928 1.1523652 1.1392934 1.1260193 1.1125867 1.0990400 1.0854240\n",
       "[50] 1.0717830 1.0581602 1.0445970 1.0311327 1.0178041 1.0046449 0.9916858\n",
       "[57] 0.9789540 0.9664732 0.9542638 0.9423426 0.9307234 0.9194164 0.9084290\n",
       "[64] 0.8977659 0.8874291 0.8774182 0.8677308 0.8583627 0.8493080 0.8405562\n",
       "[71] 0.8320964 0.8239334 0.8160494 0.8084334 0.8010742 0.7939600 0.7870789\n",
       "[78] 0.7804184 0.7739663 0.7677102 0.7616377 0.7557401 0.7499983 0.7444020\n",
       "[85] 0.7389391 0.7336009 0.7283747 0.7232485 0.7182057 0.7132512 0.7083623\n",
       "[92] 0.7035372 0.6987675 0.6940509 0.6893852 0.6847566 0.6801825 0.6756506\n",
       "[99] 0.6711644\n",
       "\n",
       "$cvlo\n",
       " [1] 1.3873913 1.3797119 1.3787634 1.3779026 1.3769603 1.3759291 1.3748010\n",
       " [8] 1.3735672 1.3722183 1.3707440 1.3691335 1.3673749 1.3654557 1.3633624\n",
       "[15] 1.3610805 1.3585950 1.3558897 1.3529476 1.3497511 1.3462817 1.3425238\n",
       "[22] 1.3385473 1.3341638 1.3294290 1.3243258 1.3188347 1.3129368 1.3066143\n",
       "[29] 1.2998503 1.2926298 1.2849398 1.2767697 1.2681119 1.2589620 1.2493194\n",
       "[36] 1.2391876 1.2285744 1.2174921 1.2059577 1.1939931 1.1816247 1.1688829\n",
       "[43] 1.1558023 1.1424214 1.1287812 1.1149253 1.1008991 1.0867488 1.0725209\n",
       "[50] 1.0582617 1.0440161 1.0298277 1.0157376 1.0017845 0.9880040 0.9744284\n",
       "[57] 0.9610866 0.9480040 0.9352024 0.9226999 0.9105114 0.8986480 0.8871181\n",
       "[64] 0.8759269 0.8650767 0.8545676 0.8443972 0.8345610 0.8250527 0.8158603\n",
       "[71] 0.8069714 0.7983926 0.7901036 0.7820918 0.7743440 0.7668466 0.7595853\n",
       "[78] 0.7525452 0.7457115 0.7390689 0.7326022 0.7262994 0.7201410 0.7141127\n",
       "[85] 0.7081995 0.7023913 0.6966732 0.6910329 0.6854504 0.6799323 0.6744509\n",
       "[92] 0.6690093 0.6635967 0.6582119 0.6528551 0.6475150 0.6422105 0.6369308\n",
       "[99] 0.6316901\n",
       "\n",
       "$nzero\n",
       " s0  s1  s2  s3  s4  s5  s6  s7  s8  s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 \n",
       " 62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62 \n",
       "s20 s21 s22 s23 s24 s25 s26 s27 s28 s29 s30 s31 s32 s33 s34 s35 s36 s37 s38 s39 \n",
       " 62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62 \n",
       "s40 s41 s42 s43 s44 s45 s46 s47 s48 s49 s50 s51 s52 s53 s54 s55 s56 s57 s58 s59 \n",
       " 62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62 \n",
       "s60 s61 s62 s63 s64 s65 s66 s67 s68 s69 s70 s71 s72 s73 s74 s75 s76 s77 s78 s79 \n",
       " 62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62 \n",
       "s80 s81 s82 s83 s84 s85 s86 s87 s88 s89 s90 s91 s92 s93 s94 s95 s96 s97 s98 \n",
       " 62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62  62 \n",
       "\n",
       "$name\n",
       "           deviance \n",
       "\"Binomial Deviance\" \n",
       "\n",
       "$glmnet.fit\n",
       "\n",
       "Call:  glmnet(x = xA, y = yA, family = \"binomial\", alpha = 0) \n",
       "\n",
       "       Df      %Dev    Lambda\n",
       "  [1,] 62 8.649e-15 272.30000\n",
       "  [2,] 62 5.906e-03 248.10000\n",
       "  [3,] 62 6.476e-03 226.00000\n",
       "  [4,] 62 7.097e-03 206.00000\n",
       "  [5,] 62 7.777e-03 187.70000\n",
       "  [6,] 62 8.520e-03 171.00000\n",
       "  [7,] 62 9.333e-03 155.80000\n",
       "  [8,] 62 1.022e-02 142.00000\n",
       "  [9,] 62 1.119e-02 129.30000\n",
       " [10,] 62 1.225e-02 117.90000\n",
       " [11,] 62 1.341e-02 107.40000\n",
       " [12,] 62 1.467e-02  97.84000\n",
       " [13,] 62 1.605e-02  89.15000\n",
       " [14,] 62 1.755e-02  81.23000\n",
       " [15,] 62 1.919e-02  74.02000\n",
       " [16,] 62 2.097e-02  67.44000\n",
       " [17,] 62 2.291e-02  61.45000\n",
       " [18,] 62 2.501e-02  55.99000\n",
       " [19,] 62 2.729e-02  51.02000\n",
       " [20,] 62 2.977e-02  46.48000\n",
       " [21,] 62 3.246e-02  42.35000\n",
       " [22,] 62 3.529e-02  38.59000\n",
       " [23,] 62 3.842e-02  35.16000\n",
       " [24,] 62 4.179e-02  32.04000\n",
       " [25,] 62 4.542e-02  29.19000\n",
       " [26,] 62 4.933e-02  26.60000\n",
       " [27,] 62 5.353e-02  24.24000\n",
       " [28,] 62 5.803e-02  22.08000\n",
       " [29,] 62 6.284e-02  20.12000\n",
       " [30,] 62 6.798e-02  18.33000\n",
       " [31,] 62 7.345e-02  16.71000\n",
       " [32,] 62 7.926e-02  15.22000\n",
       " [33,] 62 8.542e-02  13.87000\n",
       " [34,] 62 9.192e-02  12.64000\n",
       " [35,] 62 9.878e-02  11.51000\n",
       " [36,] 62 1.060e-01  10.49000\n",
       " [37,] 62 1.135e-01   9.56000\n",
       " [38,] 62 1.214e-01   8.71000\n",
       " [39,] 62 1.296e-01   7.93600\n",
       " [40,] 62 1.382e-01   7.23100\n",
       " [41,] 62 1.470e-01   6.58900\n",
       " [42,] 62 1.560e-01   6.00400\n",
       " [43,] 62 1.654e-01   5.47000\n",
       " [44,] 62 1.749e-01   4.98400\n",
       " [45,] 62 1.846e-01   4.54200\n",
       " [46,] 62 1.945e-01   4.13800\n",
       " [47,] 62 2.045e-01   3.77000\n",
       " [48,] 62 2.146e-01   3.43600\n",
       " [49,] 62 2.248e-01   3.13000\n",
       " [50,] 62 2.350e-01   2.85200\n",
       " [51,] 62 2.452e-01   2.59900\n",
       " [52,] 62 2.554e-01   2.36800\n",
       " [53,] 62 2.655e-01   2.15800\n",
       " [54,] 62 2.755e-01   1.96600\n",
       " [55,] 62 2.854e-01   1.79100\n",
       " [56,] 62 2.951e-01   1.63200\n",
       " [57,] 62 3.047e-01   1.48700\n",
       " [58,] 62 3.142e-01   1.35500\n",
       " [59,] 62 3.234e-01   1.23500\n",
       " [60,] 62 3.325e-01   1.12500\n",
       " [61,] 62 3.413e-01   1.02500\n",
       " [62,] 62 3.499e-01   0.93400\n",
       " [63,] 62 3.583e-01   0.85100\n",
       " [64,] 62 3.665e-01   0.77540\n",
       " [65,] 62 3.744e-01   0.70650\n",
       " [66,] 62 3.822e-01   0.64380\n",
       " [67,] 62 3.896e-01   0.58660\n",
       " [68,] 62 3.969e-01   0.53450\n",
       " [69,] 62 4.039e-01   0.48700\n",
       " [70,] 62 4.108e-01   0.44370\n",
       " [71,] 62 4.174e-01   0.40430\n",
       " [72,] 62 4.238e-01   0.36840\n",
       " [73,] 62 4.301e-01   0.33570\n",
       " [74,] 62 4.361e-01   0.30580\n",
       " [75,] 62 4.420e-01   0.27870\n",
       " [76,] 62 4.477e-01   0.25390\n",
       " [77,] 62 4.532e-01   0.23140\n",
       " [78,] 62 4.586e-01   0.21080\n",
       " [79,] 62 4.638e-01   0.19210\n",
       " [80,] 62 4.690e-01   0.17500\n",
       " [81,] 62 4.740e-01   0.15950\n",
       " [82,] 62 4.789e-01   0.14530\n",
       " [83,] 62 4.836e-01   0.13240\n",
       " [84,] 62 4.883e-01   0.12060\n",
       " [85,] 62 4.930e-01   0.10990\n",
       " [86,] 62 4.975e-01   0.10010\n",
       " [87,] 62 5.020e-01   0.09125\n",
       " [88,] 62 5.064e-01   0.08314\n",
       " [89,] 62 5.107e-01   0.07576\n",
       " [90,] 62 5.150e-01   0.06903\n",
       " [91,] 62 5.193e-01   0.06290\n",
       " [92,] 62 5.235e-01   0.05731\n",
       " [93,] 62 5.277e-01   0.05222\n",
       " [94,] 62 5.319e-01   0.04758\n",
       " [95,] 62 5.360e-01   0.04335\n",
       " [96,] 62 5.401e-01   0.03950\n",
       " [97,] 62 5.442e-01   0.03599\n",
       " [98,] 62 5.483e-01   0.03279\n",
       " [99,] 62 5.523e-01   0.02988\n",
       "[100,] 62 5.562e-01   0.02723\n",
       "\n",
       "$lambda.min\n",
       "[1] 0.02988029\n",
       "\n",
       "$lambda.1se\n",
       "[1] 0.04335119\n",
       "\n",
       "attr(,\"class\")\n",
       "[1] \"cv.glmnet\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model_lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate best model on test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "xTrain = as.matrix(train[,-length(train)])\n",
    "yTrain = as.matrix(train$ketBool)\n",
    "xTest = as.matrix(test[,-length(test)])\n",
    "yTest = as.matrix(test$ketBool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_use = 0.0004043\n",
    "model <- glmnet(xTrain, yTrain, family='binomial',alpha=1,lambda=lambda_use)\n",
    "pred_y = predict(model, newx = xTest)\n",
    "pred_y[pred_y>0]=1\n",
    "pred_y[pred_y<0]=0\n",
    "accuracy = sum(pred_y==yTest)/length(yTest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================================================================\n",
      "Accuracy of Best Model (Lasso with Lambda=0.0004043) on test set: 0.9 \n",
      "=====================================================================\n"
     ]
    }
   ],
   "source": [
    "cat(\"=====================================================================\\n\")\n",
    "cat(\"Accuracy of Best Model (Lasso with Lambda=0.0004043) on test set:\",accuracy,\"\\n\")\n",
    "cat(\"=====================================================================\\n\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}