中島さん卒論_backup の履歴(No.18)

• 履歴一覧
• 差分 を表示
• 現在との差分 を表示
• ソース を表示
• 中島さん卒論_backup へ行く。
  • 1 (2025-01-30 (木) 16:22:33)
  • 2 (2025-01-30 (木) 18:11:54)
  • 3 (2025-02-03 (月) 14:40:49)
  • 4 (2025-02-11 (火) 16:09:10)
  • 5 (2025-02-17 (月) 13:10:43)
  • 6 (2025-02-17 (月) 14:12:44)
  • 7 (2025-02-18 (火) 13:38:18)
  • 8 (2025-02-19 (水) 20:42:55)
  • 9 (2025-02-20 (木) 12:17:12)
  • 10 (2025-02-21 (金) 04:52:25)
  • 11 (2025-02-21 (金) 13:18:59)
  • 12 (2025-02-21 (金) 17:23:18)
  • 13 (2025-02-21 (金) 17:34:25)
  • 14 (2025-02-23 (日) 08:43:47)
  • 15 (2025-02-24 (月) 19:52:13)
  • 16 (2025-02-25 (火) 11:02:28)
  • 17 (2025-02-27 (木) 13:34:46)
  • 18 (2025-02-28 (金) 10:59:30)
  • 19 (2025-03-06 (木) 09:14:04)

---

技術資料

目次 

本研究の目的 

本研究の目的は，ヘドニック・アプローチにおける限界を克服し，価格形成におけるより 正確で信頼性の高い評価手法を提供することである．ヘドニックモデルでは，多重共線性， 欠落変数のバイアス，交互作用効果および非線形効果の取り扱いに問題がある[2]．これら は不動産市場や消費財の価格評価において誤差を生じさせる要因となる．本研究は，これ らの問題を解決するために，スパース推定法を活用した新たなアプローチを提案する． 従来のヘドニックモデルは，そのシンプルな線形構造と主要変数の関係を明確に示す点 で強力であったが，現実の市場では非線形な相互作用や高次の効果が絡むため，これらを 十分に捉えることができないことが多かった．また，異なる要因が相互作用を持つ場合，そ の効果を正確に評価するためには，より複雑なモデルが求められる．そこで本研究では， ElasticNetおよびその進化形であるAdaptiveElasticNet(AEN)を用い，従来の限界を克服 することを目指す． ElasticNetはL1およびL2正則化を組み合わせることで，多重共線性の問題に対応する と同時に，重要な変数の選択を行う能力を持つ．これにより，データ中で重要な特徴を捉 え，効果的な変数選択を実現できる．また，AENは，従来のElasticNetをさらに発展させ， データの特性に応じて正則化パラメータを動的に調整する．この適応的なアプローチは，複 雑なデータ構造に対して非常に効果的であり，非線形性や交互作用を持つデータにも適用 できる．これにより，不動産市場や消費財の価格形成における微細な変動や相互作用を捉 え，価格予測の精度を向上させることが可能となる． さらに，本研究では，従来の線形回帰モデルでは捉えきれなかった非線形性や交互作用 効果をスパース推定法を活用することで反映させることができる．非線形回帰や機械学習 技術と組み合わせることにより，より精緻な価格形成モデルが構築され，従来の方法では 見過ごされがちな要因を捉えることが可能となる．例えば，不動産市場における立地や面 積の相互作用，消費財におけるブランド力と機能性の複合的影響など，これらの複雑な関 係性を反映させることで，より現実的な市場分析が行えるようになる． また，こうした高度な統計手法や機械学習技術の活用により，従来のヘドニック・アプ ローチを大きく強化することができる．特に，従来のアプローチでは捉えきれなかった価 格形成における微細な要因を，スパース推定を用いることで明示化し，より高精度な価格 評価を可能にする．これにより，不動産市場や消費財市場における価格形成メカニズムの 理解が深まり，実務的な応用に貢献することが期待される．

使用するファイル 

本研究で使用するファイルはすべてGoogleDriveにある。
user:iie.lab.tpu.2324@gmail.com

「マイドライブ」→「学生」→「24_o3中島」→「保存期間5年」→「３）卒論（プログラム）」の中にある。
「引継ぎ用.zip」をダウンロードする。
ファイルの中身は以下のようになっている。 これを上から順に実行していく過程を以下の章に示す。

それぞれのコードを動かす際に初めに必要なファイル 前処理.ipynbを実行する前にデータセットの作成で作成された住宅地_犯罪_最短距離.csvをコピーして配置する。

次に変数選択_削除.ipynbを実行する前に前処理で作成された削除_変数作成後.csvをコピーして配置する。

hedo_削除.ipynbを実行する前に削除_変数作成後.csvと選ばれた変数_削除.csvとさらに前処理.ipynbで作成した外れ値削除.csvから取引の行われた不動産の住所と緯度と経度を手動で抜き出して位置情報_削除.csvとして保存する。この3つのcsvファイルを事前に配置する。

最後に寄与率.ipynbを実行する前に位置情報_削除.csv、削除_変数作成後.csv、選ばれた変数_削除.csvとhedo.ipynbで作成された選択された特徴量とデータ_削除.csvの4つを配置する。

コードの変更点はディレクトリを変更するだけで大丈夫.

pythonのバージョンは12.0.0を用いる.

使用したモジュールのバージョンは下になります。（すべては使用しておりません）

犯罪発生データの取得 

「get_crime_data.py」を実行すると、富山県警が公開している「犯罪発生マップ」（http://www.machi-info.jp/machikado/police_pref_toyama/index.jsp）から犯罪発生データを取得し、「crimes.csv」として保存されます。

建築年数の算出 

建築年の年号を西暦にするために以下の式を書いた列を挿入する。

   =IF(LEFT(M2,2)="令和", 2018 + MID(M2,3,LEN(M2)-3),
   IF(LEFT(M2,2)="平成", 1988 + MID(M2,3,LEN(M2)-3),
   IF(LEFT(M2,2)="昭和", 1925 + MID(M2,3,LEN(M2)-3), "エラー")))

次に取引時点から西暦を抜き出す。

   =LEFT(W2,4)

そして、取引時点の西暦から建築年の西暦を引くことで建築年数を算出する。

最後に残す列は建築年数のみ。

使用するデータ 

実際の取引データはhttps://www.reinfolib.mlit.go.jp/realEstatePrices/　このサイトからダウンロードする。最新版を用いる場合はこちらからダウンロードを。

データセットの作成\data\estats_mesh、データセットの作成\data\estats_blockのデータは下のurlからダウンロードできる。最新版が更新されればそちらを使用すればよい。

(https://www.e-stat.go.jp/gis/statmap-search?page=1&type=1&toukeiCode=00200521&toukeiYear=2020&aggregateUnit=H&serveyId=H002005112020&statsId=T001108&datum=2000&prefCode=16)

(https://www.e-stat.go.jp/gis/statmap-search?page=1&type=1&toukeiCode=00200521&toukeiYear=2020&aggregateUnit=H&serveyId=H002005112020&statsId=T001101&datum=2000&prefCode=16)
　この二つのurlからメッシュごとのデータ。

(https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200521&tstat=000001136464&cycle=0&tclass1=000001136472&tclass2=000001159889&cycle_facet=tclass1%3Acycle&tclass3val=0)
　このurlから小地域ごとのデータを。

データの作成と収集（統計データや施設データ） 

「hop_統計.ipynb」を実行すると、取引データのメッシュごとの統計データの作成が行われる. NAVITIMEからのスクレイピングや地図画像からのデータ作成を行い、メッシュごとのデータを最終的にdata_by_mesh.csvに格納する.

   df_torihiki = pd.read_csv(
   HERE / "df_housing_toyamacityメッシュ入り.csv",
   index_col = "住所",
   encoding = 'cp932')　

メッシュコードを文字列として格納　　

 meshes = df_torihiki["mesh_code"].astype(str).tolist()

作成する統計データは以下のようになる。

   import pandas as pd
   # ここで上で指定した説明変数を生成
   df_estats_block = pd.read_csv(
   DATA_DIR / "estats_block.csv",
   index_col = "KEY_CODE")
   df_estats_mesh = pd.read_csv(
   DATA_DIR / "estats_mesh.csv",
   index_col = "KEY_CODE")
   df_selected_block = pd.DataFrame()
   df_selected_mesh = pd.DataFrame()
   for line in stats_data_text.splitlines(): 
   formula, label = line.rsplit(", ", 2)  
   if len( formula.split(" ") ) == 1:        
       if formula in df_estats_block.columns:           
           df_selected_block.loc[:, label] = df_estats_block.loc[:, formula]          
       if formula in df_estats_mesh.columns:        
           df_selected_mesh.loc[:, label] = df_estats_mesh.loc[:, formula]        
   elif len( formula.split(" ") ) >= 2:
       contains_keys = re.findall("T.........", formula)  
       if contains_keys[0] in df_estats_block.columns:
           for key in contains_keys:                
               formula = formula[:formula.find(key)] + "df_estats_block.loc[:, '" + \
                   formula[formula.find(key) : formula.find(key) + len(key)] + \
                       "']" + formula[formula.find(key) + len(key):]                        
           exec(f"df_selected_block.loc[:, '{label}'] = {formula}")                
       elif contains_keys[0] in df_estats_mesh.columns:           
           for key in contains_keys:               
               formula = formula[:formula.find(key)] + "df_estats_mesh.loc[:, '" + \
                       formula[formula.find(key) : formula.find(key) + len(key)] + \
                           "']" + formula[formula.find(key) + len(key):]                            
           exec(f"df_selected_mesh.loc[:, '{label}'] = {formula}")            
   df_selected_mesh = df_selected_mesh.fillna(0)
   df_selected_block = df_selected_block.fillna(0)
   df_selected_mesh
   df_selected_block

   # 次は施設データ．
   # 統計データと同じで，カンマの左側の数字はNAVITIMEのカテゴリID
   # コンビニ一覧のURLは "https://www.navitime.co.jp/category/0201/" なので，カテゴ 
   リIDは「0201」
   building_data_text = '''\
   0201, コンビニ
   0802001, 駅
   0805, 駐車場
   0812, 駐輪場
   0501, 金融機関
   0603, 旅館
   0608, ホテル
   0202, スーパー
   0205, ショッピングモール
   0204, デパート
   0502002, 警察署／交番
   0504001, 小学校
   0504002, 中学校
   0504003, 高校
   0504004, 大学
   0504006, 保育園／幼稚園
   0101, レジャー施設'''
   # ここで，NAVITIMEからスクレイピング
   df_buildings = pd.DataFrame()
   def address2coords(address: str):   
   for trials in range(5):
           res = requests.get(
           urllib.parse.quote(
               string = f"https://msearch.gsi.go.jp/address-search/AddressSearch?q={address}",
               safe = "://?="
           )
       )        
       if res.status_code == 200:            
           break       
       else:           
           time.sleep(5 * trials)           
   else:       
       raise Exception("国土地理院のジオコーディングAPIに接続できません。")   
   result = res.json()   
   time.sleep(0.5)
   # サーバに負担をかけないように0.5秒待つ    
   if len(result) == 0:        
       return None        
   else:        
       return result[0]["geometry"]["coordinates"]
       # 複数ある場合でも，0番目を採用
   def get_building_info(category):    
   df_buildings = pd.DataFrame(columns = ["name", "address", "lng", "lat"])
   for i in range(1, sys.maxsize):
       # "while True"と同義       
       if i == 1:           
           url = f"https://www.navitime.co.jp/category/{category}/16201/"           
       else:            
           url = f"https://www.navitime.co.jp/category/{category}/16201/?page={i}"        
       for trials in range(5):           
           res = requests.get(url)           
           if res.status_code == 200:
            break            
           else:                
               time.sleep(5)           
       else:           
           break
           # NAVITIME から 5 回連続で取得できなかったら、終了               
       soup = BeautifulSoup(res.text, "html.parser")      
       buildings = []
       for element in soup.find_all(class_ = "spot-name-text"):         
           content = element.get_text()
           buildings += [{"name": content}]         
       count = 0
       for element in soup.find_all(class_ = "spot-detail-value-text"):           
           content = element.get_text()           
           matches = re.findall(r"(北海道|東京都|大阪府|京都府|神奈川県|和歌山県|鹿児島県|..県)(.{1,6}[市郡区町村])", content)
           # 完璧ではないが、これで住所かを判定            
           if len(matches) > 0:                
               address = content
               coords = address2coords(address)
               print(address, coords)               
               buildings[count]["address"] = address if coords == None:
                   buildings[count]["lng"] = None
                   buildings[count]["lat"] = None
               else:
                   buildings[count]["lng"] = coords[0]
                   buildings[count]["lat"] = coords[1]                
               count += 1                    
       for building in buildings:      
           df_buildings = add_dict_to_df(building, df_buildings)            
   return df_buildings

集めたデータを用いて不動産の取引が行われた場所の存在するメッシュごとのデータをdata_by_meshに格納する。

緯度経度からメッシュを格納 

「緯度経度.py」を実行すると、地理情報（緯度・経度）から日本独自のメッシュコードを計算し、住宅データなどのCSVファイルにそのコードを追加し、「df_housing_toyamacityメッシュ入り.csv」を生成します。

   import pandas as pd
   def latlon_to_mesh(lat, lon):
   """
   緯度経度から9桁のメッシュコード（2分の1地域メッシュ）を取得する。
   :param lat: 緯度（度）
   :param lon: 経度（度）
   :return: 9桁のメッシュコード（文字列）
   """
   # 1次メッシュ
   mesh1_lat = int(lat * 1.5)
   mesh1_lon = int(lon - 100)
   # 2次メッシュ
   mesh2_lat = int((lat * 1.5 - mesh1_lat) * 8)
   mesh2_lon = int((lon - 100 - mesh1_lon) * 8)
   # 3次メッシュ
   mesh3_lat = int(((lat * 1.5 - mesh1_lat) * 8 - mesh2_lat) * 10)
   mesh3_lon = int(((lon - 100 - mesh1_lon) * 8 - mesh2_lon) * 10)
   # 4次メッシュ（2分の1地域メッシュ）
   mesh4_lat = int((((lat * 1.5 - mesh1_lat) * 8 - mesh2_lat) * 10 - mesh3_lat) * 2)
   mesh4_lon = int((((lon - 100 - mesh1_lon) * 8 - mesh2_lon) * 10 - mesh3_lon) * 2)
   # メッシュコードの組み立て（9桁になるよう修正）
   mesh_code = f"{mesh1_lat:01}{mesh1_lon:01}{mesh2_lat:01}{mesh2_lon:01}{mesh3_lat:01}{mesh3_lon:01}{mesh4_lat:01}{mesh4_lon:01}"
   # 最後の2桁のみ取り出して、9桁にする
   mesh_code = mesh_code[:9]
   return mesh_code
   # CSVファイルの読み込み
   file_path = r"C:\Users

t011\Desktop\研究\富山取引のほう富山市のみ

   \df_housing_toyamacity.csv"
   df = pd.read_csv(file_path, encoding='utf-8')
   # メッシュコードの計算
   df["mesh_code"] = df.apply(lambda row: latlon_to_mesh(row["緯度"], row["経度"]), axis=1)
   cols = ["mesh_code"] + [col for col in df.columns if col != "mesh_code"]
   df = df[cols]
   # 結果のCSVを保存
   output_path = r"C:\Users

t011\Desktop\研究\富山取引のほう富山市のみ

   \df_housing_toyamacityメッシュ入り.csv"
   df.to_csv(output_path, index=False)
   print(f"処理完了: {output_path}")

「hop_メッシュ.ipynb」を実行すると、メッシュコードごとにデータを結合し「取引_統計データ付き.csv」に生成します。緯度・経度をもとに、500mメッシュ単位で犯罪データを分類できます。

    # メッシュGeoJsonを，分かりやすい形に直す
     _meshes = {}
     for mesh in meshes:
   meshcode = mesh["properties"]["code"]
   coords = mesh["geometry"]["coordinates"][0]
   north = coords[1][1]
   south = coords[0][1]
   east = coords[2][0]
   west = coords[0][0]
   center = [(east + west) / 2, (north + south) / 2]
   _meshes[meshcode] = {
       "north": north,
       "south": south,
       "east": east,
       "west": west,
       "center": center
   }
   meshes = _meshes
   # 犯罪データの座標から，どのグリッドセルに含まれるかを求め
   # 「mesh_code」列に格納する
   # 読み込んだメッシュGeoJsonに含まれない犯罪データは除外
   df_temp = df_crimes[0:0].copy()
   df_temp["mesh_code"] = None
   for index, row in df_crimes.iterrows():
   for meshcode, coords in meshes.items():
       if ((coords["south"] <= row["latitude"] < coords["north"]) and
           (coords["west"] <= row["longitude"] < coords["east"])):
           if not pd.isnull(row["datetime"]):
               df_temp.loc[index] = pd.concat([
                   row,
                   pd.Series({
                       "mesh_code": meshcode
                   })
               ])
           break
   print(f"{df_crimes.index.get_loc(index)} / {len(df_crimes)}")
   df_crimes = df_temp.copy()
   df_crimes["datetime"] = pd.to_datetime(df_crimes["datetime"], format = "%Y/%m/%d %H:%M")
   # CSVファイルの読み込み
   df_crimes = pd.read_csv(HERE / "犯罪_メッシュ入り.csv")
   #メッシュごとに件数をカウント
   mesh_counts = df_crimes['mesh_code'].value_counts().reset_index()
   mesh_counts.columns = ['mesh_code', '件数']
   # 結果を表示
   print(mesh_counts)
   # 結果をCSVファイルに保存する場合
   mesh_counts.to_csv(HERE / "犯罪_メッシュ件数.csv", index=False, encoding='utf-8-sig')

   # 件数データ
   df_number = pd.read_csv(
   HERE / "犯罪_メッシュ件数.csv",
   index_col = "mesh_code")

   # 取引データ
   df_torihiki = pd.read_csv(
   HERE / "df_housing_toyamacityメッシュ入り.csv",
   index_col = "mesh_code",
   encoding = 'cp932')

   # メッシュコードをキーとして結合
   df_merged = pd.merge(df_torihiki, df_number, on='mesh_code', how='left', 
   suffixes=('', '_number'))
   print(df_merged.head())
   # 結果をCSVファイルに保存する場合
   df_merged.to_csv(HERE / "取引_件数付き.csv", encoding='utf-8-sig')

   df_mesh = pd.read_csv(
   HERE / "data_by_mesh.csv",
   index_col = "mesh_code",
   encoding = "cp932")

   df_torihiki2 = pd.read_csv(
   HERE / "取引_件数付き.csv",
   index_col = "mesh_code")

   # メッシュコードでデータを結合
   merged_df = pd.merge(df_torihiki2, df_mesh, on='mesh_code', how='left')
   print(merged_df.head())
   merged_df.to_csv(HERE / "取引_統計データ付き.csv", index=False, 
   encoding='utf-8-sig')

犯罪発生率(件/250000)はhop_メッシュ.ipynbを行った後にcsvファイル上の新しい列(犯罪発生率(件/250000)を追加して行う。
計算方法は=件数/250000で行う。
件数列の空欄は0に統一する。
取引_統計データ付き.csvに手動でdf_housing_toyamacityメッシュ入り.csvからmesh_code列を一番左の列に追加する。

「distance.ipynb」を実行すると、住宅地データ (取引_統計データ付き.csv) と犯罪データ (犯罪_メッシュ入り.csv) を読み込み、各住宅地から最寄りの犯罪発生地点までの最短距離を計算するために、geopy.distanceモジュールのgeodesic関数を使用します。この関数は、2点間の地球上での距離（メートル単位）を計算します。契約された土地から最短の距離で犯罪が起きた場所までの最短距離を「取引_統計データ付き.csv」に追加した「住宅地_犯罪_最短距離.csv」を生成します。

   import pandas as pd
   from geopy.distance import geodesic
   # ファイルを読み込む
   df_housing = pd.read_csv("C:/Users/nt011/Desktop/統計正しいほう/データセット 
  の作成/取引_統計データ付き.csv")
   df_crimes = pd.read_csv("C:/Users/nt011/Desktop/統計正しいほう/データセット 
  の作成/犯罪_メッシュ入り.csv")
   print(df_housing.columns)
   # 1. 住宅地データの緯度経度情報を抽出
   df_housing['latitude'] = df_housing['緯度']
   df_housing['longitude'] = df_housing['経度']
   # 2. 犯罪データの緯度経度情報を抽出
   df_crimes['latitude'] = df_crimes['latitude']
   df_crimes['longitude'] = df_crimes['longitude']
   # 3.最短距離を計算する関数
   def calculate_min_distance(lat1, lon1, df_target, lat_col, lon_col):
   distances = df_target.apply(lambda row: geodesic((lat1, lon1), 
   (row[lat_col], row[lon_col])).meters, axis=1)
   return distances.min()
   # 4. 各住宅地から犯罪発生場所までの最短距離を計算
   df_housing['min_distance_to_crime'] = df_housing.apply(lambda row: 
   calculate_min_distance(row['latitude'], row['longitude'], df_crimes, 'latitude', 'longitude'), axis=1)
   # 距離が0になる住宅地の行を表示
   df_housing[df_housing['min_distance_to_crime'] == 0]
   # 必要なカラムを選択して新しいDataFrameに保存
   df_result = df_housing[['mesh_code', '住所', '経度', '緯度', '最寄駅：名称', '最寄駅：距離（分）', '取引価格（総額）', '面積（㎡）', '取引価格（㎡単価）', 
                       '土地の形状', '間口', '延床面積（㎡）','建築年数', '建物の構造', 
                       '用途', '前面道路：方位', '前面道路：種類', '前面道路：幅員（ｍ）', '都市計画', '建ぺい率（％）', 
                       '容積率（％）', '取引時点', '件数', '犯罪発生率(件/250000)', '人口', '18歳未満人口割合', 
                       '65歳以上人口割合', '外国人人口割合', '世帯数', '単身世帯割合', '核家族世帯割合', '正規労働者割合', 
                       '非正規労働者割合', '最終学歴が中学以下の人口割合', '最終学歴が高校の人口割合', '最終学歴が大学以上の人口割合', 
                       '居住年数5年未満の人口割合', '居住年数が20年以上の人口割合', '一戸建て割合', 'アパート・低中層マンション割合', 
                       '高層マンション割合', 'コンビニ', '駅', '駐車場', '駐輪場', '金融機関', '旅館', 'ホテル', 'スーパー', 'ショッピングモール', 'デパート', '警察署／交番', '小学校', '中学校', '高校', '大学',
                       '保育園／幼稚園', 'レジャー施設','道路の面積比率', '建物の面積比率', '水の面積比率', '空き地の面積比率', 
                       'ノード数（道路）', 'エッジ数（道路）', '密度（道路）', '平均次数（道路）', 'min_distance_to_crime']]
   # 結果の保存
   df_result.to_csv('C:/Users/nt011/Desktop/統計正しいほう/データセットの作成/ 
住宅地_犯罪_最短距離.csv', index=False, encoding='utf-8-sig')

データセットの前処理 

作成したデータセットを回帰分析がかけられる形に整える。

「前処理.ipynb」を実行し、「住宅地_犯罪_最短距離.csv」から「建築年数」と「最寄駅：距離（分）」の列に関する欠損値を含む行を削除した「住宅地_犯罪_最短距離_欠損値削除後.csv」を生成します。

   import pandas as pd
   # CSVファイルのパス
   file_path = r"C:\Users
   t011\Desktop\統計正しいほう\前処理\住宅地_犯罪_最短距離.csv"
   # CSVファイルを読み込む
   df = pd.read_csv(file_path)
   # 欠損値を含む行を削除
   df = df.dropna()
   # 「建築年数」列に `-1` または `"VALUE"` を含む行を削除
   df = df[~df["建築年数"].astype(str).isin(["-1", "#VALUE!"])]
   # 「最寄駅：距離（分）」の値を変換
   replace_dict = {
   "30分~60分": 45,
   "1H~1H30": 75,
   "1H30~2H": 105
   }
   df["最寄駅：距離（分）"] = df["最寄駅：距離（分）"].replace(replace_dict)
   # クリーンなデータを新しいCSVファイルとして保存
   output_path = r"C:\Users
   t011\Desktop\統計正しいほう\前処理\住宅地_犯罪_最短距離_欠損値削除後.csv"
   df.to_csv(output_path, index=False, encoding='utf-8-sig')
   print("データ処理完了: 欠損値と指定された異常値を削除し、駅距離の値を変換しました。")
   print(f"新しいCSVファイルが {output_path} に保存されました。")

生成した「住宅地_犯罪_最短距離_欠損値削除後.csv」から連続変数の外れ値を検出するカラムリストを定義し、データ型を数値に変換したのち検出された列を削除したリスト「外れ値削除.csv」を生成します。ダミー変数の処理を行う前にcsvファイル内の建物の構造の英語の部分を変換する必要がある。

◆不動産取引価格
https://www.reinfolib.mlit.go.jp/realEstatePrices/

このurlの制度と用語のページを参考にして変換を行う。

   import pandas as pd
   # CSVファイルのパス
   file_path = "C:/Users/nt011/Desktop/統計正しいほう/前処理/住宅地_犯罪_最短距 
   離_欠損値削除後.csv"
   # CSVファイルを読み込む
   df = pd.read_csv(file_path)
   # 外れ値を検出する連続変数のカラムリスト
   continuous_columns = [
   "最寄駅：距離（分）", "取引価格（総額）", "面積（㎡）", "取引価格（㎡単価）",
   "間口", "延床面積（㎡）", "建築年数", "前面道路：幅員（ｍ）",
   "建ぺい率（％）", "容積率（％）", "犯罪発生率(件/250000)", "人口",
   "18歳未満人口割合", "65歳以上人口割合", "外国人人口割合", "世帯数",
   "単身世帯割合", "核家族世帯割合", "正規労働者割合", "非正規労働者割合",
   "最終学歴が中学以下の人口割合", "最終学歴が高校の人口割合", "最終学歴が大学以上の人口割合",
   "居住年数5年未満の人口割合", "居住年数が20年以上の人口割合", "一戸建て割合",
   "アパート・低中層マンション割合", "高層マンション割合", "道路の面積比率",
   "建物の面積比率", "水の面積比率", "空き地の面積比率", "ノード数（道路）",
   "エッジ数（道路）", "密度（道路）", "平均次数（道路）", "min_distance_to_crime"
   ]
   # データ型を数値に変換（エラーが出たらNaNにする）
   for column in continuous_columns:
   if column in df.columns:
       df[column] = pd.to_numeric(df[column], errors='coerce')
   # 外れ値を持つ行のインデックスを格納するリスト
   outlier_indices = set()
   # 外れ値の検出と削除
   for column in continuous_columns:
   if column in df.columns:
       # 欠損値を除外
       df_clean = df[column].dropna()
       # IQR（四分位範囲）の計算
       Q1 = df_clean.quantile(0.25)
       Q3 = df_clean.quantile(0.75)
       IQR = Q3 - Q1
       # 外れ値の閾値を設定
       lower_bound = Q1 - 1.5 * IQR
       upper_bound = Q3 + 1.5 * IQR
       # 外れ値のインデックスを取得
       outliers = df[(df[column] < lower_bound) | (df[column] > upper_bound)]
       outlier_indices.update(outliers.index)
   # 外れ値の行を削除
   df_cleaned = df.drop(index=outlier_indices)
   # 結果を新しいCSVに保存
   output_path = "C:/Users/nt011/Desktop/統計正しいほう/前処理/外れ値削除.csv"
   df_cleaned.to_csv(output_path, encoding="utf-8-sig", index=False)
   print(f"外れ値を削除したデータを {output_path} に保存しました。")
   print(f"削除した行数: {len(outlier_indices)}")
   print(f"最終データの行数: {df_cleaned.shape[0]}")

それぞれのカラムの要素を必要に応じてダミー変数化した後、データフレームを保存した「削除_ダミー後.csv」を生成します。

   import pandas as pd
   # CSVファイルを読み込み
   df_data = pd.read_csv("C:/Users/nt011/Desktop/統計正しいほう/前処理/外れ値削除.csv")
   # カラム名を取得
   columns = df_data.columns.tolist()
   # 建物の構造列に「RC」がある場合、「鉄筋コンクリート造」に変換
   df_data['建物の構造'] = df_data['建物の構造'].replace('RC', '鉄筋コンクリート造')
   # 最寄り駅名を二値化（列名を「最寄り駅名_富山」に変更）
   index = columns.index("最寄駅：名称")  # 変換前の位置を取得
   df_data.insert(index + 1, '最寄り駅名_富山', df_data['最寄駅：名称'].apply(lambda x: 1 if x == '富山' else 0))
   df_data = df_data.drop(columns=['最寄駅：名称'])  # 元の列を削除
   # 都市計画を二値化（列名を「都市計画_市街化調整区域」に変更）
   index = columns.index("都市計画")
   df_data.insert(index + 1, '都市計画_市街化調整区域', df_data['都市計画'].apply(lambda x: 1 if x == '市街化調整区域' else 0))
   df_data = df_data.drop(columns=['都市計画'])
   # 土地の形状を二値化（列名を「土地の形状_不整形」に変更）
   index = columns.index("土地の形状")
   df_data.insert(index + 1, '土地の形状_不整形', df_data['土地の形状'].apply(lambda x: 1 if x == '不整形' else 0))
   df_data = df_data.drop(columns=['土地の形状'])
   # 用途を二値化（列名を「用途_住宅」に変更）
   index = columns.index("用途")
   df_data.insert(index + 1, '用途_住宅', df_data['用途'].apply(lambda x: 1 if x == '住宅' else 0))
   df_data = df_data.drop(columns=['用途'])
   # 建物の構造を二値化（列名を「建物の構造_木造」に変更）
   index = columns.index("建物の構造")
   df_data.insert(index + 1, '建物の構造_木造', df_data['建物の構造'].apply(lambda x: 1 if x == '木造' else 0))
   df_data = df_data.drop(columns=['建物の構造'])
   # 「前面道路：方位」のダミー変数化（基準値「北西」を0に設定）
   index = columns.index("前面道路：方位")  # 変換前の位置を取得
   dummies = pd.get_dummies(df_data['前面道路：方位'], prefix='前面道路：方位', dtype=int)
   # 基準値は削除
   baseline_value = '北西'
   if f'前面道路：方位_{baseline_value}' in dummies.columns:
   dummies = dummies.drop(columns=[f'前面道路：方位_{baseline_value}'])
   # 挿入（前面道路：方位のダミー変数）
   ordered_dummies = ['前面道路：方位_西', '前面道路：方位_東', '前面道路：方位_南西', '前面道路：方位_南東', 
                  '前面道路：方位_南', '前面道路：方位_北東', '前面道路：方位_北']
   # 順番に従って挿入
   for col in ordered_dummies:
   if col in dummies.columns:
       df_data.insert(index + 1, col, dummies[col])
   # 「前面道路：種類」のダミー変数化
   road_types = ['市道', '私道', '県道']
   # 順番を指定してダミー変数を挿入
   for road in road_types:
   df_data.insert(index + 1, f'前面道路：種類_{road}', (df_data['前面道路：種類'] == road).astype(int))
   # 元の「前面道路：方位」および「前面道路：種類」列を削除
   df_data = df_data.drop(columns=['前面道路：方位', '前面道路：種類'])
   # 確認用にダミー変数化後のデータフレームを表示
   print(df_data.head())
   # ダミー変数化後のデータを保存
   df_data.to_csv("C:/Users/nt011/Desktop/統計正しいほう/前処理/削除_ダミー後.csv", 
   index=False, encoding='utf-8-sig')
   print("ダミー変数変換完了！データが保存されました。")

連続変数とダミー変数を手動で分類し、連続変数のみ二次項の作成をする。 交互作用項の作成したデータフレームを保存した「削除_変数作成後.csv」を生成します。

   import pandas as pd
   import numpy as np
   from itertools import combinations
   # CSVファイルのパスを指定
   file_path = "C:/Users/nt011/Desktop/統計正しいほう/前処理/削除_ダミー後.csv"
   # UTF-8で読み込む
   df = pd.read_csv(file_path)
   # データの先頭を表示
   print(df.head())
   # 連続変数とダミー変数を手動で分類
   continuous_vars = [
   '最寄駅：距離（分）', '面積（㎡）', '間口', '延床面積（㎡）', '建築年数', '前面道路：幅員（ｍ）',
   '建ぺい率（％）', '容積率（％）', '犯罪発生率(件/250000)', '人口', '18歳未満人口割合',
   '65歳以上人口割合', '外国人人口割合', '世帯数', '単身世帯割合', '核家族世帯割合', 
   '正規労働者割合', '非正規労働者割合', '最終学歴が中学以下の人口割合', '最終学歴が高校の人口割合',
   '最終学歴が大学以上の人口割合', '居住年数5年未満の人口割合', '居住年数が20年以上の人口割合',
   '一戸建て割合', 'アパート・低中層マンション割合', 'コンビニ', '駅', '駐車場', '駐輪場', 
   '金融機関', '旅館', 'ホテル', 'スーパー', 'ショッピングモール', 
   'デパート', '警察署／交番', '小学校', '中学校', '高校', '大学', '保育園／幼稚園', 
   'レジャー施設','道路の面積比率', '建物の面積比率', '水の面積比率', '空き地の面積比率', 
   'ノード数（道路）', 'エッジ数（道路）', '密度（道路）', '平均次数（道路）', 'min_distance_to_crime'
   ]
   # ダミー変数をリストとして定義
   dummy_vars = [
   '最寄り駅名_富山','土地の形状_不整形', '建物の構造_木造', 
   '用途_住宅', '前面道路：方位_北', '前面道路：方位_北東', 
   '前面道路：方位_南', '前面道路：方位_南東', 
   '前面道路：方位_南西', '前面道路：方位_東', '前面道路：方位_西', 
   '前面道路：種類_市道', '前面道路：種類_県道', '前面道路：種類_私道', 
   '都市計画_市街化調整区域'
   ]
   # 二次項の作成（連続変数のみ）
   for var in continuous_vars:
   df[f'{var}^2'] = df[var] ** 2
   # 交互作用項の作成
   # 連続変数 × 連続変数
   for var1, var2 in combinations(continuous_vars, 2):
   df[f'{var1}×{var2}'] = df[var1] * df[var2]
   # ダミー変数 × 連続変数
   for dummy_var in dummy_vars:
   for cont_var in continuous_vars:
       df[f'{dummy_var}×{cont_var}'] = df[dummy_var] * df[cont_var]
   # ダミー変数 × ダミー変数
   for var1, var2 in combinations(dummy_vars, 2):
   df[f'{var1}×{var2}'] = df[var1] * df[var2]
   # CSVに保存
   output_path = "C:/Users/nt011/Desktop/統計正しいほう/前処理/削除_変数作成後.csv"
   df.to_csv(output_path, index=False, encoding='utf-8-sig')
   print(f"処理が完了しました。生成されたデータは {output_path} に保存されました。")

変数選択 

Elastic Netとは 

Elastic Netは、L1 正則化（Lasso）と L2 正則化（Ridge）の組み合わせにより、多重共線性があるデータに対して安定した回帰分析を行う手法です。

Lasso（L1正則化）：L1正則化は、回帰モデルの係数の一部をゼロにすることによって、特徴量選択を行います。これにより、重要な特徴量だけを残すことができますが、多重共線性が高い場合には不安定になることがあります。

Ridge（L2正則化）：L2正則化は、回帰モデルの係数を小さくすることで、過学習を防ぎます。L2正則化は係数をゼロにすることはなく、むしろ係数の大きさを均等に抑える傾向がありますが、特徴量選択の機能はありません。

StandardScaler() を用いて X（説明変数）と y（目的変数）を標準化 しているため、変数のスケールによる影響を排除し、モデルの安定性を向上させています。

Adaptive Elastic Netを用いて土地の価格に影響する重要な特徴量の係数を取得し、影響の大きい特徴量を特定します。

手順 

「変数選択__削除.ipynb」を実行すると、重要な変数のみを抽出し、「選ばれた変数.csv」に保存します。 実行を行う前に削除_変数作成後.csvのmesh_code、住所、経度、緯度と取引価格（総額）、取引時点、件数の列を削除する。そして、取引価格（㎡単価）を一番左の列に移す。

   import numbers
   import warnings
   import cvxpy
   import numpy as np
   from asgl import ASGL
   from sklearn.base import MultiOutputMixin
   from sklearn.base import RegressorMixin
   from sklearn.exceptions import ConvergenceWarning
   from sklearn.linear_model import ElasticNet
   from sklearn.linear_model._coordinate_descent import _alpha_grid
   from sklearn.utils import check_X_y
   from sklearn.utils.validation import check_is_fitted
   class AdaptiveElasticNet(ASGL, ElasticNet, MultiOutputMixin, RegressorMixin):
   """
   Objective function and parameters as described with modifications
   to allow alpha1, alpha2, l1_ratio1, and l1_ratio2 customization.
   """
   def __init__(
       self,
       alpha1=0.021544346900318846,   # First-stage ElasticNet alpha
       alpha2=0.0009443498043343188,  # Second-stage ElasticNet alpha
       *,
       l1_ratio1=0.8889,  # First-stage ElasticNet L1 ratio
       l1_ratio2=0.778,  # Second-stage ElasticNet L1 ratio
       gamma=0.5,  # Weight adjustment exponent
       fit_intercept=True,
       precompute=False,
       max_iter=10000,
       copy_X=True,
       solver=None,
       tol=None,
       positive=False,
       positive_tol=1e-3,
       random_state=None,
       eps_coef=1e-6,
       verbose=True
   ):
       params_asgl = dict(model="lm", penalization="asgl")
       if solver is not None:
           params_asgl["solver"] = solver
       if tol is not None:
           params_asgl["tol"] = tol
       super().__init__(**params_asgl)
       self.alpha1 = alpha1
       self.alpha2 = alpha2
       self.l1_ratio1 = l1_ratio1
       self.l1_ratio2 = l1_ratio2
       self.gamma = gamma
       self.fit_intercept = fit_intercept
       self.max_iter = max_iter
       self.precompute = precompute
       self.copy_X = copy_X
       self.positive = positive
       self.positive_tol = positive_tol
       self.random_state = random_state
       self.eps_coef = eps_coef
       self.verbose = verbose
       if not self.fit_intercept:
           raise NotImplementedError
   def fit(self, X, y, check_input=True):
       if check_input:
           X_copied = self.copy_X and self.fit_intercept
           X, y = self._validate_data(
               X,
               y,
               accept_sparse="csc",
               order="F",
               dtype=[np.float64, np.float32],
               copy=X_copied,
               multi_output=True,
               y_numeric=True,
           )
       # 第一段階の ElasticNet 実行
       enet_model_1 = self.elastic_net(self.l1_ratio1, alpha=self.alpha1)
       enet_model_1.fit(X, y)
       enet_coef = enet_model_1.coef_
       # 重みの計算
       weights = 1.0 / (np.maximum(np.abs(enet_coef), self.eps_coef) ** self.gamma)
       # 第二段階の最適化
       self.coef_, self.intercept_ = self._optimize_second_stage(X, y, weights)
       # モデル属性を格納
       self.enet_coef_ = enet_coef
       self.weights_ = weights
       return self
   def predict(self, X):
       check_is_fitted(self, ["coef_", "intercept_"])
       return super(ElasticNet, self).predict(X)
   def elastic_net(self, l1_ratio, **params):
       """
       Create an ElasticNet model with the specified parameters.
       """
       elastic_net = ElasticNet(l1_ratio=l1_ratio)
       for k, v in self.get_params().items():
           try:
               elastic_net = elastic_net.set_params(**{k: v})
           except ValueError:
               pass  # Ignore parameters not supported by ElasticNet
       elastic_net = elastic_net.set_params(**params)
       return elastic_net
   def _optimize_second_stage(self, X, y, weights):
       """
       Perform second-stage optimization with adaptive weights.
       Returns
       -------
       coef : np.array, shape (n_features,)
       intercept : float
       """
       n_samples, n_features = X.shape
       beta_variables = [cvxpy.Variable(n_features)]
       model_prediction = 0.0
       if self.fit_intercept:
           beta_variables = [cvxpy.Variable(1)] + beta_variables
           ones = cvxpy.Constant(np.ones((n_samples, 1)))
           model_prediction += ones @ beta_variables[0]
       # モデル予測
       model_prediction += X @ beta_variables[1]
       error = cvxpy.sum_squares(y - model_prediction) / (2 * n_samples)
       l1_coefs = self.alpha2 * self.l1_ratio2
       # 第二段階の正則化項
       l1_penalty = cvxpy.Constant(l1_coefs * weights) @ cvxpy.abs(
           beta_variables[1]
       )
       l2_penalty = (
           cvxpy.Constant(self.alpha1 * (1 - self.l1_ratio1))
           * cvxpy.sum_squares(beta_variables[1])
       )
       constraints = [b >= 0 for b in beta_variables] if self.positive else []
       # 最適化問題の定義
       problem = cvxpy.Problem(
           cvxpy.Minimize(error + l1_penalty + l2_penalty), 
       constraints=constraints
       )
       problem.solve(solver="OSQP", max_iter=self.max_iter)
       if problem.status != "optimal":
           raise ConvergenceWarning(
               f"Solver did not reach optimum (Status: {problem.status})"
           )
       beta_sol = np.concatenate([b.value for b in beta_variables], axis=0)
       beta_sol[np.abs(beta_sol) < self.tol] = 0
       intercept, coef = beta_sol[0], beta_sol[1:]
       coef = np.maximum(coef, 0) if self.positive else coef
       return coef, intercept
   from sklearn.linear_model import ElasticNetCV
   from sklearn.model_selection import GridSearchCV
   from sklearn.datasets import make_regression
   import numpy as np
   from sklearn.model_selection import train_test_split
   import pandas as pd
   from sklearn.preprocessing import StandardScaler
   # CSVファイルのパスを指定してください
   file_path = "C:/Users/nt011/Desktop/統計正しいほう/変数選択/削除_変数作成後.csv"
   # CSVファイルの読み込み
   df = pd.read_csv(file_path)  # 文字コードが異なる場合は、'utf-8' を他のエンコーディングに変更してください
   # 特徴量とターゲットの分割
   X = df.drop(columns=['取引価格（㎡単価）'])  # 取引価格（㎡単価）をyとして分離
   y = df['取引価格（㎡単価）']  # ターゲット変数
   # 🔥 標準化の実施
   scaler_X = StandardScaler()
   scaler_y = StandardScaler()
   X_scaled = scaler_X.fit_transform(X)  # Xの標準化
   y_scaled = scaler_y.fit_transform(y.values.reshape(-1, 1)).ravel()  # yの標準化
   # データ分割
   X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_scaled, 
   test_size=0.2, random_state=42)
   # 訓練データの一部を使用
   X_sample, _, y_sample, _ = train_test_split(X_train, y_train, test_size=0.8, 
   random_state=42)
   # 第一段階: ElasticNetCVを使用した最適パラメータの導出
   enet_cv = ElasticNetCV(
   l1_ratio=np.linspace(0.0001, 1, 25),  # l1_ratioの候補
   alphas=np.logspace(-5, 0, 25),  # alphaの候補
   cv=5,  # 交差検証の分割数
   random_state=42,
   n_jobs=-1,
   verbose=1
    )
   enet_cv.fit(X_train, y_train)
   # 第一段階の最適パラメータと係数を取得
   alpha1_opt = enet_cv.alpha_
   l1_ratio1_opt = enet_cv.l1_ratio_
   enet_coef = enet_cv.coef_
   print(f"第一段階の最適パラメータ: alpha1={alpha1_opt}, l1_ratio1= 
   {l1_ratio1_opt}")

   # ベストパラメータをalpha2とl1_ratio2でそれぞれ格納
   best_alpha2 = best_params['alpha2']
   best_l1_ratio2 = best_params['l1_ratio2']
   print(best_alpha2)
   print(best_l1_ratio2)

ASGLとElasticNetを組み合わせたAdaptiveElasticNetクラスを定義し、CSVから読み込んだデータを標準化後にモデルをフィッティングして、非ゼロの係数を持つ重要な特徴量を抽出し、その結果をCSVに保存する一連の処理を実行します。

   import numbers
   import warnings
   import cvxpy
   import numpy as np
   from asgl import ASGL
   from sklearn.base import MultiOutputMixin
   from sklearn.base import RegressorMixin
   from sklearn.exceptions import ConvergenceWarning
   from sklearn.linear_model import ElasticNet
   from sklearn.linear_model._coordinate_descent import _alpha_grid
   from sklearn.utils import check_X_y
   from sklearn.utils.validation import check_is_fitted
   class AdaptiveElasticNet(ASGL, ElasticNet, MultiOutputMixin, RegressorMixin):
   """
   Objective function and parameters as described with modifications
   to allow alpha1, alpha2, l1_ratio1, and l1_ratio2 customization.
   """
   def __init__(
       self,
       alpha1=0.01333521432163324,    # 第一段階で得られたalpha1を指定
       alpha2=0.034807005884284134,  # Second-stage ElasticNet alpha
       *,
       l1_ratio1=1,  # First-stage ElasticNet L1 ratio
       l1_ratio2=0.1667499999999998,  # Second-stage ElasticNet L1 ratio
       gamma=0.5,  # Weight adjustment exponent
       fit_intercept=True,
       precompute=False,
       max_iter=10000,
       copy_X=True,
       solver=None,
       tol=None,
       positive=False,
       positive_tol=1e-3,
       random_state=None,
       eps_coef=1e-6,
       verbose=True
   ):
       params_asgl = dict(model="lm", penalization="asgl")
       if solver is not None:
           params_asgl["solver"] = solver
       if tol is not None:
           params_asgl["tol"] = tol
       super().__init__(**params_asgl)
       self.alpha1 = alpha1
       self.alpha2 = alpha2
       self.l1_ratio1 = l1_ratio1
       self.l1_ratio2 = l1_ratio2
       self.gamma = gamma
       self.fit_intercept = fit_intercept
       self.max_iter = max_iter
       self.precompute = precompute
       self.copy_X = copy_X
       self.positive = positive
       self.positive_tol = positive_tol
       self.random_state = random_state
       self.eps_coef = eps_coef
       self.verbose = verbose
       if not self.fit_intercept:
           raise NotImplementedError
   def fit(self, X, y, check_input=True):
       if check_input:
           X_copied = self.copy_X and self.fit_intercept
           X, y = self._validate_data(
               X,
               y,
               accept_sparse="csc",
               order="F",
               dtype=[np.float64, np.float32],
               copy=X_copied,
               multi_output=True,
               y_numeric=True,
           )
       # 第一段階の ElasticNet 実行
       enet_model_1 = self.elastic_net(self.l1_ratio1, alpha=self.alpha1)
       enet_model_1.fit(X, y)
       enet_coef = enet_model_1.coef_
       # 重みの計算
       weights = 1.0 / (np.maximum(np.abs(enet_coef), self.eps_coef) ** self.gamma)
       # 第二段階の最適化
       self.coef_, self.intercept_ = self._optimize_second_stage(X, y, weights)
       # モデル属性を格納
       self.enet_coef_ = enet_coef
       self.weights_ = weights
       return self
   def predict(self, X):
       check_is_fitted(self, ["coef_", "intercept_"])
       return super(ElasticNet, self).predict(X)
   def elastic_net(self, l1_ratio, **params):
       """
       Create an ElasticNet model with the specified parameters.
       """
       elastic_net = ElasticNet(l1_ratio=l1_ratio)
       for k, v in self.get_params().items():
           try:
               elastic_net = elastic_net.set_params(**{k: v})
           except ValueError:
               pass  # Ignore parameters not supported by ElasticNet
       elastic_net = elastic_net.set_params(**params)
       return elastic_net
   def _optimize_second_stage(self, X, y, weights):
       """
       Perform second-stage optimization with adaptive weights.
       Returns
       -------
       coef : np.array, shape (n_features,)
       intercept : float
       """
       n_samples, n_features = X.shape
       beta_variables = [cvxpy.Variable(n_features)]
       model_prediction = 0.0
       if self.fit_intercept:
           beta_variables = [cvxpy.Variable(1)] + beta_variables
           ones = cvxpy.Constant(np.ones((n_samples, 1)))
           model_prediction += ones @ beta_variables[0]
       # モデル予測
       model_prediction += X @ beta_variables[1]
       error = cvxpy.sum_squares(y - model_prediction) / (2 * n_samples)
       l1_coefs = self.alpha2 * self.l1_ratio2
       # 第二段階の正則化項
       l1_penalty = cvxpy.Constant(l1_coefs * weights) @ cvxpy.abs(
           beta_variables[1]
       )
       l2_penalty = (
           cvxpy.Constant(self.alpha1 * (1 - self.l1_ratio1))
           * cvxpy.sum_squares(beta_variables[1])
       )
       constraints = [b >= 0 for b in beta_variables] if self.positive else []
       # 最適化問題の定義
       problem = cvxpy.Problem(
           cvxpy.Minimize(error + l1_penalty + l2_penalty), 
   constraints=constraints
       )
       problem.solve(solver="OSQP", max_iter=self.max_iter)
       if problem.status != "optimal":
           raise ConvergenceWarning(
               f"Solver did not reach optimum (Status: {problem.status})"
           )
       beta_sol = np.concatenate([b.value for b in beta_variables], axis=0)
       beta_sol[np.abs(beta_sol) < self.tol] = 0
       intercept, coef = beta_sol[0], beta_sol[1:]
       coef = np.maximum(coef, 0) if self.positive else coef
       return coef, intercept
   import pandas as pd
   from sklearn.preprocessing import StandardScaler
   # CSVファイルのパスを指定してください
    file_path = "C:/Users/monda/OneDrive/デスクトップ/統計正しいほう/変数選択/ 
   削除_変数作成後.csv"
   # CSVファイルの読み込み
   df = pd.read_csv(file_path)  # 文字コードが異なる場合は、'utf-8' を他のエン 
   コーディングに変更してください
   # 特徴量とターゲットの分割
   X = df.drop(columns=['取引価格（㎡単価）'])  # 取引価格（㎡単価）をyとして分 
   離
   y = df['取引価格（㎡単価）']  # ターゲット変数
   # 🔥 標準化の実施
   scaler_X = StandardScaler()
   scaler_y = StandardScaler()
   X_scaled = scaler_X.fit_transform(X)  # Xの標準化
   y_scaled = scaler_y.fit_transform(y.values.reshape(-1, 1)).ravel()  # yの標 
   準化
   model = AdaptiveElasticNet()
   model.fit(X_scaled, y_scaled)
   # 係数の表示（重要な特徴量を確認）
   coef = model.coef_
   print("特徴量の係数:", coef)
   # 係数がゼロでない特徴量を選定
   selected_features = X.columns[coef != 0]
   print("選定された特徴量:", selected_features)
   # 係数と特徴量名をデータフレームに変換
   selected_features_df = pd.DataFrame({
   'Feature': X.columns,
   'Coefficient': coef
   })
   # 係数がゼロでない特徴量のみをフィルタリング
   selected_features_df = selected_features_df[selected_features_df['Coefficient'] != 0]
   # CSVファイルとして保存
   output_file_path = "選ばれた変数.csv"  # 保存先を指定してください
   selected_features_df.to_csv(output_file_path, index=False, encoding='utf-8-sig')  
   # 日本語対応のエンコーディング
   print(f"選定された特徴量と係数を {output_file_path} に保存しました。")
   #ref(変数選択.png,,50%)

AENの数式は以下のようなものを扱う。そして、今回行ったハイパーパラメータは以下のものである。 AENの推計は，二段階で行われる．まず，予備推計としてEN で係数を推定する．その うえで，係数の絶対値が小さい変数に対してより強い罰則を与えるようL1 ノルムの正則化 項を変数ごとに調整したうえで，改めてEN を実施する．

選択された変数は以下のようになる

要因分析 

手順 

「hedo.ipynb」を実行すると、有意水準に基づいて有意性を判定する関数

   import pandas as pd
   # Excelファイルを読み込み
   file_path = "c:/Users/nt011/Desktop/統計正しいほう/要因分析/位置情報_削 
   除.csv"
   df = pd.read_csv(file_path)
   # 住所が含まれている列を特定して、住所ごとの出現回数を確認
   address_column = '住所'  # 住所列の名前がわからない場合はここを変更してくだ 
   さい
   if address_column in df.columns:
   address_counts = df[address_column].dropna().value_counts()  # 住所ごとの出 
   現回数を集計
   # 出現回数をデータフレームに変換
   address_counts_df = address_counts.reset_index()
   address_counts_df.columns = ['住所', '出現回数']  # 列名を設定
   # 結果をCSVファイルに書き出し
   output_file = "c:/Users/nt011/Desktop/統計正しいほう/要因分析/住所出現回数_ 
   削除.csv"
   address_counts_df.to_csv(output_file, index=False, encoding='utf-8-sig')  
   # UTF-8エンコーディングで保存
   print(f"住所ごとの出現回数がCSVファイルに保存されました: {output_file}")
   else:
   print("住所列が見つかりません")

回帰係数、t値、p値、有意性を格納したリストを「回帰分析結果_削除.csv」に生成します。

   import numpy as np
   import pandas as pd
   from sklearn.preprocessing import StandardScaler
   # ファイルパスの指定
   selected_features_path = "選ばれた変数_削除.csv"
   all_data_path = "c:/Users/nt011/Desktop/統計正しいほう/要因分析/削除_変数作 
   成後.csv"
   location_data_path = "c:/Users/nt011/Desktop/統計正しいほう/要因分析/位置情 
   報_削除.csv"
   output_data_path = "選択された特徴量とデータ_削除.csv"
  # 選定された変数のCSVファイルを読み込み
   selected_features_df = pd.read_csv(selected_features_path, encoding='utf-8-sig')
   selected_features = selected_features_df['Feature'].tolist()
   # 元データを再度読み込み
   all_data_df = pd.read_csv(all_data_path)
   # 選定された説明変数と目的変数を抽出
   df = all_data_df[selected_features + ['取引価格（㎡単価）']]
   # 位置情報データを読み込み
   location_df = pd.read_csv(location_data_path, usecols=['住所'])
   # 元データと位置情報を結合（共通のキーがある場合）
   df = pd.concat([df, location_df], axis=1)
   # 統合データをCSVファイルとして保存
   df.to_csv(output_data_path, index=False, encoding='utf-8-sig')
   print(f"選定された変数と位置情報を含むデータを {output_data_path} に保存しま 
   した。")

   import pandas as pd
   # ファイルパス
   count_file = "c:/Users/nt011/Desktop/統計正しいほう/要因分析/住所出現回数_削除.csv"
   data_file = "C:/Users/nt011/Desktop/統計正しいほう/要因分析/選択された特徴量とデータ_削除.csv"
   output_file = "C:/Users/nt011/Desktop/統計正しいほう/要因分析/ソート済みデータ_削除.csv"
   # CSVファイルの読み込み
   count_df = pd.read_csv(count_file, encoding='utf-8-sig')
   data_df = pd.read_csv(data_file, encoding='utf-8-sig')
   # 住所列の名前
   address_column = '住所'
   # 緯度・経度を除外（該当カラムがある場合）
   drop_columns = ['緯度', '経度']
   data_df = data_df.drop(columns=[col for col in drop_columns if col in 
   data_df.columns])
   # 住所の出現回数をデータに統合
   merged_df = data_df.merge(count_df, on=address_column, how='left')
   # 出現回数の多い順にソート
   sorted_df = merged_df.sort_values(by='出現回数', ascending=False).drop(columns=['出現回数'])
   # ソート済みデータをCSVに保存
   sorted_df.to_csv(output_file, encoding='utf-8-sig', index=False)
   print(f"ソート済みデータを {output_file} に保存しました。")
   import pandas as pd
   import numpy as np
   from sklearn.linear_model import LinearRegression
   from sklearn.preprocessing import StandardScaler
   import statsmodels.api as sm
   # データの読み込み
   file_path = "C:/Users/nt011/Desktop/統計正しいほう/要因分析/ソート済みデータ_削除.csv"
   data = pd.read_csv(file_path)
   # 目的変数と説明変数の設定
   y = data['取引価格（㎡単価）']  # 取引価格を目的変数に
   X = data.drop(columns=['取引価格（㎡単価）', '住所'])  # 取引価格と住所以外を説明変数に
   # 🔥 標準化の実施
   scaler_X = StandardScaler()
   scaler_y = StandardScaler()
   X_scaled = scaler_X.fit_transform(X)  # Xの標準化
   y_scaled = scaler_y.fit_transform(y.values.reshape(-1, 1)).ravel()  # yの標 
   準化
   # 線形回帰モデルの作成
   model = LinearRegression()
   model.fit(X_scaled, y)
   # 回帰分析の結果を表示
   X_scaled_with_const = sm.add_constant(X_scaled)  # 定数項を追加
   ols_model = sm.OLS(y, X_scaled_with_const).fit()
   # 必要な情報を抽出
   r_squared = ols_model.rsquared  # 決定係数
   coefficients = ols_model.params  # 回帰係数（定数項も含む）
   t_values = ols_model.tvalues  # t値（定数項も含む）
   p_values = ols_model.pvalues  # p値（有意水準を判定するため）
   variable_names = ['定数項'] + list(X.columns)
   # 有意水準に基づいて有意性を判定する関数
   def significance(p_value):
   if p_value < 0.01:
       return "1% 有意"
   elif p_value < 0.05:
       return "5% 有意"
   elif p_value < 0.10:
       return "10% 有意"
   else:
       return "有意でない"
   # 結果を格納するリスト
   results = []
   # 回帰係数、t値、p値、有意性をリストに格納
   for var_name, coef, t_val, p_val in zip(variable_names, coefficients, 
   t_values, p_values):
   significance_level = significance(p_val)
   results.append([var_name, coef, t_val, p_val, significance_level])
   # 結果をDataFrameに変換
   results_df = pd.DataFrame(results, columns=["変数", "回帰係数", "t値", "p 
   値", "有意性"])
   # 予測値の計算
   data["標準化予測取引価格（㎡単価）"] = 
   ols_model.predict(X_scaled_with_const)
   data["標準化実測取引価格（㎡単価）"] = y_scaled
   # **結果を表示**
   print("
   ==== 回帰分析結果 ====")
   print(f"決定係数 (R^2): {r_squared:.4f}
   ")
   print(results_df)
   # 予測結果を含めたデータを保存
   predicted_output_file_path = "C:/Users/nt011/Desktop/統計正しいほう/要因分 
   析/予測結果_削除.csv"
   data.to_csv(predicted_output_file_path, index=False, encoding="utf-8-sig")
   # CSVファイルとして保存
   output_file_path = "C:/Users/nt011/Desktop/統計正しいほう/要因分析/回帰分析 
   結果_削除.csv"
   results_df.to_csv(output_file_path, index=False, encoding="utf-8-sig")
   print(f"回帰分析結果が{output_file_path}に保存されました。")

寄与率の算出 

寄与率は下の数式と変数で算出を行う。
実行を行う前に寄与率フォルダ内に町ごとの寄与率フォルダを事前に作成する。
寄与率.ipynbを上から実行する。

AENで選ばれた一次項、二次項、交互作用項から一次項、二次項の説明変数を抽出して、その後、交互作用項の計算を行って、余計な変数を削除する。

   import pandas as pd
   import os

   # データの読み込み
   data_file = "C:/Users/tn011/Desktop/統計正しいほう/寄与率/町ごとの平均値_削 
   除.csv"
   location_file = "C:/Users/tn011/Desktop/統計正しいほう/寄与率/抽出した住所緯 
   度経度_削除.csv"
   output_file = "C:/Users/tn011/Desktop/統計正しいほう/寄与率/まとめたデータ_ 
   削除.csv"

   df = pd.read_csv(data_file)
   location_df = pd.read_csv(location_file)

   # 交互作用項の計算
   交互作用項 = [
   ('最寄駅：距離（分）', 'レジャー施設'),
   ('最寄駅：距離（分）', '密度（道路）'),
   ('面積（㎡）', '前面道路：幅員（ｍ）'),
   ('居住年数5年未満の人口割合', '小学校'),
   ('一戸建て割合', 'min_distance_to_crime'),
   ('アパート・低中層マンション割合', 'min_distance_to_crime'),
   ('ホテル', '建物の面積比率'),
   ('デパート', '平均次数（道路）'),
   ('中学校', '建物の面積比率'),
   ('保育園／幼稚園', 'min_distance_to_crime'),
   ('最寄り駅名_富山', '最終学歴が大学以上の人口割合'),
   ('土地の形状_不整形', '居住年数が20年以上の人口割合'),
   ('建物の構造_木造', '面積（㎡）'),
   ('建物の構造_木造', '高校'),
   ('用途_住宅', '延床面積（㎡）'),
   ('前面道路：方位_南東', 'エッジ数（道路）'),
   ('都市計画_市街化調整区域', '中学校'),
   ]

   # 交互作用項を計算し、列を追加
   for var1, var2 in 交互作用項:
   df[f'{var1}×{var2}'] = df[var1] * df[var2]

   # 不要な要素のリスト
   要素リスト = [
   '最寄駅：距離（分）', 'レジャー施設', '密度（道路）', '面積（㎡）', '前面道 
   路：幅員（ｍ）', '居住年数5年未満の人口割合',
   '小学校', '一戸建て割合', 'min_distance_to_crime', 'アパート・低中層マンショ 
   ン割合', 'ホテル',
   '建物の面積比率', 'デパート', '平均次数（道路）', '中学校', '保育園／幼稚 
   園',
   '最寄り駅名_富山', '最終学歴が大学以上の人口割合', '土地の形状_不整形', '居 
   住年数が20年以上の人口割合', '建物の構造_木造',
   '高校', '用途_住宅', '延床面積（㎡）', '前面道路：方位_南東', 'エッジ数（道 
   路）',
   '都市計画_市街化調整区域'
   ]

   # 不要な列を削除
   df.drop(columns=要素リスト, inplace=True)

   # 住所をキーにして緯度・経度をマージ
   df = df.merge(location_df, on="住所", how="left")

   # 取引価格と緯度・経度を削除
   df.drop(columns=['取引価格（㎡単価）', '緯度', '経度'], inplace=True)

   # CSVの列順を指定（1列目: 住所, 残りの変数）
   columns_order = ['住所'] + [col for col in df.columns if col != '住所']

   # 指定した列順で保存
   df = df[columns_order]
   df.to_csv(output_file, index=False, encoding='utf-8-sig')

   print("データのまとめが完了し、CSVファイルに保存しました。")

上の数式を再現したpythonコードが下になります。

   import os
   import pandas as pd
   import numpy as np
   # ======== CSVファイルの読み込み ========
   df_selected_features = pd.read_csv("C:/Users/tn011/Desktop/統計正しいほう/寄 
   与率/選ばれた変数_削除.csv")
   df_toyama_data = pd.read_csv("C:/Users/tn011/Desktop/統計正しいほう/寄与率/ 
   まとめたデータ_削除.csv")
   # ======== selected_features.csv から 回帰係数と交互作用項の係数を取得 
   ========
   beta_values = df_selected_features.set_index('Feature') 
   ['Coefficient'].to_dict()
   # 交互作用項の係数を取得（交互作用項に '×' が含まれる）
   interaction_coeffs = {key: value for key, value in beta_values.items() if 
   '×' in key}
   # 交互作用項を除外した beta_values の作成
   beta_values_without_interaction = {key: value for key, value in 
   beta_values.items() if '×' not in key}
   # ======== 寄与率の計算関数 ========
   def compute_contribution(x_values, beta_values_without_interaction, 
   interaction_coeffs, target_var):
   beta_0 = 0  # 切片は 0 で固定
   delta_x = 1  # 対象変数を 1 増加
   # 1. 基準の予測値 log(y) の計算
   log_y = beta_0
   log_y += sum(beta_values_without_interaction[var] * x_values[var] for var in 
   beta_values_without_interaction if var in x_values)
   log_y += sum(interaction_coeffs[interaction] * x_values[interaction] for 
   interaction in interaction_coeffs if interaction in x_values)
   # 2. 対象変数を 1 増加させた場合の log(y')
   x_values_new = x_values.copy()
   x_values_new[target_var] += delta_x
   log_y_new = beta_0
   log_y_new += sum(beta_values_without_interaction[var] * x_values_new[var] 
   for var in beta_values_without_interaction if var in x_values_new)
   log_y_new += sum(interaction_coeffs[interaction] * x_values_new[interaction] 
   for interaction in interaction_coeffs if interaction in x_values_new)
   # 3. 寄与率の計算
   y_original = np.exp(log_y)
   y_new = np.exp(log_y_new)
   contribution_rate = ((y_new - y_original) / y_original) * 100
   return contribution_rate
   # ======== すべての住所に対する寄与率計算 ========
   output_folder = "C:/Users/tn011/Desktop/統計正しいほう/寄与率/町ごとの寄与 
   率"
   os.makedirs(output_folder, exist_ok=True)
   # すべての住所を取得
   target_addresses = df_toyama_data['住所'].unique()
   # 各住所に対する寄与率を計算
   for target_address in target_addresses:
   # 住所ごとのデータをフィルタリング
   target_data = df_toyama_data[df_toyama_data['住所'] == target_address]
   # 数値カラムのみを抽出して平均を計算
   target_data_numeric = target_data.select_dtypes(include=[np.number]).mean()
   # 住所ごとのデータを取得
   x_values = target_data_numeric.to_dict()
   # 寄与率を計算する対象の変数を取得
   target_vars = list(x_values.keys())
   # 各住所の変数ごとの寄与率を計算
   contributions = {var: compute_contribution(x_values, 
   beta_values_without_interaction, interaction_coeffs, var) for var in 
   target_vars}
   # 結果をデータフレームにまとめる
   result_df = pd.DataFrame([{'住所': target_address, **contributions}])
   # CSVファイルとして保存
   result_file_path = os.path.join(output_folder, f"{target_address}.csv")
   result_df.to_csv(result_file_path, index=False, encoding='utf-8-sig')
   print(f"{target_address} の寄与率が保存されました: {result_file_path}")

「町ごとの寄与率」フォルダーにそれぞれの町の寄与率のcsvファイルが生成される。

寄与率のGISによる可視化 

ドライブにあるWeb-GISフォルダを統計正しいほうのフォルダに入れてもらってWeb-GIS.pyを実行する。

おわり