見出し画像

【Python】メンタルヘルスに関するアンケート結果から、精神疾患の治療必要性を予測する

メンタルヘルスに関するアンケート結果から、精神疾患の治療必要性を予測するモデルをPythonで構築しました。

開発環境

  • PC:ASUS ExpertBook

  • ブラウザ:Google Chrome

  • Pythonのバージョン:3.7.15

  • Google Colaboratory

使用したデータ

今回はkaggleのデータセット「Mental Health in Tech Survey」を使用しました。
このデータは2014 年に行われた調査から得られたもので、労働者のメンタルヘルスに対する態度と精神疾患の治療経験、職場のメンタルヘルス対策などを測定しています。

kaggle上では、データセットに含まれるデータ(変数)について以下のような説明がありました。今回はtreatmentを目的変数として、treatmentをその他の変数で予測するモデルを構築していきます。

  • Timestamp:タイムスタンプ

  • Age:年齢

  • Gender:性別

  • Country:国

  • state: 米国に住んでいる場合、どの州または準州に住んでいるか

  • self_employed : 自営業かどうか

  • family_history : 家族が精神疾患に罹患した経験があるか

  • treatment:精神疾患の治療を求めたことがありますか?

  • work_interfere : メンタルヘルスに問題がある場合、それが仕事の妨げになっていると感じるか

  • no_employees : 会社または組織には何人の従業員がいるか

  • remote_work : 少なくとも 50% の時間、リモート (オフィスの外) で仕事をしているか

  • tech_company : あなたの雇用主は主にテクノロジー企業/組織か

  • benefits:雇用主はメンタルヘルスの福利厚生を提供しているか

  • care_options : 雇用主が提供するメンタルヘルスケアのオプションを知っているか

  • wellness_program : 雇用主は、従業員の健康プログラムの一環としてメンタルヘルスについて話したことがあるか

  • seek_help : 雇用主は、メンタルヘルスの問題と助けを求める方法についてさらに学ぶためのリソースを提供しているか

  • anonymity:メンタルヘルスや薬物乱用の治療リソースを利用することを選択した場合、匿名性は保護されるか

  • leave:精神疾患のための医療休暇はどの程度取りやすいか

  • mentalhealthconsequence:メンタルヘルスの問題について雇用主と話し合うと、マイナスの結果が生じるか

  • physhealthconsequence:身体の健康問題について雇用主と話し合うことは、マイナスの結果をもたらすか

  • coworkers :メンタルヘルスの問題について同僚と話し合う気はあるか

  • supervisor:直属の上司とメンタルヘルスの問題について話し合うことはできるか

  • mentalhealthinterview:面接で潜在的な雇用主にメンタルヘルスの問題を持ち出すか

  • physhealthinterview:面接で潜在的な雇用主に身体的な健康問題を持ち出すか

  • mentalvsphysical:あなたの雇用主はメンタルヘルスをフィジカルヘルスと同じくらい真剣に考えていると思るか

  • obs_consequence : 職場でメンタルヘルスの問題を抱えた同僚が悪影響を受けていることを聞いたり、観察したりしたことがあるか

  • comments:自由コメント

1.データの中身の確認

全体像を大まかに把握

実際に今回のデータがどのようなものか、以下のコードで確認していきます。

import pandas as pd

# データセットの読み込み
survey_data = pd.read_csv("./survey.csv")

# displayで全ての列を表示させる
pd.set_option("display.max_columns", None)

# データセットの上5行を表示
display(survey_data.head())
  • データセットの読み込み:kaggleからダウンロードしたcsvを、Colaboratory左側のメニューからアップロードした上で、上記のコードを記述しました(図1)。

図1:Colaboratory上でのデータセットのアップロード
  • displayで全ての列を表示させる:データセットの読み込み後、すぐにdisplayを行ったところ、列数が多すぎて途中の列が省略(…)されるという現象が起きました(図2)。ネットで調べたところ、pd.set_option("display.max_columns", None)を使用することで、そのコード以降のdisplayで、全ての列が表示されることがわかりました。それをdisplay前に記述した上で、displayを用いました。

図2:pd.set_optionをしない状態でのデータセット
  • データセットの上5行を表示:途中までのデータを表したものが図3です。ざっと見たところ、Ageのみが数値データ、それ以外はカテゴリカルデータ(選択肢のどれを選んだかを表すデータ。Yes-No=?など、四則演算できない)であることがわかります。

図3:全列を表示したデータセット(実際には右側にデータが続いて表示されている)

要約統計量の確認

上記ではデータセットの上5行だけを表示しましたが、より全体像を把握するため、各変数の要約統計量を見ていきたいと思います。

# 数値データの統計量を表示
display(survey_data.describe())

# カテゴリカルデータの統計量を表示
display(survey_data.describe(exclude="number"))

# Genderの選択肢を確認
print(df["Gender"].value_counts())
  • 数値データの統計量を表示:先述の通り、Ageのみの統計量が表示されました。

  • カテゴリカルデータの統計量を表示

    • アウトプット(図4)を確認したところ、count行に「1259」という数字が多く並んでいます。つまりこれがデータの総数(行数、つまり回答者数)ということがわかりました。

    • 一方、state・self_employed・work_interfereなどは、1259を下回っています。つまり、これらのデータには欠損値が含まれるということです。

    • 図4でもう一つ気になるのは、uniqueの行です。Countryやstateは妥当な数ですが、Genderの数がやけに多いことがわかります。選択肢式であればもう少し少なくても良いはずです。

図4:カテゴリカルデータの要約統計量
  • Genderの回答値を確認

    • Genderのunique数が多かったことに疑問を感じたため、実際の回答値とその数を確認することにしました。結果の一部が図5です。「Male」「male」「M」「Man」のような表記ゆれや、性的マイノリティを示す表現が複数あることがわかります。つまりGenderは選択肢での回答ではなく、自由記述での回答だったのです。

    • このことから、Genderをモデルに組み込む場合、表記ゆれの修正や、性的マイノリティのグルーピング(まとめて「その他」とするなど)などのデータクレンジング・ハンドリングを行う必要が出てきました。この時点で、作業工数の多さを考え、分析からGenderを除外することとしました。

    • 余談ですが、このような問題を避けるため、アンケートの場合、性別は「男」「女」「その他」などとするのが通例です。データクレンジングで疲弊しないためにも、どのようにデータを収集するか、アンケートの設計時点で考えることが重要です。

図5:Genderの回答値とその数(さらに下に続いている)

2.データクレンジング(学習前のデータ処理)

続いてデータクレンジングとして、以下のコードを記述しました。

# 分析に使用しないTimestamp・comments・Gender列、偏りの大きいCountry列、
# 欠損の多いstate列を削除/リストワイズ削除
survey_data = survey_data.drop(["Timestamp", "comments", "Gender", "Country", 
                                "state"], axis=1).dropna()

# Ageをカテゴリカル変数化し、Ageを削除
survey_data["AgeBand"] = pd.qcut(survey_data["Age"], 4)
survey_data = survey_data.drop(["Age"], axis=1)

# One-hot encodingでtreatment以外のカテゴリカル変数をダミー変数化
survey_data = pd.get_dummies(survey_data, columns= ["self_employed", 
              "family_history",	"work_interfere", "no_employees",	
              "remote_work", "tech_company", "benefits", "care_options", 
              "wellness_program", "seek_help", "anonymity", "leave", 
              "mental_health_consequence", "phys_health_consequence",	
              "coworkers", "supervisor", "mental_health_interview", 
              "phys_health_interview", "mental_vs_physical",
              "obs_consequence", "AgeBand"])

# treatmentは、Yesを1、Noを0に置換
survey_data = survey_data.replace({"treatment": {"Yes": 1, "No": 0}})

# 一旦統計量を表示
display(survey_data.describe())
  • 分析に使用しないTimestamp・comments・Gender列、偏りの大きいCountry列、欠損の多いstate列を削除/リストワイズ削除

    • Countryは、データの半数以上がUnited Statesであることから削除(分析から除外)しました。

    • stateは、CountryがUnited Statesの人のみが回答しているため欠損が多くこれも削除しました。

    • 上記変数を削除した上で、リストワイズ削除(欠損値を含む行を全て削除)を実行しています。

  • Ageをカテゴリカル変数化し、Ageを削除:Ageのみが数値データ、かつ連続数です。外れ値が存在する場合、学習に悪影響が出る恐れがあります。そこで新変数AgeBandを作成し、四分位でAgeを4分割し、カテゴリカル変数化しました。その上で、元のAgeを削除しています。

  • One-hot encodingでtreatment以外のカテゴリカル変数をダミー変数化

    • 1でも見た通り、今回のデータセットは「Yes」「No」などカテゴリカルかつ文字の変数がほとんどです。ただし文字情報のままでは学習ができません。

    • そこで、One-hot encodingを施すことで、各変数をダミー変数化しました。例えばself_employedでYesと回答した人は、「self_employed_Yes」という変数で1になり、Yesと回答しなかった人(=Noと回答した人)は0になる、という変換が可能です。

    • この変換を、目的変数であるtreatment以外に施しました。

  • treatmentは、Yesを1、Noを0に置換:上記One-hot encodingは施せませんが、treatmentも「Yes」「No」の文字情報です。これもこのままでは学習に使用できないため、1と0に変換しました。

  • 一旦統計量を表示:上記の作業を終えた上で、再び要約統計量を表示したものが図6です。最終的な回答者数は977となりました。

図6:データクレンジング後の要約統計量

3.モデルの学習と最適なモデルの探索

最後に、以下のコードでモデルの学習と予測精度の向上を行いました。

import scipy.stats
from sklearn.linear_model import LogisticRegression
from sklearn.svm import LinearSVC
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import train_test_split
from sklearn.metrics import f1_score

# 学習データと評価データに分割
train_X, test_X, train_y, test_y = train_test_split(
    survey_data.drop("treatment", axis=1),
    survey_data["treatment"], test_size=0.3, random_state=42)

# グリッドサーチ用にモデルとパラメーターセットをまとめた辞書を用意
model_param_set_grid = {
    LogisticRegression(): {
        "C": [10 ** i for i in range(-5, 5)],
        "random_state": [42]
    },
    LinearSVC(): {
        "C": [10 ** i for i in range(-5, 5)],
        "multi_class": ["ovr", "crammer_singer"],
        "random_state": [42]
    },
    SVC(): {
        "kernel": ["linear", "poly", "rbf", "sigmoid"],
        "C": [10 ** i for i in range(-5, 5)],
        "decision_function_shape": ["ovr", "ovo"],
        "random_state": [42]
    },
    DecisionTreeClassifier(): {
        "max_depth": [i for i in range(1, 20)],
    },
    RandomForestClassifier(): {
        "n_estimators": [i for i in range(10, 20)],
        "max_depth": [i for i in range(1, 10)],
    },
    KNeighborsClassifier(): {
        "n_neighbors": [i for i in range(1, 10)]
    }
}

# ランダムサーチ用にモデルとパラメーターセットをまとめた辞書を用意
model_param_set_random = {
    LogisticRegression(): {
        "C": scipy.stats.uniform(0.00001, 1000),
        "random_state": scipy.stats.randint(0, 100)
    },
    LinearSVC(): {
        "C": scipy.stats.uniform(0.00001, 1000),
        "multi_class": ["ovr", "crammer_singer"],
        "random_state": scipy.stats.randint(0, 100)
    },
    SVC(): {
        "kernel": ["linear", "poly", "rbf", "sigmoid"],
        "C": scipy.stats.uniform(0.00001, 1000),
        "decision_function_shape": ["ovr", "ovo"],
        "random_state": scipy.stats.randint(0, 100)
    },
    DecisionTreeClassifier(): {
        "max_depth": scipy.stats.randint(1, 20),
    },
    RandomForestClassifier(): {
        "n_estimators": scipy.stats.randint(10, 100),
        "max_depth": scipy.stats.randint(1, 20),
    },
    KNeighborsClassifier(): {
        "n_neighbors": scipy.stats.randint(1, 20)
    }
}

# スコア比較用に変数を用意
max_score = 0
best_model = None
best_param = None

# グリッドサーチでパラメーターサーチ
for model, param in model_param_set_grid.items():
    clf = GridSearchCV(model, param)
    clf.fit(train_X, train_y)
    pred_y = clf.predict(test_X)
    score = f1_score(test_y, pred_y, average="micro")
    # 最高評価更新時にモデルやパラメーターも更新
    if max_score < score:
        max_score = score
        best_model = model.__class__.__name__
        best_param = clf.best_params_

# ランダムサーチでパラメーターサーチ
for model, param in model_param_set_random.items():
    clf = RandomizedSearchCV(model, param)
    clf.fit(train_X, train_y)
    pred_y = clf.predict(test_X)
    score = f1_score(test_y, pred_y, average="micro")
    # 最高評価更新時にモデルやパラメーターも更新
    if max_score < score:
        max_score = score
        best_model = model.__class__.__name__
        best_param = clf.best_params_

# 最適な学習モデル・パラメーター・ベストスコア        
print("学習モデル:{},\nパラメーター:{}".format(best_model, best_param))
print("ベストスコア:",max_score)
  • 学習データと評価データに分割:機械学習では、データセット全体のデータを使用するのではなく、モデル構築のための「学習データ」と、構築されたモデルを評価するための「評価データ」に分けて使用します。一般的な割合は、学習データ7:評価データ3です。

  • 「グリッドサーチ用にモデルとパラメーターセットをまとめた辞書を用意」以降

    • 機械学習では、学習の目的やデータの形式などによって、用いるべき手法が異なります。さらに、各手法において精度を高めるためには、「ハイパーパラメーター」と呼ばれる要素を人為的に設定しなければなりません。

    • 上記のコードでは、「グリッドサーチ」「ランダムサーチ」と呼ばれる方法を用いて、最適な手法とハイパーパラメーターを自動的に出力できる機能を実装しています。

    • グリッドサーチでは、ハイパーパラメーターの値の候補を複数指定し、ハイパーパラメーターセットを作成します。モデルの評価を繰り返すことにより、最適なハイパーパラメーターセットが作成されます。

    • 一方ランダムサーチでは、ハイパーパラメーターが取りうる値の範囲を指定し、確率で決定されたハイパーパラメーターセットを用いてモデルの評価を行うことを繰り返し、最適なハイパーパラメーターセットを探します。

    • 結果として、今回最適な学習モデルは線形SVM、ハイパーパラメーターはCが0.01, multi_classがcrammer_singer、random_stateが42となり、F値は0.8163265306122449となりました。モデルの性能を表すF値としては、かなり高いのではないでしょうか。

終わりに

コードまとめ

今回記述したコードをまとめると、以下のようになります。なお、私の環境では、コードの実行に8分9秒かかりました。分割して実行したところ、特に3の部分で時間がかかっているようです。

import pandas as pd
import scipy.stats
from sklearn.linear_model import LogisticRegression
from sklearn.svm import LinearSVC
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import train_test_split
from sklearn.metrics import f1_score



# 1.データの中身の確認

# データセットの読み込み
survey_data = pd.read_csv("./survey.csv")

# displayで全ての列を表示させる
pd.set_option("display.max_columns", None)

# データセットの上5行を表示
display(survey_data.head())

# 数値データの統計量を表示
display(survey_data.describe())

# カテゴリカルデータの統計量を表示
display(survey_data.describe(exclude="number"))

# Genderの選択肢を確認
print(df["Gender"].value_counts())



# 2.データクレンジング

# 分析に使用しないTimestamp・comments・Gender列、偏りの大きいCountry、
# 欠損の多いstate列を削除/リストワイズ削除
survey_data = survey_data.drop(["Timestamp", "comments", "Gender", "Country", 
                                "state"], axis=1).dropna()

# Ageをカテゴリカル変数化し、Ageを削除
survey_data["AgeBand"] = pd.qcut(survey_data["Age"], 4)
survey_data = survey_data.drop(["Age"], axis=1)

# One-hot encodingでtreatment以外のカテゴリカル変数をダミー変数化
survey_data = pd.get_dummies(survey_data, columns= ["self_employed", 
              "family_history",	"work_interfere", "no_employees",	
              "remote_work", "tech_company", "benefits", "care_options", 
              "wellness_program", "seek_help", "anonymity", "leave", 
              "mental_health_consequence", "phys_health_consequence",	
              "coworkers", "supervisor", "mental_health_interview", 
              "phys_health_interview", "mental_vs_physical",
              "obs_consequence", "AgeBand"])

# treatmentは、Yesを1、Noを0に置換
survey_data = survey_data.replace({"treatment": {"Yes": 1, "No": 0}})

# 一旦統計量を表示
display(survey_data.describe())



# 3.モデルの学習と最適なモデルの探索

# 学習データと評価データに分割
train_X, test_X, train_y, test_y = train_test_split(
    survey_data.drop("treatment", axis=1),
    survey_data["treatment"], test_size=0.3, random_state=42)

# グリッドサーチ用にモデルとパラメーターセットをまとめた辞書を用意
model_param_set_grid = {
    LogisticRegression(): {
        "C": [10 ** i for i in range(-5, 5)],
        "random_state": [42]
    },
    LinearSVC(): {
        "C": [10 ** i for i in range(-5, 5)],
        "multi_class": ["ovr", "crammer_singer"],
        "random_state": [42]
    },
    SVC(): {
        "kernel": ["linear", "poly", "rbf", "sigmoid"],
        "C": [10 ** i for i in range(-5, 5)],
        "decision_function_shape": ["ovr", "ovo"],
        "random_state": [42]
    },
    DecisionTreeClassifier(): {
        "max_depth": [i for i in range(1, 20)],
    },
    RandomForestClassifier(): {
        "n_estimators": [i for i in range(10, 20)],
        "max_depth": [i for i in range(1, 10)],
    },
    KNeighborsClassifier(): {
        "n_neighbors": [i for i in range(1, 10)]
    }
}

# ランダムサーチ用にモデルとパラメーターセットをまとめた辞書を用意
model_param_set_random = {
    LogisticRegression(): {
        "C": scipy.stats.uniform(0.00001, 1000),
        "random_state": scipy.stats.randint(0, 100)
    },
    LinearSVC(): {
        "C": scipy.stats.uniform(0.00001, 1000),
        "multi_class": ["ovr", "crammer_singer"],
        "random_state": scipy.stats.randint(0, 100)
    },
    SVC(): {
        "kernel": ["linear", "poly", "rbf", "sigmoid"],
        "C": scipy.stats.uniform(0.00001, 1000),
        "decision_function_shape": ["ovr", "ovo"],
        "random_state": scipy.stats.randint(0, 100)
    },
    DecisionTreeClassifier(): {
        "max_depth": scipy.stats.randint(1, 20),
    },
    RandomForestClassifier(): {
        "n_estimators": scipy.stats.randint(10, 100),
        "max_depth": scipy.stats.randint(1, 20),
    },
    KNeighborsClassifier(): {
        "n_neighbors": scipy.stats.randint(1, 20)
    }
}

# スコア比較用に変数を用意
max_score = 0
best_model = None
best_param = None

# グリッドサーチでパラメーターサーチ
for model, param in model_param_set_grid.items():
    clf = GridSearchCV(model, param)
    clf.fit(train_X, train_y)
    pred_y = clf.predict(test_X)
    score = f1_score(test_y, pred_y, average="micro")
    # 最高評価更新時にモデルやパラメーターも更新
    if max_score < score:
        max_score = score
        best_model = model.__class__.__name__
        best_param = clf.best_params_

# ランダムサーチでパラメーターサーチ
for model, param in model_param_set_random.items():
    clf = RandomizedSearchCV(model, param)
    clf.fit(train_X, train_y)
    pred_y = clf.predict(test_X)
    score = f1_score(test_y, pred_y, average="micro")
    # 最高評価更新時にモデルやパラメーターも更新
    if max_score < score:
        max_score = score
        best_model = model.__class__.__name__
        best_param = clf.best_params_

# 最適な学習モデル・パラメーター・ベストスコア        
print("学習モデル:{},\nパラメーター:{}".format(best_model, best_param))
print("ベストスコア:",max_score)

振り返り

当初の目的である「精神疾患の治療必要性を予測する」という点については、今回のアンケート結果から、かなり高い水準で実現可能であることがわかりました。今後新たな調査を実施する際には、今回のアンケートの設問設計が参考になりそうです。

また、今回の作業では、特にデータクレンジングの部分に時間を費やしました。「機械学習ではデータクレンジングに多くの時間が割かれる」という話がありますが、まさにそれを実感したところです。

メンタルヘルス、あるいは人事組織領域では、まだまだ機械学習を用いた分析が進んでいませんが、今後もより良い労働・職場環境づくりを支援するため、機械学習を用いたアプローチについて研鑽を積みたいと思います。

この記事が気に入ったらサポートをしてみませんか?