Elixir Livebook で地域経済分析システム RESAS から産業構造を取得して分析する

## はじめに

RESAS は日本の地域経済データを取得できるサービスです

https://resas.go.jp/#/44/44201

RESAS-API に利用登録すると、 REST API でデータを取得することもできます

https://opendata.resas-portal.go.jp/

前回は人口構成を可視化してみましたが、今回は産業構造を分析してみたいと思います

https://qiita.com/RyoWakabayashi/items/2c9d208e5055fe20f1e3

今回実装したノートブックの全量はこちら

https://github.com/RyoWakabayashi/elixir-learning/blob/main/livebooks/resas/industry.livemd

## 出典

この記事はRESAS（地域経済分析システム）のデータを加工して作成しています

## 実行環境

このリポジトリーの Docker コンテナ上で実行しました

https://github.com/RyoWakabayashi/elixir-learning

## RESAS-API の使い方

前回の記事を参照してください

https://qiita.com/RyoWakabayashi/items/2c9d208e5055fe20f1e3#resas-api-%E3%81%AE%E4%BD%BF%E3%81%84%E6%96%B9

## セットアップ

Livebook を起動してノートブックを開きます

以下のコードをノートブック上で実行してください

```elixir
Mix.install([
  {:httpoison, "~> 1.8"},
  {:json, "~> 1.4"},
  {:explorer, "~> 0.3"},
  {:kino, "~> 0.7"},
  {:kino_vega_lite, "~> 0.1.4"}
])
```

必要なライブラリをインストールします

- httpoison: REST API を呼び出す
- json: JSON をエンコード、デコードする
- explorer: データ分析
- kino: 実行結果を可視化する
- kino_vega_lite: データをグラフ化する

以下のコードを実行するとテキストエリアが表示されるので、 Tellus で作っておいたトークンを入力します

```elixir
# RESAS のAPIキーを入力する
api_key_input = Kino.Input.password("API_KEY")
```

ベースAPIを設定します

```elixir
base_url = "https://opendata.resas-portal.go.jp"
```

エイリアスを指定しておきます

```elixir
alias Explorer.DataFrame
alias Explorer.Series
alias VegaLite, as: Vl
```

## RESAS の認証

RESAS の認証は API キーをヘッダーに入れるだけです

認証ヘッダーを設定しておきます

```elixir
auth_header = {"X-API-KEY", Kino.Input.read(api_key_input)}
```

## 都道府県一覧の取得

都道府県一覧は以下のようにして呼び出します

```elixir
prefectures_url = "#{base_url}/api/v1/prefectures"

prefectures_df =
  prefectures_url
  |> HTTPoison.get!([auth_header])
  |> then(&JSON.decode!(&1.body))
  |> then(&DataFrame.new(&1["result"]))

prefectures_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 10.33.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/0c977656-7248-a015-cbec-388d96ddeab1.png)

大分県の都道府県コードを取得します

```elixir
oita_pref_code =
  prefectures_df
  |> DataFrame.filter_with(&Series.equal(&1["prefName"], "大分県"))
  |> DataFrame.pull("prefCode")
  |> Series.first()
```

![スクリーンショット 2022-11-22 10.35.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/a7353a0d-9b9d-2c35-0c5e-3dcedcbeb0cd.png)

## 産業一覧の取得

産業は大分類、中分類、小分類でそれぞれ分類されています

```elixir
# 大分類
industries_broad_df =
  "#{base_url}/api/v1/industries/broad"
  |> HTTPoison.get!([auth_header])
  |> then(&JSON.decode!(&1.body))
  |> then(&DataFrame.new(&1["result"]))

industries_broad_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 10.36.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/0903acc5-5cbf-5b06-b27d-169b3fec8918.png)

中分類は大分類コードでループして取得します

```elixir
# 中分類
get_industries_middle = fn sic_code ->
  "#{base_url}/api/v1/industries/middle?sicCode=#{sic_code}"
  |> HTTPoison.get!([auth_header])
  |> then(&JSON.decode!(&1.body))
  |> then(&DataFrame.new(&1["result"]))
end

industries_middle_df =
  industries_broad_df
  |> DataFrame.pull("sicCode")
  |> Series.to_list()
  |> Enum.map(&get_industries_middle.(&1))
  |> Enum.filter(&(DataFrame.n_rows(&1) > 0))
  |> DataFrame.concat_rows()

industries_middle_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 10.37.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/43afe333-313f-bd8c-0eb7-6115df1d2e87.png)

小分類は中分類コードでループして取得します

```elixir
# 小分類
get_industries_narrow = fn simc_code ->
  "#{base_url}/api/v1/industries/narrow?simcCode=#{simc_code}"
  |> HTTPoison.get!([auth_header])
  |> then(&JSON.decode!(&1.body))
  |> then(&DataFrame.new(&1["result"]))
end

industries_narrow_df =
  industries_middle_df
  |> DataFrame.pull("simcCode")
  |> Series.to_list()
  |> Enum.map(&get_industries_narrow.(&1))
  |> Enum.filter(&(DataFrame.n_rows(&1) > 0))
  |> DataFrame.concat_rows()

industries_narrow_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 10.38.57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/4708c27a-5f7c-c201-a9e8-e5be9ae67567.png)

## 地域別特化係数の取得

産業構造を数値化するものとして、特化係数を取得します

特化係数は、ある地域のある産業が、全国と比べてどれだけ特化しているか、を示す値です

1.0 を超えていれば、その地域はその産業に特化した地域である、と言えます

今回は付加価値額を基準に算出した特化係数を使用します

$$
特化係数（付加価値額） = \frac{域内における当該産業の付加価値額÷域内における全産業の付加価値額}{全国の当該産業の付加価値額÷全国の全産業の付加価値額}
$$

地域別特化係数は以下のようにして取得します

```elixir
get_industry_power = fn year, pref_code, area_type, disp_type, sic_code, simc_code ->
  query =
    "?year=#{year}" <>
      "&prefCode=#{pref_code}" <>
      "&areaType=#{area_type}" <>
      "&dispType=#{disp_type}" <>
      "&sicCode=#{sic_code}" <>
      "&simcCode=#{simc_code}"

  url = "#{base_url}/api/v1/industry/power/forArea#{query}"

  url
  |> HTTPoison.get!([auth_header])
  |> then(&JSON.decode!(&1.body))
  |> then(& &1["result"])
end
```

大分県の 2016 年時点、農業の地域別特化係数（付加価値額）を取得する場合は以下のように指定します

```elixir
industry_power = get_industry_power.("2016", oita_pref_code, "1", "1", "A", "01")
```

![スクリーンショット 2022-11-22 10.47.51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8212dea0-67ef-0027-a257-73c6c31738f8.png)

このとき、大分県内の各市町村の特化係数と同時に、都道府県別の特化係数も取得できます

今回はこちらの都道府県別データを使っていきます

![スクリーンショット 2022-11-22 10.48.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/69fa8d13-0fa3-6489-b944-9170fafbe63b.png)

都道府県別データをデータフレーム化します

```elixir
prefectures_industry_power_df =
  industry_power["prefectures"]
  |> DataFrame.new()

prefectures_industry_power_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 10.49.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/c10f9c0b-cd00-f580-17a1-0bbd1ac60073.png)

都道府県毎にどれだけ農業に特化しているか、のデータが一覧化できました

グラフ化するときに使う関数を定義しておきます

```elixir
get_values = fn df, col ->
  df
  |> DataFrame.pull(col)
  |> Series.to_list()
end
```

棒グラフを表示する関数です

```elixir
# 都道府県別棒グラフ
prefectures_bar = fn df, col, value_label ->
  x = get_values.(df, "prefName")
  y = get_values.(df, col)

  Vl.new(width: 800, height: 400, title: "都道府県別#{value_label}")
  |> Vl.data_from_values(x: x, y: y)
  |> Vl.mark(:bar)
  |> Vl.encode_field(
    :x,
    "x",
    type: :nominal,
    title: "都道府県"
  )
  |> Vl.encode_field(
    :y,
    "y",
    type: :quantitative,
    title: value_label
  )
end
```

鹿児島、青森、宮崎、北海道が農業に特化しています

![スクリーンショット 2022-11-22 11.04.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/cb0b2484-5b38-4dd8-6ab0-31736439bee8.png)

中分類単位で全産業の特化係数を取得しましょう

```elixir
prefectures_industry_power_df =
  # 中分類コードリストを取得
  industries_middle_df
  |> DataFrame.select(["sicCode", "simcCode"])
  |> DataFrame.to_rows()
  # 中分類毎に特化係数を取得
  |> Enum.map(fn industry ->
    get_industry_power.(
      "2016",
      oita_pref_code,
      "1",
      "1",
      industry["sicCode"],
      industry["simcCode"]
    )
    |> then(& &1["prefectures"])
    # 該当データがないものは除去する
    |> Enum.filter(&(&1 != nil))
    # データを整形する
    |> Enum.map(fn datum ->
      Map.merge(datum, %{"sicCode" => industry["sicCode"], "simcCode" => industry["simcCode"]})
    end)
  end)
  # 中分類毎にデータフレーム化する
  |> Enum.map(&DataFrame.new(&1))
  # 該当データがないものは除去する
  |> Enum.filter(&(DataFrame.n_rows(&1) > 0))
  # 全中分類のデータフレームを結合する
  |> DataFrame.concat_rows()
  # 大分類、中分類名を表示するために結合する
  |> DataFrame.join(industries_broad_df)
  |> DataFrame.join(industries_middle_df)

prefectures_industry_power_df
|> Kino.DataTable.new(sorting_enabled: true)
```

グラフタイトルに産業中分類を表示するための関数を用意します

```elixir
# 中分類名取得
get_simc_name = fn simc_code ->
  industries_middle_df
  |> DataFrame.filter_with(&Series.equal(&1["simcCode"], simc_code))
  |> DataFrame.pull("simcName")
  |> Series.first()
end
``

```elixir
# 都道府県別産業構造中分類棒グラフ
prefectures_industry_bar = fn simc_code ->
  label = get_simc_name.(simc_code)

  prefectures_industry_power_df
  |> DataFrame.filter_with(&Series.equal(&1["simcCode"], simc_code))
  |> prefectures_bar.("value", "#{label}特化係数")
end
```

中分類コード `81` の学校教育特化係数を表示してみます

```elixir
prefectures_industry_bar.("81")
```

京都は教育関連の産業に特化しているようです

![スクリーンショット 2022-11-22 11.06.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/2e217a0b-5c1b-389f-149d-23ea198eb4ef.png)

中分類コード `50` の各種商品卸売業特化係数を表示してみます

```elixir
prefectures_industry_bar.("50")
```

卸売は東京に一極集中しています

北海道はマイナスなので、卸売では2016年時点で赤字のようです

![スクリーンショット 2022-11-22 11.07.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/130f1aea-6b5f-bb59-85e6-92ca71f41bbe.png)

## ピボット

特化係数と他の値の相関を見るため、データを形を変えます

今１列になっている中分類コードを、中分類コード別の列に変換します

このような操作をピボットと言います

変換前

- prefCode
- prefName
- simcCode
- value

変換後

- prefCode
- prefName
- 01
- 02
- 03
...

`DataFrame.pivot_wider`　で、行数を減らして列数を増やします

```elixir
pivot_df =
  prefectures_industry_power_df
  |> DataFrame.select(["prefCode", "prefName", "simcCode", "value"])
  |> DataFrame.pivot_wider("simcCode", "value")

pivot_df
|> Kino.DataTable.new(sorting_enabled: true)
```

今まで都道府県数 * 中分類数（一部欠落） になっていた行数が、都道府県数（47）になりました

![スクリーンショット 2022-11-22 11.12.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/f13d15dd-6c23-4f5f-4956-441d164e68cd.png)

ただし、 `03` の列を見てみると、栃木県、群馬県の値が nil になっています

おそらくデータが取得できなかったものと思われます

このような値が存在しない部分を欠損値と言います

今後の分析のため欠損値を 0 で補完します

```elixir
# 補完対象の列
cols =
  prefectures_industry_power_df
  |> DataFrame.distinct(["simcCode"])
  |> DataFrame.pull("simcCode")
  |> Series.to_list()
```

![スクリーンショット 2022-11-22 11.16.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8f527c7f-6fd2-7676-f939-8ef413828427.png)

`Series.fill_missing` を各列に適用することで補完できます

```elixir
# 欠損値の補完
pivot_df =
  cols
  |> Enum.reduce(pivot_df, fn col, df ->
    DataFrame.mutate_with(df, &%{col => Series.fill_missing(&1[col], 0.0)})
  end)

pivot_df
|> Kino.DataTable.new()
```

![スクリーンショット 2022-11-22 11.16.48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/fd70af8d-0a30-9d50-6ad2-7b3715dbbf78.png)

## 散布図の表示

中分類毎の特化係数がそれくらい相関しているか、散布図で見てみます

ついでに回帰直線も表示しましょう

```elixir
scatter = fn df, x_col, y_col, size ->
  x = get_values.(df, x_col)
  y = get_values.(df, y_col)

  Vl.new(width: size, height: size)
  |> Vl.data_from_values(x: x, y: y)
  |> Vl.encode_field(:x, "x",
    type: :quantitative,
    scale: [domain: [Enum.min(x), Enum.max(x)]],
    title: x_col
  )
  |> Vl.encode_field(:y, "y",
    type: :quantitative,
    scale: [domain: [Enum.min(y), Enum.max(y)]],
    title: y_col
  )
  |> Vl.layers([
    Vl.new()
    |> Vl.mark(:point),
    Vl.new()
    |> Vl.mark(:line)
    |> Vl.transform(regression: "x", on: "y"),
  ])
end
```

農業と林業の散布図は以下のようになります

```elixir
scatter.(pivot_df, "01", "02", 300)
```

　点が回帰直線の周辺にあり、そこそこ相関がありそうなことが分かります

つまり、農業に特化した地域は林業にも特化している傾向がある、ということです

![スクリーンショット 2022-11-22 11.19.48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/7bf2ab19-7ae7-d94b-228f-e56edd607a6f.png)

農業と娯楽業だとこのようになります

```elixir
scatter.(pivot_df, "01", "80", 300)
```

あまり関係性はないか、負の相関、つまり農業に特化していれば娯楽業には特化していない、という傾向がありそうです

![スクリーンショット 2022-11-22 11.22.40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/fbae5b94-fe99-a210-2f40-1058bfbc9536.png)

## 一人当たり地方税の取得

都道府県の豊かさをどう示すか難しいですが、税収を比べてみましょう

地方税は地方公共団体が課税する税金で、以下のようなものが含まれます

- 住民税
- 事業税
- 固定資産税
- 地方消費税
- 自動車税

一人当たりの地方税を見れば、どれだけたくさん稼いだり資産を持っていたりするのか分かりそうです

一人当たり地方税を以下のようにして取得します

```elixir
get_local_tax = fn pref_code ->
  query =
    "?prefCode=#{pref_code}&cityCode=-"

  url = "#{base_url}/api/v1/municipality/taxes/perYear#{query}"

  url
  |> HTTPoison.get!([auth_header])
  |> then(&JSON.decode!(&1.body))
  |> then(& &1["result"]["data"])
end
```

全都道府県について取得します

```elixir
local_tax_df =
  prefectures_df
  |> DataFrame.pull("prefCode")
  |> Series.to_list()
  |> Enum.flat_map(fn pref_code ->
    pref_code
    |> get_local_tax.()
    |> Enum.map(&Map.merge(&1, %{"prefCode" => pref_code}))
  end)
  |> DataFrame.new()
  |> DataFrame.join(prefectures_df)

local_tax_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 11.33.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/59888ab4-bf41-8c83-54af-e3a880c86226.png)

大分県の一人当たり地方税推移を見てみましょう

```
target_df = DataFrame.filter_with(local_tax_df, &Series.equal(&1["prefCode"], oita_pref_code))
```

Smart Cells の Chart を使います

2009 年から 2014 年に落ち込み、 2015 年以降は上昇しています

![スクリーンショット 2022-11-22 11.38.23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/888f13a0-008d-1204-cb0d-8e2c71209430.png)

2019 年の全国の一人当たり地方税を見てみましょう

```elixir
local_tax_df
|> DataFrame.filter_with(&Series.equal(&1["year"], 2019))
|> prefectures_bar.("value", "一人当たり地方税")
```

東京が突出して高いのが分かります

![スクリーンショット 2022-11-22 11.40.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/ba911932-343d-2dce-626f-d7f0075f63a8.png)

## 人口構成の取得

産業構造と人口構成の関係も見てみましょう

```elixir
get_population_composition = fn pref_code ->
  query = "?prefCode=#{pref_code}&cityCode=-"

  url = "#{base_url}/api/v1/population/composition/perYear#{query}"

  url
  |> HTTPoison.get!([auth_header])
  |> then(&JSON.decode!(&1.body))
  |> then(& &1["result"]["data"])
end
```

```
population_composition_df =
  prefectures_df
  |> DataFrame.pull("prefCode")
  |> Series.to_list()
  |> Enum.flat_map(fn pref_code ->
    pref_code
    |> get_population_composition.()
    |> Enum.flat_map(fn data ->
      data["data"]
      |> Enum.map(fn datum ->
        Map.merge(datum, %{"label" => data["label"], "prefCode" => pref_code})
      end)
    end)
  end)
  |> DataFrame.new()
  |> DataFrame.join(prefectures_df)

population_composition_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 11.43.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/7e80e065-a8e5-a411-bcdb-419e2477a30c.png)

前回と同じではありますが、大分県の人口構成推移です

```elixir
target_df =
  population_composition_df
  |> DataFrame.filter_with(&Series.equal(&1["prefCode"], oita_pref_code))
```

![スクリーンショット 2022-11-22 11.44.44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/422a8113-ad81-6bd3-7c10-a2c1cd2136c2.png)


全都道府県の 2015 年時点の人口構成を見てみましょう

```elixir
target_df =
  population_composition_df
  |> DataFrame.filter_with(&Series.equal(&1["year"], 2015))
  |> DataFrame.filter_with(&Series.not_equal(&1["label"], "総人口"))
```

東京、神奈川、大阪の人口が多いですね

![スクリーンショット 2022-11-22 11.45.55.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/3acc147f-5c89-5b4d-a3ec-5560a7150535.png)


人口構成データを産業構造データと同じ都道府県毎のデータにするため 2015 年に絞ってピボットします

また、人口構成を数値化するため、老年率を計算します

```
population_df =
  population_composition_df
  |> DataFrame.filter_with(&Series.equal(&1["year"], 2015))
  |> DataFrame.mutate(value: &Series.cast(&1["value"], :float))
  |> DataFrame.select(["prefCode", "prefName", "label", "value"])
  |> DataFrame.pivot_wider("label", "value")
  |> DataFrame.mutate_with(fn df ->
    [
      老年率: Series.divide(df["老年人口"], df["総人口"])
    ]
  end)

population_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 11.49.34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/0ee00438-d0c9-afa9-5490-0ad25e473bbb.png)

都道府県別の老年率を見てみましょう

```elixir
population_df
|> prefectures_bar.("老年率", "老年率")
```

秋田が高く、沖縄が低いです

![スクリーンショット 2022-11-22 11.50.48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/458686cd-ea1f-eaa4-705a-5c8cd293273b.png)

## 相関係数の取得

明確に相関を表すため、相関係数を算出します

産業構造、地方税、人口構成のデータを結合します

```elixir
target_df =
  local_tax_df
  |> DataFrame.filter_with(&Series.equal(&1["year"], 2015))
  |> DataFrame.mutate(value: &Series.cast(&1["value"], :float))
  |> DataFrame.select(["prefCode", "value"])
  |> DataFrame.rename(["prefCode", "localTax"])
  |> DataFrame.join(population_df)
  |> DataFrame.join(pivot_df)

target_df
|> Kino.DataTable.new(sorting_enabled: true)
```

![スクリーンショット 2022-11-22 11.53.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8a2d06e1-49dd-ec59-246b-a175e773fede.png)

いくつか散布図を見てみましょう

農業と地方税

```elixir
scatter.(target_df, "localTax", "01", 300)
```

右下の東京が地方税で突出しているため、全体が見えにくいです

![スクリーンショット 2022-11-22 11.55.10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/d77cdda6-9b4f-92fe-5a8b-89a5098b6600.png)

東京を外れ値として除外して表示しましょう

```elixir
target_df
|> DataFrame.filter_with(&Series.not_equal(&1["prefCode"], 13))
|> scatter.("localTax", "01", 300)
```

見やすくなりましたが、点がバラバラなので相関はなさそうです

![スクリーンショット 2022-11-22 11.56.38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b823d166-5806-b4be-005a-7b23960f8293.png)

地方税と総人口で見てみます

```elixir
scatter.(target_df, "localTax", "総人口", 300)
```

これまた右上に東京が突出しています

![スクリーンショット 2022-11-22 11.57.55.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/35293535-8698-17b2-07a7-ddb4bc1bf8bf.png)

東京以外で見てみましょう

```elixir
target_df
|> DataFrame.filter_with(&Series.not_equal(&1["prefCode"], 13))
|> scatter.("localTax", "総人口", 300)
```

これもそんなに相関がなさそうです

![スクリーンショット 2022-11-22 11.59.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/465acb8a-3c5e-cf83-30f4-7d9dc03b3fb4.png)

以前やったのと同じように相関行列を計算します

https://qiita.com/RyoWakabayashi/items/a06ae639337726be2e9d#%E7%9B%B8%E9%96%A2%E8%A1%8C%E5%88%97

相関行列に含ませる列を指定します

```elixir
cols =
  prefectures_industry_power_df
  |> DataFrame.distinct(["simcCode"])
  |> DataFrame.pull("simcCode")
  |> Series.to_list()

cols = ["localTax", "総人口", "年少人口", "生産年齢人口", "老年人口", "老年率"] ++ cols
```

データを標準化します

```
standardize = fn df, column ->
  mean =
    df
    |> DataFrame.pull(column)
    |> Series.mean()

  std =
    df
    |> DataFrame.pull(column)
    |> Series.std()

  df
  |> DataFrame.mutate_with(fn in_df ->
    %{column => Series.subtract(in_df[column], mean)}
  end)
  |> DataFrame.mutate_with(fn in_df ->
    %{column => Series.divide(in_df[column], std)}
  end)
end
```

```elixir
standardized_df =
  cols
  |> Enum.reduce(target_df, fn col, df ->
    standardize.(df, col)
  end)

standardized_df
|> Kino.DataTable.new()
```

![スクリーンショット 2022-11-22 12.02.12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/4477a1a5-c5d9-a88b-7354-1cefbc79dede.png)

行列化して転置したものと内積を取ります

```elixir
df_to_tensor = fn df ->
  df
  |> DataFrame.names()
  |> Enum.map(fn col ->
    standardized_df
    |> DataFrame.pull(col)
    |> Series.to_tensor()
  end)
  |> Nx.concatenate()
  |> Nx.reshape({DataFrame.n_columns(df), DataFrame.n_rows(df)})
end

standardized_tensor =
  standardized_df
  |> DataFrame.select(cols)
  |> df_to_tensor.()
  |> Nx.transpose()
```

```elixir
covariance_tensor =
  standardized_tensor
  |> Nx.transpose()
  |> Nx.dot(standardized_tensor)
  |> Nx.divide(DataFrame.n_rows(standardized_df))
```

![スクリーンショット 2022-11-22 12.04.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/659aef2d-7216-e4d2-8389-2fe4ca1f99b5.png)

相関行列をデータフレームに変換します

```elixir
add_cols_label = fn list, cols_ ->
  [{"x", cols_} | list]
end

covariance_df =
  cols
  |> Stream.with_index()
  |> Enum.map(fn {col, index} ->
    {col, Nx.to_flat_list(covariance_tensor[index])}
  end)
  |> add_cols_label.(cols)
  |> DataFrame.new()

covariance_df
|> Kino.DataTable.new(keys: ["x" | cols])
```
![スクリーンショット 2022-11-22 12.05.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/227564a1-9f75-fb06-8bfa-46a008ec1239.png)

地方税との相関係数をヒートマップ化します

```elixir
local_tax_heatmap =
  cols
  |> Stream.with_index()
  |> Enum.map(fn {col_1, index_1} ->
    covariance =
      covariance_tensor[0][index_1]
      |> Nx.to_number()

    %{
      x: "localTax",
      y: col_1,
      covariance: covariance
    }
  end)
  |> List.flatten()
```

```elixir
Vl.new(width: 100, height: 1600)
|> Vl.data_from_values(local_tax_heatmap)
|> Vl.mark(:rect)
|> Vl.encode_field(:x, "x", type: :nominal)
|> Vl.encode_field(:y, "y", type: :nominal)
|> Vl.encode_field(
  :fill,
  "covariance",
  type: :quantitative,
  scale: [
    domain: [-1, 1],
    scheme: :blueorange
  ]
)
```

赤いところは正の相関が強く、青いところは負の相関が強いです

![visualization.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8c55aac7-2f6f-ef10-c670-297c6c987c53.png)

中分類コード `37` 通信業の相関は強そうです

散布図を見てみましょう

```elixir
# 通信業
target_df
|> scatter.("localTax", "37", 300)
```

通信業は東京一強なので、相関が強くなっているようです

![スクリーンショット 2022-11-22 12.10.29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/a07dea21-213a-b6f1-6e68-19e3ebd39b26.png)

```elixir
prefectures_industry_bar.("37")
```

![スクリーンショット 2022-11-22 12.12.26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/f85120a2-df97-48a0-4379-794885f88cb6.png)

逆に社会保険・社会福祉・介護事業は真っ青です

```elixir
# 社会保険・社会福祉・介護事業
target_df
|> scatter.("localTax", "85", 300)
```

![スクリーンショット 2022-11-22 12.14.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/bb8bd609-ef29-4ad1-f638-76290010f419.png)

```elixir
prefectures_industry_bar.("85")
```

![スクリーンショット 2022-11-22 12.15.42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/c9a8e5fc-79c4-6436-b558-a3f2f101b641.png)

お金を稼いでいる東京ですが、医療や福祉では全く稼げていないことが分かります

続いて老齢率の相関ヒートマップも見てみましょう

```elixir
elderly_heatmap =
  cols
  |> Stream.with_index()
  |> Enum.map(fn {col_1, index_1} ->
    covariance =
      covariance_tensor[5][index_1]
      |> Nx.to_number()

    %{
      x: "老年率",
      y: col_1,
      covariance: covariance
    }
  end)
  |> List.flatten()
```

```elixir
Vl.new(width: 100, height: 1600)
|> Vl.data_from_values(elderly_heatmap)
|> Vl.mark(:rect)
|> Vl.encode_field(:x, "x", type: :nominal)
|> Vl.encode_field(:y, "y", type: :nominal)
|> Vl.encode_field(
  :fill,
  "covariance",
  type: :quantitative,
  scale: [
    domain: [-1, 1],
    scheme: :blueorange
  ]
)
```

![visualization (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/eae3a46d-89c2-908d-3d22-de73db9c944b.png)

そんなに相関が強いものはなく、不動産取引で負の相関があるくらいでしょうか

```elixir
# 不動産取引業
target_df
|> scatter.("老年率", "68", 300)
```

![スクリーンショット 2022-11-22 12.19.01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/02338147-a4a4-baa9-b124-f8da9017e770.png)

```elixir
prefectures_industry_bar.("68")
```

![スクリーンショット 2022-11-22 12.20.04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/d46d39ef-5f1c-69a2-d9c9-55b2c19da0d8.png)

お年寄りが多いと、あまり引っ越しをしない、ということでしょうか

## まとめ

RESAS には他にもいろいろな API があるので、色々なデータの組み合わせを見ると面白そうです