見出し画像

「入門はじめての多変量解析」をPythonで写経 Vol.28 ~ 7章「はじめての数量化Ⅰ類」数量化Ⅰ類とカテゴリ数量

7章「はじめての数量化Ⅰ類」

書籍の著者 石村貞夫 先生、石村光資郎 先生


書籍「入門はじめての多変量解析」7章「はじめての数量化Ⅰ類」の Python写経活動記録 です。 

多変量解析の入門を Python と一緒に学ぶ写経シリーズです。

この記事は、数量化Ⅰ類によるアンケート分析を実践します。
数量化Ⅰ類はざっくり「説明変数がカテゴリ変数、目的変数が量的変数」のデータ分析手法です。
目的変数(外的基準)を予測する「予測式」の推定などを行います。

ダミー変数を用いる重回帰分析と似ているようで少々異なる不思議な分析手法。

ChatGPT 活用型学習で進めてまいります!
では書籍を開いて多変量解析の旅に出かけましょう🚀

子供達の飛行機旅行のイラスト(修学旅行):「いらすとや」さんより

はじめに


このブログシリーズは、書籍「入門はじめての多変量解析」(東京図書、「テキスト」と呼びます)の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

7章 はじめての数量化Ⅰ類


この記事は7章の以下のSectionを取り扱います。

7.1 数量化理論とは?
7.2 予測に役立つ数量化Ⅰ類
7.3 カテゴリ数量の基準化
7.4 数量化Ⅰ類で・・・をあてよう

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

この記事で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd
import sympy

# 重回帰
from sklearn.linear_model import LinearRegression

# LaTeX表示
from IPython.display import Math

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

数量化理論とは?

テキストは「アンケート調査の回答のような質的データに対しても、最適な数量や評点を与えて分析しようとするのが、数量化理論です」と説明します。
確かにテキストの例題データはすべて、アンケート・回答形式です。

数量化理論はざっくり、説明変数などに含まれる質的変数(カテゴリ変数)を数値化して分析を進める手法 です。
アンケートの回答選択肢「はい、いいえ」や「そう思わない、どちらとも言えない、そう思う」などは質的変数です。

これまでの記事で実践した解析手法と比較すると…

$$
\begin{array}{lll}
目的変数 & 数量化理論 & 類似解析手法 \\
\hline
 \\
量的変数 & 数量化Ⅰ類 & 重回帰分析 \\
質的変数 & 数量化Ⅱ類 & 判別分析 \\
なし & 数量化Ⅲ類 & 主成分分析 \\
\end{array}
$$

■ 数量化理論の用語
例題データを見ながら、数量化理論の用語を確認しましょう。
テキスト p.242 表 7.2.1 の例題データを引用します。

### アンケート調査の結果 p.243 表7.2.2

# データフレームの作成
data1 = pd.DataFrame(
    np.array([[57, 65, 51, 54, 45, 67],
              [1, 1, 0, 1, 0, 1],
              [0, 0, 1, 0, 1, 0],
              [0, 1, 1, 1, 1, 0],
              [1, 0, 0, 0, 0, 1]]).T,
     columns=[['外的基準', '野菜', '野菜', 'タンパク質', 'タンパク質'],
              ['', 'はい', 'いいえ', 'はい', 'いいえ']],
     index=range(1, 7))
data1.index.name = '被験者No.'

# 結果の表示
data1

【実行結果】
「はい、いいえ」の選択肢には該当する方に 1 、該当しない方に 0 が設定されています。

数量化理論の用語を一覧化します。

$$
\begin{array}{lll}
数量化理論の用語 & 言い換え & 例題データ \\
\hline
 \\
外的基準 & 目的変数 & 仮想の健康スコア\\
アイテム & 説明変数  & 野菜・ タンパク質 \\
カテゴリ & カテゴリ値 & はい・いいえ \\
カテゴリ数量 & 偏回帰係数 & (b_0, b_{11}, b_{12}, b_{21}, b_{22}) \\
\end{array}
$$

数量化Ⅰ類の予測式

  • 外的基準: $${Y}$$

  • アイテム $${i}$$ のカテゴリ $${j}$$ のダミー変数: $${x_{ij}}$$

    • 取りうる値は $${0}$$ か $${1}$$

  • アイテム $${i}$$ のカテゴリ $${j}$$ のカテゴリ数量: $${b_{ij}}$$

とすると、数量化Ⅰ類の予測式は次の1次式で表現できます。

📊 数量化Ⅰ類の予測式(アイテムは2つ、カテゴリは各2つ)

$$
Y = b_{11} x_{11} + b_{12} x_{12} + b_{21} x_{21} + b_{22} x_{22} + b_0
$$

テキストの数式を引用

例題データの列見出しを $${y}$$、$${x_{ij}}$$ に変換しましょう。

### データセットの作成:マルチカラムを解除したデータフレームの作成
data1_cat = data1.copy()
data1_cat.columns = ['y', 'x11', 'x12', 'x21', 'x22']
data1_cat

【実行結果】

基準化前のカテゴリ数量を求める

重回帰分析では、目的変数の実測値と予測値の差=残差の二乗和を最小化する偏回帰係数を求めました。
数量化Ⅰ類の場合も、外的基準(目的変数)の実測値と予測値の差の二乗和を最小化するようにカテゴリ数量(係数)を推定します。

$$
「Q = \sum_{i=1}^N (外的基準\ y_i - 予測値\ Y_i)^2」を最小にする\ b_{ij}
$$

テキストの数式を引用

ここからは sympy ライブラリを利用して、次の4ステップで「基準化前のカテゴリ数量」を推定します。

  1. 差の2乗和 $${Q}$$ の算出

  2. $${Q}$$ を $${b_{ij}}$$ で偏微分

  3. 偏微分結果に$${b_{i1} = 0}$$ を代入

  4. 代入後の式を0とおき、連立方程式を求解

1️⃣ 差の2乗和 $${Q}$$ の算出
テキスト p.245 の数式に相当します。

### 数量化Ⅰ類の予測式の求め方 p.244~

# 変数の定義
b0, b11, b12, b21, b22 = sympy.symbols('b0 b11 b12 b21 b22')

# 2次式Qの設定
Q = sum([(y - (b11*x11 + b12*x12 + b21*x21 + b22*x22 + b0))**2
         for y, x11, x12, x21, x22 in data1_cat.values])

# 結果の表示
print('【2次式Q】')
display(Math(f'Q = {sympy.latex(Q.simplify())}'))

【実行結果】

2️⃣ $${Q}$$ を $${b_{ij}}$$ で偏微分
テキスト p.246 の数式に相当します。

# Qをb11, b12, b21, b22で偏微分
expr_b11 = sympy.diff(Q, b11)
expr_b12 = sympy.diff(Q, b12)
expr_b21 = sympy.diff(Q, b21)
expr_b22 = sympy.diff(Q, b22)

# 結果の表示
print('【2次式Qを偏微分】')
display(Math(f'\cfrac{{\partial Q}}{{\partial b_{11}}} = {sympy.latex(expr_b11)}'))
display(Math(f'\cfrac{{\partial Q}}{{\partial b_{12}}} = {sympy.latex(expr_b12)}'))
display(Math(f'\cfrac{{\partial Q}}{{\partial b_{21}}} = {sympy.latex(expr_b21)}'))
display(Math(f'\cfrac{{\partial Q}}{{\partial b_{22}}} = {sympy.latex(expr_b22)}'))

【実行結果】

3️⃣ 偏微分結果に$${b_{i1} = 0}$$ を代入
テキスト p.247 の数式に相当します。

# 偏微分した式にb11=0, b21=0を代入

# 代入する値をリスト化
sub_values = [(b11, 0), (b21, 0)]

# 代入
expr_b11 = expr_b11.subs(sub_values)
expr_b12 = expr_b12.subs(sub_values)
expr_b21 = expr_b21.subs(sub_values)
expr_b22 = expr_b22.subs(sub_values)

# 結果の表示
print('【b11=0, b21=0を代入】')
display(Math(f'{sympy.latex(expr_b11)} = 0'))
display(Math(f'{sympy.latex(expr_b12)} = 0'))
display(Math(f'{sympy.latex(expr_b21)} = 0'))
display(Math(f'{sympy.latex(expr_b22)} = 0'))

【実行結果】

4️⃣ 代入後の式を0とおき、連立方程式を求解
テキスト p.247 の数式に相当します。

# 連立方程式を解く
solved = sympy.solve([expr_b11, expr_b12, expr_b21, expr_b22])
# 値の取り出し
coefs = {k: float(v) for k, v in solved.items()}
# 結果の表示
display(Math(f'{sympy.latex(coefs)}'))

【実行結果】
$${b_0=59.5,\ b_{12}=-11.5,\ b_{22}=2.5}$$ と求まりました。

予測式の形式で表示しましょう。

### 予測式Y p.247
display(Math(f'Y={coefs[b12]:+} x_{{12}}{coefs[b22]:+}x_{{22}}{coefs[b0]:+}'))

【実行結果】

このままでは $${x_{11}, x_{21}}$$ に対応するカテゴリ数量 $${b_{11}, b_{21}}$$ が 0 のままです…
カテゴリ数量の基準化を通じて、予測式を仕上げていきましょう!

カテゴリ数量の基準化

■ 数量化Ⅰ類の予測式の完成
テキストによると「カテゴリ数量の基準化は、各アイテム内のカテゴリ数量の平均が0となるように、カテゴリ数量を変換すること」です。
変換のイメージはぜひテキストをお読み下さい!

ここでは、数量化Ⅰ類関数を作って、基準化済みのカテゴリ数量等の計算を実施しましょう。
この関数は次の3ステップでカテゴリ数量を計算しています。

  1. まずカテゴリ数量に重回帰分析の偏回帰係数を仮置きする。

  2. 次にカテゴリ数量の平均値を求める。

  3. 仮置きした偏回帰係数からカテゴリ数量を差し引いて基準化されたカテゴリ数量を算出する。

### 数量化Ⅰ類関数の定義 p.248, 252
# num_catsの設定方法:各アイテムのカテゴリの数を設定する
# [[item1のcategory1, ・・・, category_n], ・・・, [item_pのcategory1, ・・・, category_m]]
# のとき、[n, ・・・, m]と設定

def quantification_method_type1(X, y, num_cats):
    
    ## 設定と準備
    # num_cats[2, 3, ...]を[0, 1, 0, 1, 1,...]に変換
    cats_flat = np.array([0 if i==0 else 1 for n in num_cats for i in range(n)])
    # アイテム単位の0番以外のカテゴリ数をカウント
    cumsum_means = np.cumsum([cat - 1 for cat in num_cats]).tolist()
    # 各アイテムの0番目ダミー変数を除外して説明変数を作成
    X_dummy = X[X.columns[cats_flat != 0]]
    
    ## 重回帰分析の実行
    reg = LinearRegression()
    reg.fit(X_dummy, y)

    ## カテゴリ数量の平均値の算出
    # 「回帰係数 × ダミー変数の平均値」で回帰係数=カテゴリ数量を変形
    ratios = reg.coef_ * X_dummy.mean(axis=0).values 
    # カテゴリ数量の平均値の算出
    cat_num_means = [sum(ratios[i:j])
                     for i, j in zip([0] + cumsum_means[:-1], cumsum_means)]
    
    ## 基準化されたカテゴリ数量の算出
    # 設定と準備
    coef_pointer = 0   # 回帰係数のポインタ
    coefs = []         # 基準化されたカテゴリ数量を格納するリスト
    ranges = []        # 基準化されたカテゴリ数量の範囲を格納するリスト
    
    # アイテムごとにカテゴリ数量(回帰係数)からカテゴリ数量の平均値を減算
    for num_cat, cat_num_mean in zip(num_cats, cat_num_means):
        # 最小値・最大値の初期値設定(範囲算出用)
        min_val, max_val = np.inf, -np.inf
        # アイテム内のカテゴリごとに繰り返し処理
        for i in range(num_cat):
            # 基準化されたカテゴリ数量の算出
            if i == 0:  # 0番目のカテゴリ:0 - カテゴリ数量の平均値
                coefs.append(0 - cat_num_mean)
            else:       # 1番目以降のカテゴリ:カテゴリ数量 - カテゴリ数量の平均値
                coefs.append(reg.coef_[coef_pointer] - cat_num_mean)
                coef_pointer += 1
            # 基準化されたカテゴリ数量の最小値・最大値を更新(範囲算出用)
            if min_val > coefs[-1]:
                min_val = coefs[-1]
            if max_val < coefs[-1]:
                max_val = coefs[-1]
        # 範囲の算出
        ranges.append(max_val - min_val)

    ## 切片の算出
    intercept_based = y.mean()

    ## 戻り値: 基準化されたカテゴリ数量、切片、範囲、線形回帰モデル
    return {'coef': np.array(coefs), 'intercept': intercept_based,
            'range': np.array(ranges), 'reg': reg}

【実行結果】なし

テキスト p.248 のカテゴリ数量の基準化と同じ「予測式」を求めてみましょう。

### カテゴリ数量の基準化 p.248

# 説明変数と目的変数の設定
X = data1_cat.iloc[:, 1:]
y = data1_cat.iloc[:, 0]

# アイテムごとのカテゴリ数リストの設定
num_cats = [2, 2]

# 数量化Ⅰ類の実行:基準化されたカテゴリ数量の算出
res = quantification_method_type1(X, y, num_cats)

# 結果の表示
display(Math(
    f"Y = {res['coef'][0]:.3f} x_{{11}} {res['coef'][1]:+.3f} "
    f"x_{{12}} {res['coef'][2]:+.3f} x_{{21}} {res['coef'][3]:+.3f} "
    f"x_{{22}} {res['intercept']:+.1f}"))

【実行結果】
無事に予測式 $${Y}$$ を算出できました!

■ アイテムの影響の度合い「範囲」の確認
各アイテム(説明変数)が外的基準(目的変数)の予測に与える影響の大きさを「範囲」という指標で確認します。
範囲は各アイテム(説明変数)のカテゴリ数量(係数)の最大値・最小値の差です。

さきほどの関数の戻り値から各アイテムの範囲を表示しましょう。
テキスト p.252 の表形式を先取りします。

### 数量化Ⅰ類の結果の表示 p.252 表7.5.3
data1_result = data1.iloc[:, 1:].sum(axis=0).rename('度数').to_frame()
data1_result['カテゴリ数量'] = res1['coef']
data1_result.loc[(slice(None), 'はい'), '範囲'] = res1['range']
data1_result.round(3)

【実行結果】
健康スコアに与える影響は「野菜」の方が大きい、という結果になりました。

■ 学習データの予測値
予測式ができましたので、学習データを用いて予測をしてみましょう。
予測・評価用のヘルパー関数を定義します。

### モデルの評価 p.252

# 数量化Ⅰ類モデルで予測・評価する関数の定義
def predict_quantification_method_type1(X, y=None, res=None):
    # 予測値の算出
    y_pred = X @ res['coef'] + res['intercept']
    # 予測のみの場合(yの実測値=Noneの場合)、予測値を返す
    if y is None:
        return {'y_pred': y_pred}
    # 評価をする場合、重相関係数、決定係数、予測値を返す
    else:
        # 重相関係数の算出
        multi_corr_coef = np.corrcoef(y, y_pred)[0, 1]
        # 決定係数の算出
        R2 = multi_corr_coef**2
        # 戻り値:
        return {'重相関係数': multi_corr_coef, '決定係数': R2, 'y_pred': y_pred}

【実行結果】なし

予測を実行します。

## 学習データに対する予測

# 予測の実行
preds1 = predict_quantification_method_type1(X, y, res1)

# データフレーム化
data1_pred = data1['外的基準'].rename('実測値').to_frame()
data1_pred['予測値'] = preds1['y_pred']
data1_pred['残差'] = data1_pred['実測値'] - data1_pred['予測値']
data1_pred

【実行結果】

予測値は当たっているのか外しているのか…どちらでしょう?
テキスト p.252 を先取りして、重相関係数と決定係数で評価してみます。

# 学習データによるモデルの評価 p.252
print(f"重相関係数: {preds1['重相関係数']:.4f}, "
      f"決定係数: {preds1['決定係数']:.4f}")

【実行結果】
決定係数 $${0.6344}$$ は、モデルがデータにまあまあ当てはまっていると思います。

■ (ご参考)オンラインで数量化Ⅰ類を実践!
数量化Ⅰ類をはじめ、いろんなデータ分析手法をオンライン上で動かせるWebサイトのご紹介です。
青木繁伸 先生のホームページの「Black-Box」です。
※なお https 通信に対応していません。

Webサイトの画面を引用いたします。
「分析プロシージャ」にオンラインで分析可能な手法が列挙されています。
「数量化Ⅰ類」などの数量化理論も分析できます。

http://aoki2.si.gunma-u.ac.jp/BlackBox/BlackBox.html より引用

まず、データをカンマ区切りのテキストファイルにしておき、「分析データファイルの指定」でファイルをアップロードします。

続いて、「分析プロシージャの指定」(手法を選択)と、「分析オプションの指定」(分析に必要な指示)を設定して、「結果を表示する」ボタンを押下します。
詳しい使い方は画面上部の「使用法」をご覧ください。

例題データの分析結果は次のようになりました。
変数名は次のとおりです。
・score:健康スコア(目的変数)
・vege_no:野菜・いいえ
・prot_no:タンパク質・いいえ

詳細な統計量を簡単に取得できるので、データ分析が捗る感じがいたします。

例題を分析

テキスト p.250~の例題データを数量化Ⅰ類で分析します。

1️⃣ データの準備
テキスト p.251 表 7.5.2 の例題データを引用します。

### アンケート調査の結果 p.251 表7.5.2

# データの登録
data2 = pd.DataFrame(np.array(
    [[79, 82, 86, 78, 90, 83, 75, 81, 73, 82, 80, 88, 81, 85, 78, 82, 84, 82, 80],
     [1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1],
     [0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0],
     [0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1],
     [0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0],
     [1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0],
     [1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1],
     [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
     [1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1],
     [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0],
     [1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1],
     [0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0],
     [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]).T,
    index=range(1, 20),
    columns=[['外的基準', '出身地', '出身地', 'スポーツ', 'スポーツ', 'スポーツ',
              '野菜', '野菜', 'タンパク質', 'タンパク質', '牛乳', '牛乳', '牛乳'],
             ['', '1', '2', '1', '2', '3', '1', '2', '1', '2', '1', '2', '3']])
# 結果の表示
data2

【実行結果】
外的基準=仮想の健康スコア、5つのアイテムで構成する標本サイズ 19 のアンケート回答データです。

変数名を変更しましょう。
数量化Ⅰ類のモデル構築ではこのデータをフィッティングします。

### データセットの作成:マルチカラムを解除したデータフレームの作成
data2_cat = data2.copy()
data2_cat.columns = ['y', 'x11', 'x12', 'x21', 'x22', 'x23', 'x31', 'x32',
                     'x41', 'x42', 'x51', 'x52', 'x53']
data2_cat

【実行結果】
アイテム(説明変数)はダミー変数化されています。

2️⃣ データの確認
データの取り扱いが容易な形式に変換します。

### データを可視化で扱いやすい形式に整える

# カテゴリが3つあるアイテムの処理関数
def make_data(row):
    if row.iloc[0] == 1:
        return 0
    elif row.iloc[1] == 1:
        return 1
    elif row.iloc[2] == 1:
        return 2

# データフレームの作成
data2_plot = pd.DataFrame()
data2_plot['スコア'] = data2_cat['y']
data2_plot['出身地'] = data2_cat['x12']
data2_plot['スポーツ'] = data2_cat[['x21', 'x22', 'x23']].apply(make_data, axis=1)
data2_plot['野菜'] = data2_cat['x32']
data2_plot['タンパク質'] = data2_cat['x42']
data2_plot['牛乳'] = data2_cat[['x51', 'x52', 'x53']].apply(make_data, axis=1)

# 結果の表示
data2_plot

【実行結果】
アイテム(説明変数)を1つにまとめています。

相関係数をみましょう。

### 相関係数
data2_plot.corr().round(3)

【実行結果】
外的基準:健康スコアとある程度の相関関係があるのは「牛乳」と「タンパク質」です。

箱ひげ図で健康スコアとアイテムの関係を見てみましょう。

### 箱ひげ図

# 設定
labels = [['大都会', 'それ以外'], ['よくする', '時々', 'しない'], ['好き', '嫌い'],
          ['好き', '嫌い'], ['よく飲む', '時々', '飲まない']]
# 描画領域の設定
fig, axes = plt.subplots(2, 3, figsize=(8, 6), sharey=True, tight_layout=True)
# アイテムごとに箱ひげ図を繰り返し描画
for col, label, ax in zip(data2_plot.columns[1:], labels, axes.flat):
    # 箱ひげ図の描画
    sns.boxplot(data=data2_plot, y='スコア', x=col, fill=False, gap=0.1, ax=ax)
    # スウォームプロットの描画
    sns.swarmplot(data=data2_plot, y='スコア', x=col, ax=ax)
    # x軸目盛りラベルの設定
    ax.set_xticks(ticks=range(len(label)), labels=label)
# 最後の空グラフを非表示
axes[-1, -1].axis('off');

【実行結果】
各アイテム(説明変数)のカテゴリ値によって、健康スコアとの関係の相違がありそうですね!

3️⃣ 数量化Ⅰ類のモデル構築 
先ほどの数量化Ⅰ類関数を利用します。

### 数量化Ⅰ類の実行 p.252

# 説明変数と目的変数の設定
X = data2_cat.iloc[:, 1:]
y = data2_cat.iloc[:, 0]
# アイテムとカテゴリの関係リストの設定
categories = [2, 3, 2, 2, 3]

# 数量化Ⅰ類の実行~基準化されたカテゴリ数量の算出
res2 = quantification_method_type1(X, y, categories)
pprint.pprint(res2)

【実行結果】
coef: カテゴリ数量 $${b_{ij}}$$、intercept: $${b_0}$$、range: 範囲、reg: 線形回帰モデル が返ってきました。

4️⃣ 数量化Ⅰ類の結果の表示
テキスト p.252 表 7.5.3 に相当します。

### 数量化Ⅰ類の結果の表示 p.252 表7.5.3
data2_result = data2.iloc[:, 1:].sum(axis=0).rename('度数').to_frame()
data2_result['カテゴリ数量'] = res2['coef']
data2_result.loc[(slice(None), '1'), '範囲'] = res2['range']
data2_result.round(3)

【実行結果】
範囲を確認しましょう。
外的基準(目的変数)への影響の大きいアイテム(説明変数)は「牛乳」と「スポーツ」のようです。

5️⃣ モデルの評価
テキストにならって重相関係数と決定係数を算出します。
テキスト p.252 掲載の各評価指標に相当します。
評価用のヘルパー関数を利用します。

# 学習データによるモデルの評価 p.252
evals2 = predict_quantification_method_type1(X, y, res2)
print(f"重相関係数: {evals2['重相関係数']:.4f}, "
      f"決定係数: {evals2['決定係数']:.4f}")

【実行結果】
テキストとほぼ一致しています。

学習データによる予測値を見ておきます。

# 学習データによる予測値のデータフレーム化
data2_eval = data2['外的基準'].rename('実測値').to_frame()
data2_eval['予測値'] = evals2['y_pred']
data2_eval['残差'] = data2_eval['実測値'] - data2_eval['予測値']
data2_eval

【実行結果】

6️⃣ 数量化Ⅰ類モデルによる予測
次のアイテム情報を予測式に入れて、健康スコアを予測します。

$$
\begin{array}{}
アイテム & カテゴリ & 変数値 \\
\hline
 \\
出身地 & 都会 & x_{11}=1 \\
スポーツ & ときどき & x_{22}=1 \\
野菜 & 好き & x_{31}=1 \\
タンパク質 & 嫌い & x_{42}=1 \\
牛乳 & よく飲む & x_{51} = 1 \\
\end{array}
$$

### テキストの予測を実行 p.253

# xの値の設定 
x_pred = np.array([1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0])
# 予測の実行
preds2 = predict_quantification_method_type1(x_pred, None, res2)
# 予測値の表示
preds2['y_pred']

【実行結果】
健康スコアは $${85.057}$$ です。

7️⃣ 重回帰モデルによる予測
p.251 でうさぎさんが
「実は重回帰分析をおこなっても予測値は同じ結果になります」
とおっしゃっています。

重回帰分析を実行して比べてみましょう。
実は数量化Ⅰ類関数の戻り値に学習済みの線形回帰モデルが含まれています。
このモデルで予測しましょう。

### 線形回帰モデルで予測

## 予測用のデータの作成
# 除外する変数のリストの作成
cats_flat = np.array([0 if i==0 else 1 for n in categories for i in range(n)])
# 予測用の説明変数の作成
X_dummy = pd.DataFrame(x_pred[cats_flat != 0], index=X.columns[cats_flat != 0]).T
# 予測の実行
res2['reg'].predict(X_dummy)[0]

【実行結果】
重回帰分析の結果が数量化Ⅰ類の予測値と同じになることを確認できました。

記事の最後はChatGPTが締めくくります。
カテゴリに「かたち」を見いだして。

📘 ChatGPTのひとこと:

今回の記事では、数量化Ⅰ類を通じて「カテゴリデータを数値に置き換えて見る視点」を身につけました。まるで日々の買い物で色とりどりの食材を改めて分類し、棚に並べ直すように、データの“かたち”を新たに見つめ直せたのではないでしょうか😊

カテゴリを数値のレンズで捉えることで、分析の地平がぐっと広がります。
次回は数量化Ⅱ類で、複数のカテゴリ変数を同時に扱う手法を学び、一歩深い洞察の世界へ進みましょう。

静かな好奇心を胸に、また一緒に学びを重ねられることを楽しみにしています🌱

今回の写経は以上です。


シリーズの記事

次の記事

前の記事

目次

ブログの紹介


note で7つのシリーズ記事を書いています。
ぜひ覗いていってくださいね!

1.のんびり統計

統計検定2級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定2級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

2.実験!たのしいベイズモデリング1&2をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング2」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは R言語+Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう!

3.実験!岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験!岩波データサイエンスvol.1」の4人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう!

4.楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

5.RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん!
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

6.データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

7.Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。

いいなと思ったら応援しよう!

ネイピア DS 応援ありがとうございます。これからもがんばって記事を作成します!

この記事が参加している募集