2  第2回 顧客層別の売上

2.1 はじめに

年代ごとの客層に違いがあるのかや、自社と他社の顧客層に違いがあるのかを調べる方法として、 クロス表の作成方法やカイ二乗検定の方法を学びます。

2.2 この章で使うファイルとパッケージ

第2回ファイルで使うデータは

  • chp2.xlsx

です。 ではまず必要なパッケージを読み込みます。

必要なパッケージの読み込み
options(scipen = 100) # 科学的表記を避ける
pacman::p_load(
  tidyverse,
  readxl,
  ggthemes,
  gt,
  gtExtras,
  janitor
  )

ここで新しいパッケージjanitorを読み込んでいます。 janitorはデータのクリーニングや集計を簡単にするためのパッケージで、特にクロス集計表を作成するための関数が充実しています。

ではreadxlパッケージを使ってデータchp2.xlsxを読み込んでみましょう。 前章と同様にこのExcelファイルのシートの一覧を表示してみます。

Excelファイルのシートの確認
readxl::excel_sheets("data/chp2.xlsx")
[1] "いつものPOSデータ"    "ピボットテーブル"     "表2-2・図2-2"        
[4] "表2-3・図2-4"         "表2-4・表2-5・図2-12" "表2-6・図2-13"       
[7] "表2-7・図2-14"        "表2-8・表2-10・検定"  "表2-11・表2-12・検定"

最初のシート「いつものPOSデータ」を読み込みたいので、シート番号を指定せずにread_excel()関数を使ってExcelファイルを読み込みます。 ついでに文字列をカテゴリー変数として因子型に変換しておきます

データの読み込みと変数の型変換
df <- readxl::read_excel("data/chp2.xlsx") # データを読み込む
df <- df |> # 文字列を因子型に変換
  mutate(曜日 = factor(曜日, levels = c("月", "火", "水", "木", "金", "土", "日"), ordered = TRUE),
         性別 = factor(性別, levels = c("男性", "女性"), ordered = TRUE),
         年代 = factor(年代, levels = c("20歳未満", "20代", "30代", "40代", "50代", "60歳以上"), ordered = TRUE),
         メーカー = as.factor(メーカー),
         商品名 = as.factor(商品名))
head(df) # 先頭6行を表示

2.3 自社商品における購入者の属性を調べる

2.3.1 セグメンテーション

顧客を属性で分類することをセグメンテーションといいます。 年代や性別といったセグメントごとの顧客の特徴を調べてみます。

まずは自社商品におおける年代別と性別ごとの販売個数を計算してみましょう。

自社商品の年代別・性別ごとの販売個数
df_jisha_age_sex_num <- df |>
  filter(メーカー == "自社") |> # 自社商品のみ抽出
  summarise(
    .by = c(年代, 性別),# 年代と性別ごとに
    販売個数合計 = sum(個数) # 個数を合計
  )

df_jisha_age_sex_num

janitorパッケージを使ってクロス集計表を作成することもできます。

自社商品の年代別・性別ごとの販売個数
df |>
  filter(メーカー == "自社") |>
  tabyl(年代, 性別) |> # 年代と性別のクロス集計表を作成
  adorn_rounding(digits = 0) |>
  gt() |>
  tab_header(title = "表2-1 性別・年代ごとの購入回数(クロス集計表)")
表2-1 性別・年代ごとの購入回数(クロス集計表)
年代 男性 女性
20歳未満 1583 4902
20代 4365 15130
30代 7568 25638
40代 6125 19134
50代 6982 19679
60歳以上 4064 10857

表が出力されましたが,いわゆる「ロング型」となっているため,表示に適した教科書のようなワイド型に変換してみましょう。 ついでに教科書のように,年代別合計を列に,性別ごとの合計を行に追加します。これが以外に面倒です。

自社商品の年代別・性別ごとの販売個数の表
total_sex <- df_jisha_age_sex_num |>
  summarise(
    .by = c(性別), # 性別ごとに
    合計 = sum(販売個数合計) # 個数を合計
    ) |>
  pivot_wider(
    names_from = 性別,
    values_from = 合計
  ) |> # ワイドに
  mutate(年代 = "合計") # 年代に"合計"を追加

df_jisha_age_sex_num |> # ワイドに
  pivot_wider(
    names_from = 性別,
    values_from = 販売個数合計
  ) |>
  rbind(total_sex) |> # 合計の行を追加
  mutate( # 合計の列を追加
    合計 = 女性 + 男性
  ) |>
  gt() |>
  fmt_number(columns = 2:4, decimals = 0) |>
  tab_header(title = "表2-1 性別・年代ごとの購入回数をまとめてクロス集計表")
表2-1 性別・年代ごとの購入回数をまとめてクロス集計表
年代 女性 男性 合計
50代 31,457 11,000 42,457
20代 24,164 6,979 31,143
60歳以上 17,291 6,464 23,755
30代 40,730 12,075 52,805
40代 30,477 9,826 40,303
20歳未満 7,838 2,540 10,378
合計 151,957 48,884 200,841

年代と性別の多重クロス表が完成しました。 販売個数ではなく,売上金額の表にしたい場合は,summarise(売上合計 = sum(金額))とすればよいです。

df_jisha_sale_age_sex <- df |>
  filter(メーカー == "自社") |>
  summarise(
    .by = c(年代, 性別), # 年代と性別ごとに
    売上合計 = sum(金額) # 金額を合計
  )

df_jisha_sale_age_sex |>
  pivot_wider(names_from = 性別, values_from = 売上合計) |>
  gt() |>
  fmt_number(columns = 2:3, decimals = 0) |>
  tab_header(title = "表2-2 性別・年代ごとの自社商品の売上金額")
表2-2 性別・年代ごとの自社商品の売上金額
年代 女性 男性
50代 4,718,550 1,650,000
20代 3,624,600 1,046,850
60歳以上 2,593,650 969,600
30代 6,109,500 1,811,250
40代 4,571,550 1,473,900
20歳未満 1,175,700 381,000

2.3.2 棒グラフによる可視化

可視化するため,棒グラフにしてみる。 geom_col()関数で棒グラフを描画し,position = "dodge"で性別ごとに横並びの棒グラフにします。

df_jisha_sale_age_sex |>
  ggplot() + aes(x = 年代, y = 売上合計, fill = 性別) +
  geom_col(position = "dodge") +
  labs(title = "図2-2 性別・年代ごとの自社商品の売上金額",
       x = "年代",
       y = "売上金額") +
  theme_bw(base_family = "HiraKakuPro-W3") +
  scale_fill_tableau(name = "Tableau 20")

積み上げ棒グラフにするなら,position = "stack"を指定します。

df_jisha_sale_age_sex |>
  ggplot() + aes(x = 年代, y = 売上合計, fill = 性別) +
  geom_col(position = "stack") +
  labs(title = "図2-2 自社商品の年代別・性別別売上金額",
       x = "年代",
       y = "売上金額") +
  theme_bw(base_family = "HiraKakuPro-W3") +
  scale_fill_tableau(name = "Tableau 20")

合計を1とした割合を示す積み上げ棒グラフにするなら,position = "fill"を指定します。

df_jisha_sale_age_sex |>
  ggplot() + aes(x = 年代, y = 売上合計, fill = 性別) +
  geom_col(position = "fill") +
  labs(title = "図2-2 自社商品の年代別・性別別売上金額",
       x = "年代",
       y = "売上金額") +
  theme_bw(base_family = "HiraKakuPro-W3") +
  scale_fill_tableau(name = "Tableau 20")

売上個数のグラフも作成してみます。

df |>
  filter(メーカー == "自社") |>
  summarise(
    .by = c(年代, 性別), # 年代と性別ごとに
    販売個数 = sum(個数) # 個数を合計
  ) |>
  ggplot() + aes(x = 年代, y = 販売個数, fill = 性別) +
  geom_col(position = "dodge") +
  labs(title = "図2-4 性別・年代ごとの自社商品の売上個数",
       x = "年代",
       y = "販売個数") +
  theme_bw(base_family = "HiraKakuPro-W3") +
  scale_fill_tableau(name = "Tableau 20")

2.3.3 割合を見る

表2-2の表を見ると,年代ごとの売上金額の合計が異なるため,単純に比較することができません。そこで売上高合計に占める割合を見てみます。

df_jisha_sale_age_sex |>
  mutate( # 割合を計算し,小数点以下2桁で丸める
    割合 = 売上合計 / sum(売上合計) * 100
  ) |>
  select(年代, 性別, 割合) |> # 必要な列だけ抽出
  pivot_wider(names_from = 性別, values_from = 割合) |>
  mutate( # 性別ごとの合計を追加
    合計 = 女性 + 男性
  ) |>
  gt() |>
  fmt_number(columns = 2:4, decimals = 2) |>
  tab_header(title = "表2-4 自社商品の売上金額における性別・年代の比率(%)")
表2-4 自社商品の売上金額における性別・年代の比率(%)
年代 女性 男性 合計
50代 15.66 5.48 21.14
20代 12.03 3.47 15.51
60歳以上 8.61 3.22 11.83
30代 20.28 6.01 26.29
40代 15.17 4.89 20.07
20歳未満 3.90 1.26 5.17

次に売上全体に対する割合ではなく,性別ごとの合計に対する割合を示す表を作成してみます。

df_row_ratio <- df_jisha_sale_age_sex |>
  pivot_wider(names_from = 性別, values_from = 売上合計) |>
  mutate(
    女性 = 女性 / sum(女性) * 100,
    男性 = 男性 / sum(男性) * 100
  ) |>
  select(年代, 女性, 男性)
df_row_ratio |>
  gt() |>
  fmt_number(columns = 2:3, decimals = 2) |>
  tab_header(title = "表2-5 各性別を100%とした場合の性別に対する年代の比率 (%)")
表2-5 各性別を100%とした場合の性別に対する年代の比率 (%)
年代 女性 男性
50代 20.70 22.50
20代 15.90 14.28
60歳以上 11.38 13.22
30代 26.80 24.70
40代 20.06 20.10
20歳未満 5.16 5.20

性別ごとの年代の比率から棒グラフを作成します。

df_row_ratio |>
  pivot_longer(cols = c(女性, 男性), names_to = "性別", values_to = "割合") |>
  ggplot() + aes(x = 年代, y = 割合, fill = 性別) +
  geom_col(position = "fill") +
  labs(title = "図2-9 帯グラフの作成",
       x = "年代",
       y = "割合") +
  theme_bw(base_family = "HiraKakuPro-W3") +
  scale_fill_tableau(name = "Tableau 20") +
  # グラフを横向き
  coord_flip()

次に,性別ごとの年代別の割合を見てみます。

df_row_ratio |>
  pivot_longer(cols = c(女性, 男性), names_to = "性別", values_to = "割合") |>
  ggplot() + aes(x = 性別, y = 割合, fill = 年代) +
  geom_col(position = "fill") +
  labs(title = "図2-10 帯グラフの作成",
       x = "性別",
       y = "割合") +
  theme_bw(base_family = "HiraKakuPro-W3") +
  #カラーパレットを変更
  scale_fill_tableau(name = "Tableau 20") +
  coord_flip()

2.3.4 年代別・性別ごとの購入者の属性を調べる

df_sales_total <- df |>
  dplyr::filter(メーカー == "自社") |>
  summarise(
    .by = c(年代, 性別),
    売上合計 = sum(金額)
  ) |>
  pivot_wider(names_from = 性別, values_from = c(売上合計))

df_sales_total |>
  rowwise() |>  # 行ごとに処理
  mutate(
    # 年代以外の列で,行を合計で割る
    across(-年代, ~ . / sum(c_across(-年代), na.rm = TRUE))
    ) |>
  gt() |>
  fmt_number(columns = 2:3, decimals = 2) |>
  tab_header(title = "表2-3 自社商品の年代別・性別別売上金額")
表2-3 自社商品の年代別・性別別売上金額
年代 女性 男性
50代 0.74 0.26
20代 0.78 0.22
60歳以上 0.73 0.27
30代 0.77 0.23
40代 0.76 0.24
20歳未満 0.76 0.24

2.4 2-2 競合商品との属性の違い

競合する商品の売り上げ個数を,性別・年代ごとに集計してみる。

df_maker_age_sex <- df |>
  summarise(
    .by = c(メーカー, 年代, 性別), # メーカー、年代、性別ごとに
    販売個数 = sum(個数) # 個数を合計
    )

競合A社の性別・年代ごとの売上個数を見てみる。

df_maker_age_sex |>
  filter(メーカー == "競合A")|>
  pivot_wider(names_from = 性別, values_from = 販売個数) |>
  select(-メーカー) |>
  mutate(
    合計 = 男性 + 女性
  ) |>
  gt() |>
  fmt_number(columns = 2:4, decimals = 0) |>
  tab_header(title = "表2-6 年代・性別ごとの競合A社の売上個数")
表2-6 年代・性別ごとの競合A社の売上個数
年代 女性 男性 合計
30代 35,578 12,788 48,366
60歳以上 33,680 13,267 46,947
50代 51,643 21,789 73,432
20代 22,661 8,235 30,896
40代 38,297 15,132 53,429
20歳未満 9,813 3,964 13,777
自社商品の年代別・性別ごとの購入回数のクロス集計表
df |>
  filter(メーカー == "自社") |>
  tabyl(年代, 性別) |>
  adorn_totals(where = c("row", "col")) |> # 行と列に合計を追加
  adorn_rounding(digits = 0) |>
  gt() |>
  tab_header(title = "表2-1 性別・年代ごとの購入回数(クロス集計表)")
表2-1 性別・年代ごとの購入回数(クロス集計表)
年代 男性 女性 Total
20歳未満 1583 4902 6485
20代 4365 15130 19495
30代 7568 25638 33206
40代 6125 19134 25259
50代 6982 19679 26661
60歳以上 4064 10857 14921
Total 30687 95340 126027

棒グラフにしてみる。

df_maker_age_sex |>
    filter(メーカー == "競合A")|>
    ggplot() + aes(x = 年代, y = 販売個数, fill = 性別) +
    geom_col(position = "dodge") +
    labs(title = "図2-13 競合A社の年代別・性別別売上個数",
         x = "年代",
         y = "販売個数") +
    theme_bw(base_family = "HiraKakuPro-W3") +
    scale_fill_tableau(name = "Tableau 20")

自社と競合A社を比較してみます。

# 自社とA社の比較
df_ji_A_item <- df_maker_age_sex |>
  filter(メーカー %in% c("自社", "競合A")) |> # 自社と競合A社のみ抽出
  pivot_wider(names_from = メーカー, values_from = 販売個数) |> # ワイド
  select(年代, 性別, 自社, 競合A) # 必要な変数のみ選択

# 全体売上個数を計算する
df_total_item <- df |>
  summarise(
      .by = c(年代, 性別), # 年代と性別ごとに
    合計 = sum(個数) # 売上個数を合計
  )

# 自社とA社のデータフレームと全体のデータフレームを結合
df_ji_A_item <- df_ji_A_item |>
  left_join(df_total_item, by = c("年代", "性別"))
# 作表
df_ji_A_item |>
  arrange(性別) |> # 性別で並び替えて
  group_by(性別) |> # 性別でグループ分けして
  gt() |> # 表にする。
  fmt_number(columns = 3:5, decimals = 0) |>
  tab_header(title = "表2-7 自社と競合A社の年代・性別別売上個数")
表2-7 自社と競合A社の年代・性別別売上個数
年代 自社 競合A 合計
男性
60歳以上 6,464 13,267 28,181
50代 11,000 21,789 46,809
20代 6,979 8,235 21,738
40代 9,826 15,132 35,443
30代 12,075 12,788 35,757
20歳未満 2,540 3,964 9,264
女性
30代 40,730 35,578 108,567
50代 31,457 51,643 119,413
60歳以上 17,291 33,680 72,066
20代 24,164 22,661 67,154
40代 30,477 38,297 97,851
20歳未満 7,838 9,813 24,904

割合棒グラフにしてみます。

df_ji_A_item |>
  pivot_longer(cols = -c(年代, 性別), names_to = "メーカー", values_to = "販売個数") |>
  mutate(
    年代性別 = paste0(年代, 性別),
    年代性別 = factor(年代性別, levels = rev(unique(年代性別))),
    メーカー = factor(メーカー, levels = c("自社", "競合A", "合計"))
    ) |>
  select(年代性別, メーカー, 販売個数) |>
  ggplot() + aes(x = メーカー, y = 販売個数, fill = 年代性別) + geom_col(position = "fill") +
  theme_bw(base_family = "HiraKakuPro-W3") +
  #カラーパレットを変更
  scale_fill_brewer(palette = "Paired") +
  coord_flip()

2.5 統計的推定

統計学において,関心のある調査対象全体を母集団(population)といいます。 通常,母集団は非常に大きいデータであるため,全てのデータを調査することはコストや手間の面で非常に困難です。 母集団の平均や分散といった統計量を母数(parameter)といい,通常母数は観察不可能な値と考えます。

そこで,母集団から一部を抽出して調査することにします。 この母集団から抽出された,母集団より数ケタ小さな一式のデータを標本(sample)といいます。 この標本の平均や分散といった集約した値を統計量(statistic)といい,標本平均や標本分散は計算可能な値です。 統計学では,全体を見ることなく,その一部である標本を見ることで,母集団の特性を推定することを目的としています。

2.6 カテゴリー変数間の関連性の検定

次に,ある特性をもつ2つの集団(カテゴリー)が関連しているかどうかを検定する方法を学びます。 カテゴリーに属するかということを表す変数をカテゴリー変数といいます。 たとえば,男なら1,女なら0といった性別カテゴリーの変数や,20代や50代といった年代に属するかどうかを表すカテゴリー変数があります。 このような2つのカテゴリー変数の関連性を調べるためにはクロス集計表を作ることが有益です。

例えば、1学年400名の生徒うち、男が245名、女が155名いるとしましょう。 その学年の中で、メガネをかけている男が41名、メガネを掛けている女が81名いました。 このクロス集計表は次のようなものになります。

メガネをかけている メガネをかけていない 合計
41 204 245
81 74 155
合計 122 278 400

この表から、メガネをかけている人の割合は、男性の中で 41\div 245 = 0.1673469、女性の中で 81 \div 155 = 0.52であることから、女子学生の方がメガネをかけている割合が高いことが分かりました。

2.6.1 カイ二乗検定

このクロス集計表から読み取れる関係が、統計的に意味があるのかどうかを調べるためには、\chi ^2(カイ二乗)検定を行います。 \chi^2検定は次のステップで実行します。

  1. 帰無仮説として、カテゴリー変数間に関連性はないと仮定
  2. その仮定のもとで、観測されたクロス集計表の度数が、理論的に予測される度数と大きく異なるかどうかを検定
  3. 予測される度数と観測された度数の差が大きいほど、帰無仮説が棄却され,対立仮説が採択される

\chi^2検定で用いられる統計量は、\chi^2統計量と呼ばれ、次の式で計算されます。

\chi^2 = \sum_{i=1}^n \sum_{j=1}^m \frac{(O_{ij} - E_{ij})^2}{E_{ij}} ここで、O_{ij}は観測された度数(観測度数)、E_{ij}は理論的に予測される度数(期待度数)です。nmはカテゴリー変数のカテゴリー数です。 つまり、2つのカテゴリー変数の関連性を調べる場合、\chi^2統計量は次のように計算されます。

\begin{aligned} \chi^2 &= \frac{(O_{11} - E_{11})^2}{E_{11}} + \frac{(O_{12} - E_{12})^2}{E_{12}} \\ &+ \frac{(O_{21} - E_{21})^2}{E_{21}} + \frac{(O_{22} - E_{22})^2}{E_{22}} \end{aligned}

ここで、期待度数Eをどうやって求めるのか、が問題となります。 期待度数の「期待」の意味は、帰無仮説のもとで期待される度数です。

E_{ij} = \frac{O_{i\cdot} \times O_{\cdot j}}{O_{\cdot \cdot}} ここで、O_{i\cdot}i行目の合計(横の合計)、O_{\cdot j}j列目の合計(縦の合計)、O_{\cdot \cdot}は全体の合計です

先のメガネの例で計算してみます。 観察度数Oは次のようになります。

メガネをかけている メガネをかけていない 合計
41 204 245
81 74 155
合計 122 278 400

男の行合計O_{男\cdot}は245、女の行合計O_{女\cdot}は155、メガネ有りの列合計O_{\cdot メガネ有}は122、メガネなしの列合計O_{\cdot メガネ無}は278、全体の合計O_{\cdot \cdot}は400となります。 ここから、期待度数は次のように計算されます。

\begin{aligned} E_{男, メガネ} &= \frac{245 \times 122}{400} = 74.725 \\ E_{男, メガネ無} &= \frac{245 \times 278}{400} = 170.275 \\ E_{女, メガネ} &= \frac{155 \times 122}{400} = 47.275 \\ E_{女, メガネ無} &= \frac{155 \times 278}{400} = 107.725 \end{aligned}

よって期待度数Eは次のようになります。

メガネをかけている メガネをかけていない 合計
74.725 170.275 245
47.275 107.725 155
合計 122 278 400

ここから、定義通りに、\chi^2統計量を計算します。

\begin{aligned} \chi^2 = \frac{(41 - 74.725)^2}{74.725} + \frac{(204 - 170.275)^2}{170.275} + \frac{(81 - 47.275)^2}{47.275} + \frac{(74 - 107.725)^2}{107.725} = 56.51731 \end{aligned}

ここで計算した\chi^2統計量56.51731は、自由度1の\chi^2分布に従うということが知られています。 この自由度は、カテゴリー変数のカテゴリー数から1を引いたものです。 ここでは、2カテゴリー同士のクロス集計表なので、自由度は(2-1) \times (2-1) = 1となります。

自由度1のカイ二乗分布の確率密度関数は次のようになります。

x = c(1:2500) / 250 # 数列を作成
y1 = dchisq(x,1) # カイ二乗分布の確率密度

df_2 <- data.frame(x, y1) # データフレームを作成
p <- ggplot(df_2) + aes(x = x,y = y1) + geom_line(linewidth = 1) # 折れ線グラフ
p <- p +
    scale_y_continuous(expand = c(0,0), limits = c(0,1)) +
    scale_x_continuous(expand = c(0,0), limits = c(-0.1,10)) +
    ylab("密度") + xlab("カイ二乗値") + ggtitle("自由度1のχ2分布の確率密度") +
    theme_bw(base_family = "HiraKakuPro-W3") +
    scale_fill_tableau(name = "Tableau 20")
print(p)
Warning: Removed 34 rows containing missing values or values outside the scale range
(`geom_line()`).

参考までに、自由度が変わると\chi^2分布の形状は次のようなものになります。

df_chi_graph <- data.frame(
  x = c(1:2500) / 250,
  y1 = dchisq(x,1), # 自由度1のカイ二乗分布の確率密度
  y3 = dchisq(x,3), # 自由度3のカイ二乗分布の確率密度
  y5 = dchisq(x,5) # 自由度5のカイ二乗分布の確率密度
)
df_chi_graph <- df_chi_graph %>%
  pivot_longer(names_to = "y",values_to = "value",cols = -x) # データを長く整形
p <- ggplot(df_chi_graph) + aes(x = x,y = value, group = y, color = y) # 作図
p <- p + geom_line(linewidth = 1) # 折れ線グラフ
p <- p + ylab("密度") + xlab("カイ二乗値") +
  scale_y_continuous(expand = c(0,0), limits = c(0,1)) +
  scale_x_continuous(expand = c(0,0), limits = c(-0.1,10))
p <- p + ggtitle("自由度1,3,5のχ2分布の確率密度") +
  theme_bw(base_family = "HiraKakuPro-W3") +
  scale_fill_tableau(name = "Tableau 20")
print(p)
Warning: Removed 34 rows containing missing values or values outside the scale range
(`geom_line()`).

自由度1の\chi^2分布における有意水準5%の値を調べるにはqchisq()関数を使います。引数は、pに確率、dfに自由度を指定します。

alpha <- 0.05  # 有意水準(ここでは5%)
degf <- 1        # 自由度
qchisq(1 - 0.05, 1)
[1] 3.841459
# qchisq(1 - alpha, degf)

自由度1のカイ二乗分布における有意水準5%の値は3.84であることが分かりました。 この値を超えると、有意水準5%で帰無仮説を棄却することになります。

では先程計算した\chi^2統計量は、有意水準5%で帰無仮説を棄却するかどうかを調べてみましょう。

chi2 <- 56.51731
qchisq(1 - alpha, degf) < chi2
[1] TRUE

より、 \chi^2 統計量 56.51731 は有意水準5%で帰無仮説を棄却することが分かりました。 ちなみに、自由度1のカイ二乗分布の確率密度関数と\chi^2統計量の位置を重ねてみると次のようになります。

df_chi_graph |>
  pivot_wider(names_from = "y", values_from = "value") |> # ワイド型に変換
  # 折れ線グラフを作成
  ggplot() + aes(x = x, y = y1) + geom_line(linewidth = 1) +
  geom_vline(xintercept = chi2, linetype = "dashed", color = "red") +
  annotate("text", x = chi2, y = 0.1, label = "χ2 statistics", color = "red") +
  # グラフの見た目の設定
  ggtitle("自由度1のχ2分布の確率密度") +  # タイトル
  ylab("密度") + xlab("カイ二乗値") + # 軸ラベル
  scale_y_continuous(expand = c(0,0), limits = c(0,1)) +
  scale_x_continuous(expand = c(0,0), limits = c(-0.1,60)) +
  theme_bw(base_family = "HiraKakuPro-W3") + # テーマの設定
  scale_fill_tableau(name = "Tableau 20") # カラーパレットの設定

このように、\chi^2統計量は、自由度1のカイ二乗分布のもとで生じる確率は、

1 - pchisq(chi2, df = 1)
[1] 0.0000000000000557332

となり、非常に小さな値であることが分かりました。 つまり、2つのカテゴリー変数の間に関係がない、という帰無仮説の下で、観測された度数が発生することはほぼありえない、ということが言えるので、帰無仮説は棄却され、2つのカテゴリー変数には関係があると結論付けられます。

2.6.2 Rでχ2検定

Rではchisq.test()関数を使ってχ2検定を行うことができます。 引数は、xに度数表、correctに補正を行うかどうか、pに期待度数を指定します。

先ほどの男女とメガネの例をここでも使ってみます。 まずmatrix()関数を使ってクロス集計表を行列として作成します。

O <- matrix(c(41, 81, 204, 74), nrow = 2, ncol = 2)
row.names(O) <- c("男性", "女性")
colnames(O) <- c("メガネ", "メガネなし")

E <- matrix(c(
  sum(O[,1])*sum(O[1,])/sum(O), #男眼鏡
  sum(O[,1])*sum(O[2,])/sum(O), #男眼鏡無
  sum(O[,2])*sum(O[1,])/sum(O), #女眼鏡
  sum(O[,2])*sum(O[2,])/sum(O)  #女眼鏡無
), nrow = 2, ncol = 2)

print(O)
     メガネ メガネなし
男性     41        204
女性     81         74
print(E)
       [,1]    [,2]
[1,] 74.725 170.275
[2,] 47.275 107.725

この観察度数と期待度数から、定義通りに\chi^2統計量を計算してみます。

chi <-  ((O[1,1] - E[1,1])^2 / E[1,1]) + #男眼鏡
        ((O[1,2] - E[1,2])^2 / E[1,2]) + #男眼鏡無
        ((O[2,1] - E[2,1])^2 / E[2,1]) + #女眼鏡
        ((O[2,2] - E[2,2])^2 / E[2,2])   #女眼鏡無
print(chi)
[1] 56.51731

\chi^2 統計量が 56.5173098 となりました。 この \chi^2 統計量が自由度1の \chi^2 分布にしたがう場合,この統計量が得られる確率は次のようになります。

prop <- 1 - pchisq(chi, df = 1)
print(prop)
[1] 0.0000000000000557332

この確率は0.0000000000000055733となり,ほぼゼロであることが分かりました。 よって、2つのカテゴリー変数は無関係である,という帰無仮説は棄却され、2つのカテゴリー変数には関係があると結論付けられます。

ちなみに,上記のようなめんどくさい処理をしなくても,Rにはchisq.test()という関数が用意されています。 chisq.test()は引数として、xに度数表、correctに補正を行うかどうか、pに期待度数を指定します。 補正は行わないので,correctFALSEとします。 pはデフォルトで等確率となっているので,今回は省略します。

chisq.test(O, correct = FALSE)

    Pearson's Chi-squared test

data:  O
X-squared = 56.517, df = 1, p-value = 0.00000000000005571

となり,先ほどの結果と一致しました。

各マスに入る度数が少ない場合には、フィッシャーの直接確率検定を使います。

fisher.test(O)

    Fisher's Exact Test for Count Data

data:  O
p-value = 0.0000000000001204
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.1127341 0.2982245
sample estimates:
odds ratio 
 0.1845096 

2.6.3 自社商品の売上個数における性別と年代の差

表2−3のクロス集計表を標本とみなして、この標本から母集団の性別・年代ごとに売上個数の差があるのかどうかを検証します。 ここでは、

  • H_0: 性別と年代による売上個数の差はない
  • H_1: 性別と年代による売上個数の差はある

この仮説を検定するために、\chi^2検定を行います。 年代別・性別ごとの売上個数のクロス集計表を作成します。

df_chi2 <- df |>
  filter(メーカー == "自社") |>
  summarise(
    .by = c(年代, 性別), # 年代と性別ごとに
    販売個数 = sum(個数)
  ) |>
  pivot_wider(names_from = 性別, values_from = 販売個数)

df_chi2data.frame型なので、そのままではchisq.test()関数に渡すことができません。 そのため、matrix型に変換します。

mat_chi2 <- as.matrix(df_chi2[,2:3])
mat_chi2
      女性  男性
[1,] 31457 11000
[2,] 24164  6979
[3,] 17291  6464
[4,] 40730 12075
[5,] 30477  9826
[6,]  7838  2540

作成した行列型mat_chi2を使って\chi^2検定を行います。

chisq.test(mat_chi2, correct = FALSE)

    Pearson's Chi-squared test

data:  mat_chi2
X-squared = 288.42, df = 5, p-value < 0.00000000000000022

\chi^2 値が288.42となり、 p 値が0.00000000000000022とほぼゼロとなっていることから、 「自社と競合A社の年代・性別別売上個数に差がない」という帰無仮説は棄却されました

トップに戻る