{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rを利用した細胞診 (FNA) 結果からの乳がん診断\n",
    "\n",
    "## 概要\n",
    "\n",
    "このノートブックでは、細胞診結果にもとづく乳がん診断をRを使って行います。細胞の画像から直接診断するのではなく、細胞の均一さなどの特徴から診断します。このような特徴量と診断結果のデータセットが、[Breast Cancer Wisconsin (Diagnostic) Data Set](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic))として公開されているので、これを利用します。\n",
    "\n",
    "\n",
    "## データのダウンロードと確認\n",
    "\n",
    "まずデータをダウンロードします。ダウンロードする `breast-cancer-wisconsin.data`はヘッダのない csv 形式のファイルです。`read.csv`すると、適当に設定されたヘッダ名とともにデータを見ることができます。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url <- \"https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data\"\n",
    "data <- read.csv(url)\n",
    "head(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "このままだとデータがわかりにくいので正しいヘッダを付与します。\n",
    "- id: サンプルのID\n",
    "- CT: 細胞の塊の厚さ\n",
    "- UCSize: 細胞サイズの均一さ\n",
    "- UCShape: 細胞の形の均一さ\n",
    "- MA: 辺縁接着\n",
    "- SECS: 単一上皮細胞サイズ\n",
    "- BN: 裸核\n",
    "- BC: 無刺激性クロマチン\n",
    "- NN: 通常の核\n",
    "- M: 分裂\n",
    "- diagnosis: 良性なら2、悪性なら4\n",
    "\n",
    "この問題は良性の2と悪性の4を認識する問題ですが、2種類に分ける際、通常は0か1に分けることが多いです。そこで、 outcome という列を作って、diagnosis が 4 なら outcome を 1 に、diagnosis が 2 なら outcome を 0 にするようにします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data <- read.csv(file = url, header = FALSE,\n",
    "                col.names = c(\"id\",\"CT\", \"UCSize\", \"UCShape\", \"MA\", \"SECS\", \"BN\", \"BC\", \"NN\",\"M\", \"diagnosis\") )\n",
    "\n",
    "data$outcome[data$diagnosis==4] = 1\n",
    "data$outcome[data$diagnosis==2] = 0\n",
    "data$outcome = as.integer(data$outcome)\n",
    "head(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## データの前処理と分析\n",
    "\n",
    "データのうち id は予測と無関係な値であることが自明ですので除去しましょう。また、BNには欠損値が含まれていますので、今回は削除します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(dplyr)\n",
    "library(tidyverse)\n",
    "data2 <- data %>% select(-id, -BN)\n",
    "\n",
    "data2$outcome[data2$diagnosis==4] = 1\n",
    "data2$outcome[data2$diagnosis==2] = 0\n",
    "data2$outcome = as.integer(data2$outcome)\n",
    "head(data2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PerformanceAnalytics のライブラリをインストールすることでデータ間の相関を分析することができます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "install.packages(\"PerformanceAnalytics\")\n",
    "library(PerformanceAnalytics)\n",
    "chart.Correlation(data2[, c(1:10)], histogram=TRUE, col=\"grey10\", pch=1, main=\"Cancer Means\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "データセットで良性と悪性の分布をグラフ化します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(data2, aes(x = outcome)) +\n",
    "  geom_bar(aes(fill = \"blue\")) +\n",
    "  ggtitle(\"Distribution of diagnosis for the entire dataset\") +\n",
    "  theme(legend.position=\"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "データの8割を学習データにして、残りをテストデータにします。学習データのみを学習に利用して、テストデータでモデルを評価します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sample_size = floor(0.8 * nrow(data2))\n",
    "\n",
    "set.seed(0)\n",
    "train_set = sample(seq_len(nrow(data2)), size = sample_size)\n",
    "\n",
    "training = data2[train_set, ]\n",
    "\n",
    "testing = data2[-train_set, ]\n",
    "\n",
    "head(training)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## モデルの学習\n",
    "\n",
    "良性か悪性かの2つに分類する問題のため、単純な方式であるロジスティック回帰を利用します。ロジスティック回帰とは、以下のような回帰式です。\n",
    "\n",
    "$P=\\frac{e^y}{1+e^y}$  \n",
    "$y = w_0 + w_1 x_1 + w_2 x_2 + \\cdots w_n x_n$\n",
    "\n",
    "ここで CT や UCSize の値が $x_1$や $x_2$になります。ただし、yを outcomeそのものとして線形回帰を行うと、0や1以外の値を取りうる可能性が出てくるので、シグモイド関数を取り入れて 0から1の間になるようにします。シグモイド関数を計算して得られる P は、例えば、良性や悪性の確率とみなすことができます。\n",
    "\n",
    "R でロジスティック回帰を行う際は glm を利用します。これは一般化線形モデルのライブラリで、2項分布 binomial を仮定すればロジスティック回帰になります。どの説明変数を使うか、どれを目的変数とするかは、数式のような形で指定します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_algorithm = model = glm(outcome ~ CT + \n",
    "                                UCSize +\n",
    "                                UCShape +\n",
    "                                MA +\n",
    "                                SECS + \n",
    "                                BC  +\n",
    "                                NN  +\n",
    "                                M ,\n",
    "                               family=binomial(link='logit'), control = list(maxit = 50),data=training)\n",
    "\n",
    "print(summary(model_algorithm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "学習が終わったらテストしてみましょう。悪性の確率が0.5を超えるものは悪性とみなします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prediction_testing = predict(model_algorithm,testing, type = \"response\")\n",
    "prediction_testing = ifelse(prediction_testing > 0.5, 1, 0)\n",
    "error = mean(prediction_testing != testing$outcome)\n",
    "print(paste('Model Accuracy',1-error))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上記では0.5を悪性か良性かの判断に使いましたが、0.5よりも小さくすると多くを悪性と見なす傾向になり、悪性の見落としは減りますが、良性も悪性とみなすため診断効率は下がります。良い診断というのは、どのようなしきい値でも良い精度が得られるというものです。\n",
    "\n",
    "これを評価するために、しきい値を変えながら悪性の正しい判定（True Positive）と良性の誤検出 (False Positive) を計算し、カーブをかいたものがROC曲線とよばれているものです。カーブの右上はしきい値が最も低い状況で、悪性も全て検出しますが（True Positive =1.0)、良性も悪性と検出します（False Positive =1.0)。\n",
    "\n",
    "このROC曲線の内側の領域 (AUC: Area under curve) が広いほうが良い評価となります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "install.packages(\"ROCR\")\n",
    "library(ROCR)\n",
    "p = predict(model_algorithm, testing, type=\"response\")\n",
    "pr = prediction(p, testing$outcome)\n",
    "prf = performance(pr, measure = \"tpr\", x.measure = \"fpr\")\n",
    "plot(prf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SageMaker を利用して学習と推論を行う\n",
    "\n",
    "より大規模かつ複雑なデータの場合に、ロジスティック回帰よりも高度なアルゴリズムを使いたくなるかもしれません。SageMakerでは、表形式のデータに対して高い分類性能を誇るXGBoostを利用できます。詳細は以下の例も参考にしてください。\n",
    "\n",
    "https://github.com/aws/amazon-sagemaker-examples/blob/master/r_examples/r_batch_transform/r_xgboost_batch_transform.ipynb\n",
    "\n",
    "まずは必要なライブラリをインストールします。SageMaker のPythonパッケージを使うために、reticulate をインストールしてからSageMakerをインポートします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "install.packages('tictoc', repos='http://cran.us.r-project.org')\n",
    "library(tictoc)\n",
    "# Install reticulate library and import sagemaker\n",
    "library(reticulate)\n",
    "sagemaker <- import('sagemaker')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SageMaker に必要なセッションの情報や保存先 (bucket)などを指定します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "session <- sagemaker$Session()\n",
    "bucket <- session$default_bucket()\n",
    "prefix <- 'r-batch-transform'\n",
    "role_arn <- sagemaker$get_execution_role()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "さきほどロジスティック回帰で利用したデータをここでも利用します。**注意点**として、SageMakerのXGBoostを利用する際は、**一列目にラベルが必要**ですので列の順番を入れ替えます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_data <- data2\n",
    "xgb_data <- xgb_data[, c(10, 1, 2, 3, 4, 5, 6, 7, 8)]\n",
    "head(xgb_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "学習データとバリデーションデータを用意します。用意したら SageMaker の学習用に CSV ファイルとして保存します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = floor(0.8 * nrow(xgb_data))\n",
    "\n",
    "set.seed(0)\n",
    "train_valid_set = sample(seq_len(nrow(xgb_data)), size = sample_size)\n",
    "train_valid_data = xgb_data[train_valid_set, ]\n",
    "testing = xgb_data[-train_valid_set, ]\n",
    "\n",
    "sample_size = floor(0.8 * nrow(train_valid_data))\n",
    "train_set = sample(seq_len(nrow(train_valid_data)), size = sample_size)\n",
    "training = train_valid_data[train_set, ]\n",
    "validation = train_valid_data[-train_set, ]\n",
    "\n",
    "write_csv(training, 'train.csv', col_names = FALSE)\n",
    "write_csv(validation, 'valid.csv', col_names = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SageMaker の upload_data の関数でそれぞれのファイルを S3 にアップロードします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s3_train <- session$upload_data(path = 'train.csv', \n",
    "                                bucket = bucket, \n",
    "                                key_prefix = 'data')\n",
    "s3_valid <- session$upload_data(path = 'valid.csv', \n",
    "                                bucket = bucket, \n",
    "                                key_prefix = 'data')\n",
    "s3_train_input <- sagemaker$inputs$TrainingInput(s3_data = s3_train,\n",
    "                                     content_type = 'csv')\n",
    "s3_valid_input <- sagemaker$inputs$TrainingInput(s3_data = s3_valid,\n",
    "                                     content_type = 'csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SageMaker のXGBoost のアルゴリズムは、AWS が管理するコンテナイメージで提供されます。そのイメージの場所を知っておく必要があるため、get_image_uriで取得します。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "container <- sagemaker$amazon$amazon_estimator$get_image_uri(session$boto_region_name, 'xgboost',repo_version='1.0-1')\n",
    "container"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "取得が完了すれば、そのコンテナイメージに加えて、学習用のインスタンスタイプや数を指定して学習を行います。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s3_output <- paste0('s3://', bucket, '/output')\n",
    "estimator <- sagemaker$estimator$Estimator(image_uri = container,\n",
    "                                           role = role_arn,\n",
    "                                           train_instance_count = 1L,\n",
    "                                           train_instance_type = 'ml.m5.large')\n",
    "\n",
    "\n",
    "tic(\"Model Fitting\")\n",
    "estimator$set_hyperparameters(num_round = 100L,objective = 'binary:logistic')\n",
    "job_name <- paste('sagemaker-train-xgboost', format(Sys.time(), '%H-%M-%S'), sep = '-')\n",
    "input_data <- list('train' = s3_train_input,\n",
    "                   'validation' = s3_valid_input)\n",
    "estimator$fit(inputs = input_data, job_name = job_name)\n",
    "toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "学習が終わるとモデルをホストして、推論を行う環境を構築することができます。CSV形式でデータを送信するので、CSVSerializer を指定します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tic(\"Model Deployment\")\n",
    "model_endpoint <- estimator$deploy(initial_instance_count = 1L,\n",
    "                                   instance_type = 'ml.t2.medium')\n",
    "model_endpoint$serializer <- sagemaker$predictor$CSVSerializer()\n",
    "toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "テストデータを取り出してエンドポイントにリクエストを送って結果を取得します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_data<- testing[-1]\n",
    "num_predict_rows <- nrow(testing)\n",
    "test_sample <- as.matrix(test_data[1:num_predict_rows, ])\n",
    "dimnames(test_sample)[[2]] <- NULL\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tic(\"Invoke Endpoint\")\n",
    "library(stringr)\n",
    "predictions <- model_endpoint$predict(test_sample)\n",
    "predictions <- str_split(predictions, pattern = ',', simplify = TRUE)\n",
    "predictions <- as.numeric(predictions)\n",
    "predictions = ifelse(predictions > 0.5, 1, 0)\n",
    "toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "取得した結果を表示したり、精度・ROCを確認します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert predictions to Integer\n",
    "result <- cbind(predicted_outcome = as.integer(predictions), \n",
    "                      testing[1:num_predict_rows, ])\n",
    "head(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error = mean(result$predicted_outcome != result$outcome)\n",
    "print(paste('Model Accuracy',1-error))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pr = prediction(result$predicted_outcome, result$outcome)\n",
    "prf = performance(pr, measure = \"tpr\", x.measure = \"fpr\")\n",
    "plot(prf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最後にエンドポイントを削除します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "session$delete_endpoint(model_endpoint$endpoint)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}