ロジスティック回帰分析による喫煙者の肺がん罹患予測

目次
1. [はじめに](#1-はじめに)
2. [分析概要](#2-分析概要)
3. [分析前の準備](#3-分析前の準備)
4. [分析](#4-分析)
5. [再分析](#5-再分析)
6. [コードまとめ](#6-コードまとめ)
7. [おわりに](#7-おわりに)
## 1. はじめに

機械の商社で営業事務として働いています。

昨夏から所属している課でサポートすることになった業務の一つに、主に工場ラインなどに使用される制御システムのソフトウェアライセンスの販売があり、DXへの興味関心が高まりました。
DX化によって業務効率化を図る顧客が数多いる中で、私の勤め先では社内運営やビジネス展開にデータ分析があまり活用されていないことにも気がつきました。

せっかく興味を持ったDXの世界について、まず身の回りでできる学習としてPythonとデータ分析を学びたいと思い、自分のペースで学習を進められるAidemy premiumのデータ分析講座の受講を決意しました。

そんな私の現在の座席は喫煙者に囲まれています。
右隣、そのまた右隣、真後ろ、右斜め前、左斜め前はみな愛煙家で、タバコ休憩の度にさまざまな情報を仕入れてくる優秀な営業マンたちなのですが、やはり気になるのは彼らの健康です。

今回、Aidemy premiumの受講修了のための課題を公開するということで、愛煙家営業マンたちに「数年以内に禁煙しよう！」と思ってもらえるような分析に取り組みたいと思い、喫煙者が肺がんになる確率を予測するモデルを作ってみることにしました。

## 2. 分析概要

今回使用するデータセットと環境は、kaggleで公開されている「Smoking related lung cancers」です。

[Smoking related lung cancers](https://www.kaggle.com/datasets/raddar/smoking-related-lung-cancers)

ファイル名は「lung_cancer.csv」です。


US National Lung Screening Trial (NLST)により行われた調査で、喫煙者・喫煙経験者を対象に肺がんを発症したかどうかを7年間追跡し、性別・人種・調査開始時の年齢・肺がん発症時のステージ・調査開始から何年後に発症したかという項目のデータが記載されています。

本分析では、性別、調査開始時の年齢、肺がん発症時のステージのデータに着目し、肺がん罹患予測を行います。


## 3. 分析前の準備
3.1データの読み込み

```python 
#使用するライブラリのインポート
import numpy as np
import pandas as pd
from pandas import DataFrame

#データの読み込み
train_df = pd.read_csv('/kaggle/input/smoking-related-lung-cancers/lung_cancer.csv')
#データ内容の確認
train_df.info()
```
出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/e5dfc569-690f-319d-02e1-ebb1a0647ea2.png)

→データフレーム型で、53,427のデータ数があることを確認できました。

3.2 各項目の特徴

3.2.1 肺がんステージ
まず、分析時にターゲットとする予定のstage_of_cancerの特徴を確認します。

```python
train_df['stage_of_cancer'].unique()
```
出力結果
array([nan, 'IA', 'IIB', 'IV', 'IIIB', 'IIIA', 'IIA', 'IB'], dtype=object)

stage_of_cancerには欠損値があり、欠損値の部分は肺がんにかからなかった人と見なせるため、欠損値を"0"とし、ステージIを1、ステージIIとIIIを2、ステージIVを3に分けて、0-3の人数を可視化します。

```python
import matplotlib.pyplot as plt
import seaborn as sns

#'stage_of_cancer'の欠損値にゼロを入れる
train_df = train_df.fillna({'stage_of_cancer':0,})

#'stage_of_cancer'でがんのステージ別にIを"1"、IIとIIIを"2"、IVを"3"と定義する。
train_df['stage_of_cancer'] = train_df['stage_of_cancer'].replace({'IA': 1, 'IB' : 1,'IIA':2,
                                                       'IIB':2,'IIIA':2,'IIIB':2,'IV':3}).infer_objects(copy=False)

#'stage_of_cancer'の中でがんを発症した人の人数と発症していない人の数をステージ別に確認
train_df.groupby('stage_of_cancer')['stage_of_cancer'].count()
train_df.info()
count = train_df.groupby('stage_of_cancer')['stage_of_cancer'].count()
print(count)

# 'stage_of_cancer'を文字型に変換
from sklearn.preprocessing import LabelEncoder
label_encoder = LabelEncoder()
train_df['stage_of_cancer'] = label_encoder.fit_transform(train_df['stage_of_cancer'].astype(str))

# stage_of_cancerごとにヒストグラムを作成
plt.figure(figsize=(10, 6))
train_df['stage_of_cancer'].value_counts().plot(kind='bar')
plt.xlabel('Stage of Cancer')
plt.ylabel('Number of Patients')
plt.title('Number of Patients by Stage of Cancer')
plt.xticks(rotation=0)
plt.show()
```
各ステージごとの患者数の棒グラフ
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/9774b2e8-03ff-fe7d-61d2-879976831fee.png)


各ステージごとの患者数
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/81df0697-ed65-e090-b767-891df144bdbf.png)


→実際に肺がんを発症した人数はかなり少ないことがわかります。

3.2.2 性別
'gender'の特徴を確認します。

```python
# 'gender'ごとに患者数をカウントして円グラフを作成
labels = ['Male', 'Female']
sizes = [59, 41]
plt.figure(figsize=(3, 3))
plt.pie(sizes, labels=labels, autopct='%1.1f%%', startangle=140)
plt.axis('equal') 
plt.show()
```

性別の割合
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/42f3465f-d120-0d3c-b56f-c8856e09e107.png)


→男性の方がやや多いことがわかります。


3.2.3 喫煙者・喫煙経験者
'smoker'の特徴を確認します。

```python
# 'smoker'の特徴の確認
print(train_df['smoker'].unique())

# 'smoker'ごとに患者数をカウントして円グラフを作成
labels = ['Current', 'Former']
sizes = [52, 48]
plt.figure(figsize=(3, 3))
plt.pie(sizes, labels=labels, autopct='%1.1f%%', startangle=140)
plt.axis('equal') 
plt.show()
```
喫煙者・喫煙経験者の割合
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/04595d6d-5f8a-1c35-6968-673181785451.png)


→喫煙者・喫煙経験者の割合に差はありません。


3.2.4 年齢
最後に'age'の特徴を確認します。

```python
# 'age'の特徴の確認
print(train_df['age'].unique())

# 'age'ごとに患者数をカウントしてヒストグラムを作成
plt.figure(figsize=(10, 6))
sns.countplot(data=train_df, x='age', palette='viridis')
plt.xlabel('Age')
plt.ylabel('Count')
plt.title('Count of Age')
plt.show()
```
年齢別人数
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/4a0fc1cf-9b1e-bc05-fbdf-b1d78c7d4bf9.png)
				
						

					
3.3 性別、喫煙者・喫煙経験者、年齢における肺がんステージの割合を可視化				

					
```python
#'gender','smoker','age'別に、各肺がんステージの割合を可視化
col = ['age', 'gender', 'smoker']
fig = plt.figure(figsize=(15,15))

for i in range(len(col)):
    plt.subplot(3,2,i+1)
    plt.title(col[i])
    sns.histplot(data=train_df,y=train_df[col[i]],hue='stage_of_cancer')
plt.tight_layout()
plt.show()
```
出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/95e30f44-6f08-fd99-bdfa-c2d14647926e.png)

→肺がんが発症していない人の割合がほとんどを占めている中で、喫煙経験者より現在の喫煙者の肺がん罹患割合がわずかに多いことはわかります。


3.4 データの準備
分析前に、カテゴリカル変数である'gender','smoker','age','stage_of_cancer'をOne-Hot-Encodingによりモデルで取り扱えるように変換します。

```python
# 'gender','smoker','age','stage_of_cancer'をOne-Hot-Encodingで変換
train_df1 = pd.get_dummies(train_df[['gender','age','smoker','stage_of_cancer']])
# データを確認
train_df1.head()
```

出力内容
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/27b0fe16-9973-1d80-95e5-1e1b6a3eee1e.png)



## 4. 分析 
4.1 多項ロジスティック回帰モデルによる分析
ターゲットとする肺がんステージが0-3の4段階のため、多項ロジスティック回帰による分析を行います。

```python
from sklearn.model_selection import train_test_split
target = train_df1['stage_of_cancer']

train = train_df1.drop(columns=['stage_of_cancer']) 

#訓練用と検証用に分ける
X_train,X_val,y_train,y_val = train_test_split(train, target, test_size=0.2, shuffle=True, random_state=0)

from sklearn.linear_model import LogisticRegression

#モデルを定義し学習
multi_logreg = LogisticRegression(multi_class='multinomial', solver='lbfgs')
multi_logreg.fit(X_train, y_train)

#正答率を算出
print('Training Accuracy: ', multi_logreg.score(X_train, y_train))
print('Testing Accuracy: ', multi_logreg.score(X_val, y_val))
```

出力結果
Training Accuracy:  0.9619802999461875
Testing Accuracy:  0.9619127830806663

→高い正答率を出していますが、比較のため他のモデルによる分析も行います。

4.2 Random Forestによる分析

```python
from sklearn.ensemble import RandomForestClassifier
​
clf = RandomForestClassifier(random_state=0)
#モデルを定義し学習
clf.fit(X_train,y_train)
y_pred = clf.predict(X_val)
​
#正答率を算出
print('Training Accuracy: ', clf.score(X_train, y_train))
print('Testing Accuracy: ', clf.score(X_val, y_val))
```
出力結果
Training Accuracy:  0.9619802999461875
Testing Accuracy:  0.9619127830806663


4.3 サポートベクターマシーン(SVM)による分析

```python
from sklearn.svm import SVC

# SVMモデルのインスタンス化
svm_model = SVC(kernel='linear', C=1.0, random_state=42)
# モデルを定義し学習
svm_model.fit(X_train, y_train)

#正答率を算出
print('Training Accuracy: ', clf.score(X_train, y_train))
print('Testing Accuracy: ', clf.score(X_val, y_val))
```

出力結果
Training Accuracy:  0.9619802999461875
Testing Accuracy:  0.9619127830806663


⇒いずれのモデルでも高く同程度の正答率となっていますが、これは分析しているデータのほとんどが肺がんに罹患していない("0")ことによると考えられます。そのため、下記条件に変更の上再度分析してみることにします。
・肺がんに罹患していない人("0")の人数を減らす。
・がんステージI-IVを"1"にまとめ、"0"or"1"の二値分類を行う。


## 5. 再分析

5.1 変更データの読み込みから可視化

```python
import numpy as np
import pandas as pd
from pandas import DataFrame

#データの読み込み
train_df = pd.read_csv('/kaggle/input/lung-cancer-2-1/lung_cancer_2.csv')
#データの表示
train_df

#データ内容の確認
train_df.info()
```

出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/94a334f2-bd0a-5cc3-afd4-07d9ee17a60b.png)

→肺がんに罹患していない人を50,913人減らし、2514人分のデータに変更しました。

次に、がんステージの分類について、I-IVを全て"1"と定義します。

```python
train_df['stage_of_cancer'] = train_df['stage_of_cancer'].replace({'IA': 1, 'IB' : 1,'IIA':1,
                                                       'IIB':1,'IIIA':1,'IIIB':1,'IV':1}).infer_objects(copy=False)
```

あとは最初の分析と同様に各内容を可視化します。

肺がん発症・未発症患者の人数の棒グラフ
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/1569c270-f47d-a2f9-b15e-ff703ce2f127.png)


肺がん発症・未発症患者の人数
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/1e535a3c-20c8-24a8-b4b2-7126bbce522f.png)

'gender','smoker','age'別に、肺がん発症・未発症の割合を可視化
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3842015/a798bc52-e987-7059-84fc-77f36a35cb1b.png)

→当初データと比較すると、女性より男性、喫煙経験者より現在喫煙者、調査時60歳台の発症の割合が高いことがわかります。


5.2 ロジスティック回帰による学習と正答率の算出
ロジスティック回帰の二値分類を実行します。

```python
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score

# 反復回数を1,000回に設定し、モデルを定義し学習
model = LogisticRegression(max_iter = 1000)
model.fit(X_train, y_train)

y_pred = model.predict(X_train)

#正答率を算出
print('Training Accuracy: ', model.score(X_train, y_train))
print('Testing Accuracy: ', model.score(X_val, y_val))
```
出力結果
Training Accuracy:  0.8055693684733963
Testing Accuracy:  0.8190854870775348

⇒当初分析より正答率が下がりましたが、極端に偏ったデータの分析ではなくなりました。




## 6. コードまとめ

<details><summary>当初コード</summary>

```python
import numpy as np
import pandas as pd
from pandas import DataFrame
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.model_selection import train_test_split

#データの読み込み
train_df = pd.read_csv('/kaggle/input/smoking-related-lung-cancers/lung_cancer.csv')
#データの表示
train_df
#データ内容の確認
train_df.info()

#'stage_of_cancer'の特徴の確認
train_df['stage_of_cancer'].unique()
#'stage_of_cancer'の欠損値にゼロを入れる
train_df = train_df.fillna({'stage_of_cancer':0})

#'stage_of_cancer'でがんのステージ別にIを"1"、IIとIIIを"2"、IVを"3"と定義する。
train_df['stage_of_cancer'] = train_df['stage_of_cancer'].replace({'IA': 1, 'IB' : 1,'IIA':2,
                                                       'IIB':2,'IIIA':2,'IIIB':2,'IV':3}).infer_objects(copy=False)
                                                       
#'stage_of_cancer'の中でがんと診断された人の人数と診断されていない人の数を確認
train_df.groupby('stage_of_cancer')['stage_of_cancer'].count()
#データ内容の確認
train_df.info()

#'stage_of_cancer'の中でがんと診断された人の人数と診断されていない人の数を確認
count = train_df.groupby('stage_of_cancer')['stage_of_cancer'].count()
print(count)

# 'stage_of_cancer'を文字型に変換
from sklearn.preprocessing import LabelEncoder
label_encoder = LabelEncoder()
train_df['stage_of_cancer'] = label_encoder.fit_transform(train_df['stage_of_cancer'].astype(str))

# stage_of_cancerごとにヒストグラムを作成
plt.figure(figsize=(10, 6))
train_df['stage_of_cancer'].value_counts().plot(kind='bar')
plt.xlabel('Stage of Cancer')
plt.ylabel('Number of Patients')
plt.title('Number of Patients by Stage of Cancer')
plt.xticks(rotation=0)
plt.show()

# 'gender'の特徴の確認
print(train_df['gender'].unique())

# 'gender'ごとに患者数をカウントして円グラフを作成
labels = ['Male', 'Female']
sizes = [59, 41]
plt.figure(figsize=(3, 3))
plt.pie(sizes, labels=labels, autopct='%1.1f%%', startangle=140)
plt.axis('equal') 
plt.show()

# 'smoker'の特徴の確認
print(train_df['smoker'].unique())

# 'smoker'ごとに患者数をカウントして円グラフを作成
labels = ['Current', 'Former']
sizes = [52, 48]
plt.figure(figsize=(3, 3))
plt.pie(sizes, labels=labels, autopct='%1.1f%%', startangle=140)
plt.axis('equal') 
plt.show()

# 'age'の特徴の確認
print(train_df['age'].unique())

# 'age'ごとに患者数をカウントしてヒストグラムを作成
plt.figure(figsize=(10, 6))
sns.countplot(data=train_df, x='age', palette='viridis')
plt.xlabel('Age')
plt.ylabel('Count')
plt.title('Count of Age')
plt.show()

#'age','gender','smoker'別に、各肺がんステージの割合を可視化
col = ['age', 'gender', 'smoker']
fig = plt.figure(figsize=(15,15))

for i in range(len(col)):
    plt.subplot(3,2,i+1)
    plt.title(col[i])
    sns.histplot(data=train_df,y=train_df[col[i]],hue='stage_of_cancer')
plt.tight_layout()
plt.show()

# 'gender''age','gender','smoker','stage_of_cancer'をOne-Hot-Encodingで変換
train_df1 = pd.get_dummies(train_df[['gender','age','smoker','stage_of_cancer']])
# データを確認
train_df1.head()

from sklearn.model_selection import train_test_split
target = train_df1['stage_of_cancer']

train = train_df1.drop(columns=['stage_of_cancer'])  

#訓練用と検証用に分ける
X_train,X_val,y_train,y_val = train_test_split(train, target, test_size=0.2, shuffle=True, random_state=0)

#多項ロジスティック回帰モデルを定義し学習
from sklearn.linear_model import LogisticRegression
multi_logreg = LogisticRegression(multi_class='multinomial', solver='lbfgs')
multi_logreg.fit(X_train, y_train)

print('Training Accuracy: ', multi_logreg.score(X_train, y_train))
print('Testing Accuracy: ', multi_logreg.score(X_val, y_val))

from sklearn.ensemble import RandomForestClassifier

clf = RandomForestClassifier(random_state=0)
#ランダムフォレストモデルを定義し学習
clf.fit(X_train,y_train)
y_pred = clf.predict(X_val)

print('Training Accuracy: ', clf.score(X_train, y_train))
print('Testing Accuracy: ', clf.score(X_val, y_val))

from sklearn.svm import SVC

svm_model = SVC(kernel='linear', C=1.0, random_state=42)
#SVM回帰モデルを定義し学習
svm_model.fit(X_train, y_train)

print('Training Accuracy: ', clf.score(X_train, y_train))
print('Testing Accuracy: ', clf.score(X_val, y_val))
```
</details>


<details><summary>データ変更後コード</summary>

```python
import numpy as np
import pandas as pd
from pandas import DataFrame
import matplotlib.pyplot as plt
import seaborn as sns

#データの読み込み
train_df = pd.read_csv('/kaggle/input/lung-cancer-2-1/lung_cancer_2.csv')
#データの表示
train_df
#データ内容の確認
train_df.info()


#'stage_of_cancer'の特徴の確認
train_df['stage_of_cancer'].unique()

#'stage_of_cancer'と'days_to_cancer'の欠損値にゼロを入れる
train_df = train_df.fillna({'stage_of_cancer':0,'days_to_cancer':0})

#がんのステージI-IVを全て"1"と定義
train_df['stage_of_cancer'] = train_df['stage_of_cancer'].replace({'IA': 1, 'IB' : 1,'IIA':1,
                                                       'IIB':1,'IIIA':1,'IIIB':1,'IV':1}).infer_objects(copy=False)

#'stage_of_cancer'の中でがんと診断された人の人数と診断されていない人の数を確認
train_df.groupby('stage_of_cancer')['stage_of_cancer'].count()
train_df.info()

#'stage_of_cancer'の中でがんと診断された人の人数と診断されていない人の数を確認
count = train_df.groupby('stage_of_cancer')['stage1_of_cancer'].count()
print(count)

#'stage_of_canecer'を文字型に変換
from sklearn.preprocessing import LabelEncoder
label_encoder = LabelEncoder()
train_df['stage_of_cancer'] = label_encoder.fit_transform(train_df['stage_of_cancer'].astype(str))
#データ内容の確認
train_df.info()

# stage_of_cancerごとにヒストグラムを作成
plt.figure(figsize=(10, 6))
train_df['stage_of_cancer'].value_counts().plot(kind='bar')
plt.xlabel('Stage of Cancer')
plt.ylabel('Number of Patients')
plt.title('Number of Patients by Stage of Cancer')
plt.xticks(rotation=0)
plt.show()

# 'stage_of_cancer'を文字型に変換
label_encoder = LabelEncoder()
train_df['stage_of_cancer'] = label_encoder.fit_transform(train_df['stage_of_cancer'].astype(str))

# データ内容の確認
train_df.info()

#'age','gender','smoker'別に、各肺がんステージの割合を可視化
col = ['age', 'gender', 'smoker']
fig = plt.figure(figsize=(15,15))

for i in range(len(col)):
    plt.subplot(3,2,i+1)
    plt.title(col[i])
    sns.histplot(data=train_df,y=train_df[col[i]],hue='stage_of_cancer')
plt.tight_layout()
plt.show()

# 'gender''age','gender','smoker','stage_of_cancer'をOne-Hot-Encodingで変換
train_df1 = pd.get_dummies(train_df[['gender','age','smoker','stage_of_cancer']])
# データを確認
train_df1.head()

from sklearn.model_selection import train_test_split
target = train_df1['stage_of_cancer']

train = train_df1.drop(columns=['stage_of_cancer'])  

#訓練用と検証用に分ける
X_train,X_val,y_train,y_val = train_test_split(train, target, test_size=0.2, shuffle=True, random_state=0)

#ロジスティック回帰二値分類モデルを定義し学習
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score

model = LogisticRegression(max_iter = 1000)
model.fit(X_train, y_train)

y_pred = model.predict(X_train)

print('Training Accuracy: ', model.score(X_train, y_train))
print('Testing Accuracy: ', model.score(X_val, y_val))

```
</details>



## 7. おわりに


想定より喫煙者の肺がん罹患割合は高くなかったのですが、それでも禁煙した人の罹患割合は低かったので、愛煙家の営業マンたちに少しは注意喚起できる分析ができたかなと思います。

約3か月間の教材による学習と、slackや技術カウンセリングによる先生方のご指導により、複数のモデルを使用した分析の方法を身につけることができました。

統計学の知識を増やすことも重要だと認識しているので、コード記載の練習と平行して統計学の学習も行いたいと思います。

会社方針にて勤め先のDX化は近い未来ではありませんが、DX化の際には主担当に抜擢されるよう学習を継続していきたいと思います。