回帰によるブルーベリー（野生種）の収穫量予測

はじめに

2025年2月から半年間、AIデータサイエンスに関して学んできました。最終課題は何を分析して知見を得てみようかと考えていたところ、ブルーベリーの収穫量予測というkaggleコンペティションを発見。自宅で父がベリー系統の家庭菜園をしていることもあり、知識の一助になるやもと考え回帰による分析結果をまとめてみました。

解決したい社会課題

現代日本は第一次産業に従事する若い労働者が足りていない状況下にあると言われている中、収穫量や収入が予測できるデータを作れば、第一次産業に参入する人々への足掛かりになると考え（今回はブルーベリーとピンポイントですが）データを分析してみることにしました。

実行環境

パソコン：Inspiron 14 5440
開発環境：Kaggle Notebook
言語：Python
ライブラリ：Matplotlib, Seaborn, Pandas, Numpy,
Warnings, OS, Scikit-learn, XGB

分析するデータ

以下のKaggelコンペティションから取得したデータを使って分析を実行します。
Regression with a Wild Blueberry Yield Dataset

分析の流れ

データセットの取得と確認
データセットの前処理
予測モデルを作成し、平均絶対誤差で評価
作成した予測モデルの特徴量重要度の可視化

データセットの取得と確認

まず最初に必要なライブラリと上記のkaggleコンペが用意したCSVを読み込みます。

実行したコード

# Kaggle Notebookの最初のセルを実行しセッションスタート
import numpy as np 
import pandas as pd 
import os
for dirname, _, filenames in os.walk('/kaggle/input'):
    for filename in filenames:
        print(os.path.join(dirname, filename))
        
# グラフに日本語を使用するためjapanize-matplotlibをインストール
!pip install japanize-matplotlib

# 残りの必要なライブラリをインポートする
import warnings
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib

# CSVファイルを読み込む
train = pd.read_csv("../input/playground-series-s3e14/train.csv")
test = pd.read_csv("../input/playground-series-s3e14/test.csv")

続いて読み込んだデータの状態をそれぞれ確認します。

実行したコード

# 学習データの状態を確認する
train.info()
train.head()
train.tail()

# テストデータの状態を確認する
test.info()
test.head()
test.tail()

実行結果（学習データ）

※一部省略

実行結果（テストデータ）

※一部省略

確認してわかったこと

info()の情報からどちらも欠損値が無く、データのカテゴリも全て数値型のため複雑な前処理は必要なさそうです。

データセットの前処理

次に、前処理が必要な特徴量を把握するためにヒストグラムと散布図、箱ひげ図で可視化します。

実行したコード

# FutureWarningが出現するので無視する
warnings.filterwarnings("ignore", category=FutureWarning)

# 今回のデータは全て数値型なのでId列のみを除いて取得
numerical_cols = train.drop("id", axis=1, errors="ignore").columns

# グラフの列数と行数を定義
n_cols = 3
n_rows = (len(numerical_cols) + n_cols -1) // n_cols

# グラフの大きさを指定し、ループ処理をしやすくするためにaxesを一次配列にする
fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 5, n_rows * 4))
axes = axes.flatten()

for i, col in enumerate(numerical_cols):
    sns.histplot(train[col], bins=30, kde=True, ax=axes[i])
    axes[i].set_title(f"{col}の頻度数", fontsize=10)
    axes[i].set_xlabel(col, fontsize=8)
    axes[i].set_ylabel("頻度", fontsize=8)
    axes[i].grid(axis="y", alpha=0.75)

plt.tight_layout()
plt.show()

# 収穫量(yield)とId列を除いた特徴量を獲得
features = train.drop(["id","yield"], axis=1, errors="ignore" ).columns

# 全体の図と各サブプロットの軸を作成
fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 6, n_rows * 5))
axes = axes.flatten()

# 散布図を作成
for i, col in enumerate(features):
    sns.regplot(x=train[col], y=train["yield"], ax=axes[i], line_kws={"color": "red"})
    axes[i].set_title(f"{col}と収穫量の相関関係", fontsize=10)
    axes[i].set_xlabel(col, fontsize=8)
    axes[i].set_ylabel("収穫量", fontsize=8)

plt.tight_layout()
plt.show()

# 全体の図と各サブプロットの軸を作成
fig, axes = plt.subplots(n_rows, n_cols, figsize=(n_cols * 6, n_rows * 5))
axes = axes.flatten()

# 箱ひげ図を作成
for i, col in enumerate(features):
    sns.boxplot(train[col], ax=axes[i])
    axes[i].set_title(f"{col}のデータ分布", fontsize=10)
    axes[i].set_ylabel(f"{col}の値", fontsize=8)

plt.tight_layout()
plt.show()

実行結果

図1. 各特徴量のヒストグラム

ヒストグラムから得られた情報と考察

・分布の確認と前処理の判断
目的変数であるyield列、説明変数のfruitset、fruitmass、seedsの各列は歪みがなく釣鐘型に近いため、前処理は必要ないと判断しました。
一方で、honeybee列は極端に大きいデータが含まれており、右に歪んだ分布を示しています。

・気温データの考察
最高気温と最低気温に関してはおそらく華氏で記録されていると推測されます。当てはめると、華氏95度は摂氏約35度、華氏24度は摂氏約-20度となります。
ワイルドブルーベリーの主な産地がアメリカ北部からカナダにかけてである点を踏まえると、これらの気温は現地の気候（同緯度にはスウェーデンなど亜寒帯の国が多い）と捉えることができ、現実的な記録だと解釈できます。

図2. 各特徴量と収穫量の散布図

散布図から得られた情報と考察

・clonesizeの相関と考察
clonesizeは大きいほど収穫量は下がる傾向があります。これはブルーベリーが自己不和合性(同じ品種同士だと受粉が成功しない)という特性からなり、クローン株が大きいほど別の品種に受粉する可能性が下がるため、負の相関があるのは正しいと判断できます。

・ポリネーター（受粉媒介者）の相関と考察
ポリネーター（花粉を運ぶことで受粉活動の手助けをする動物の総称）の場合は密度が高いほど収穫量は挙がる傾向が殆どですが、ミツバチだと下がる傾向があります。
2012年の鹿児島大学農学部の実験で、一安定多収栽培を目指したミツバチ放飼効果の検討というものがあり、ミツバチを放逐したことにより前年の約2倍もの収穫量を記録したという実験結果が記述されています。この実験からもミツバチとブルーベリーは高い正の相関があると私も考えていたのですが、可視化した結果は逆となりました。
理由として、このデータは日本ではなくおよそアメリカ北部からカナダにかけてのデータであると解釈したことから、気候が亜寒帯に属する地域のため気候の差が影響を与えていると推測されます。
事実、セイヨウミツバチは気温12℃を下回ると越冬行動（生物が冬の環境を生き延びるために行う様々な行動）に入り、活動が鈍くなる特徴を持っています。最低気温のデータに氷点下の気温も記録されているため、気候差による影響で相関が下がっていると考えられます。

・気温データと降雨日数記録の相関と考察
気温データに関しては若干負の相関を確認できる程度の情報です。何か他の特徴量を抽出した方がいいかもしれません。
降雨日数と平均降雨日数は数値が高いと収穫量は下がる傾向があります。これはブルーベリーが雨による裂果（水分を急激に吸収し実が爆ぜる現象）を起こしやすい果物であることからデータとしては正しいと言えます。

・その他の特徴量の相関と考察
残りのデータに関しては収穫量に直結するデータ群であるため正しい正の相関があると判断できます。

図3. 各特徴量の箱ひげ図

箱ひげ図から得られた情報と考察

・ポリネーターの外れ値処理の検討
honeybee列はかなり極端な分布をしていることが分かりました。
しかし、農業的背景を踏まえると農園と並行して養蜂も同時に営んでいる場合、ミツバチの巣が敷地内にあれば空間密度が高くなるのは容易に想像できます。
農園の規模に対してミツバチの巣の量が多ければ高密度になる可能性は高く、異常値と判断せず、ユニークな情報を持った数値として扱うべきだと判断します。
honeybee列のデータだけ大きい数値があるのも、養蜂はセイヨウミツバチが主流のため説明できます。
そのため、削除は行わず対数変換による処理を検討します。
他のポリネーターのデータに関しても、その地域ではあまり見ない、この地域では出現頻度が高いなど、分布地域により変わる可能性は十分にあり得るため異常値はないと判断します。

・降雨日数の外れ値処理の検討
降雨日に関しても地域差による違いは高確率でありうると考えられます。乾燥した地域に農園が位置すればありうる数値のため、外れ値の削除はしない方向で進みます。

その他特徴量の外れ値処理の検討
fruitset、fruitmass、seedsの各列に関しては、農園の規模による差ではないかと推測されます。大規模な農園もあれば、メインで育てている果樹ではなく、もののついでに育てている可能性も否定できません。
以上の理由から、他の農園より数値が低いことも情報の損失をしないために、各列の外れ値の処理は必要ないと判断します。

honeybee列の対数変換の実行

# ミツバチ(honeybee)列を対数変換する
train["honeybee_log"] = np.log1p(train["honeybee"])

# 対数変換後のhoneybee列を散布図で可視化する
plt.figure(figsize=(8, 6))
sns.regplot(x=train["honeybee_log"], y=train["yield"], line_kws={"color": "red"})
plt.title("対数変換後のミツバチと収穫量の相関関係")
plt.xlabel("ミツバチの密度(対数変換後)")
plt.ylabel("収穫量")
plt.grid(True, linestyle="--", alpha=0.7)
plt.show()

実行結果

完全に歪みは解消されてはいませんが、対数変換を実行したことによってデータの特徴を把握しやすくなりました。以降honeybee列は対数変換した数値を使用します。

ヒートマップの作成

特徴量のヒートマップを作成し、多重共線性が発生してないかの確認をします。

実行したコード

# ヒートマップを作成する
# この時点で必要のない列を削除しておく
train = train.drop(["id", "honeybee", ], axis=1, errors="ignore")

# 相関行列を計算する
corr_matrix = train.corr()

# ヒートマップを作成する
plt.figure(figsize=(12, 10))
sns.heatmap(corr_matrix, annot=True, cmap="coolwarm", fmt=".2f", linewidth=.5)
plt.title("学習データの相関関係ヒートマップ")
plt.show()

実行結果

・clonesizeとhoneybee_logの相関関係と対策
clonesize、honeybee_logの正の相関が大きいです。大規模な農園ほどclonesizeが大きくなり、養蜂の規模も大きくなるという事でしょうか。
この二つに関しては積を計算し一つにまとめ、新たにinteraction列にします。
気温データの相関関係と対策
気温データはお互いが気温記録に対して高い正の相関を持っています。
どれか一つを残し、それ以外の特徴量を削除することも考えましたが、今回は最高平均気温AverageOfUpperTRangeと最低平均気温AverageOfLowerTRangeの差を計算し、平均気温差AvgTempRangeという新しい特徴量を作成します。寒暖差は果物の糖度などに大きく影響を与える傾向があります（日中気温が高いほど光合成で糖分が合成され、夜間気温が低いほど呼吸量が減り、果実に糖分が蓄積される）。収穫量に影響を与える可能性が高いと判断し、今回の処理は平均気温差で学習します。
・降雨データの相関関係と対策
降雨データも互いに高い正の相関があります。降雨日数の方が収穫量への影響があると判断し、降雨日数記録を残し、平均降雨日数のデータは削除します。
・その他の相関関係と対策
fruitset、fruitmass、seeds各列も互いに高い正の相関があります。これらに関してはデータが持つ情報の重要度が高いため、前処理は必要ないと判断します。

新たな特徴量を追加する

上記の前処理以外に、新しい特徴量として最も密度の高いポリネーターを判別するDominantPollinator列を追加します。
ポリネーターは花粉を花に運び受粉させる重要な役割を持った動物です。しかしながら、受粉活動がみな同じような方法で行っているわけではありません。例えば、bumblesマルハナバチは振動受粉（口吻を花に刺してぶら下がり蜜を集める手法）をする種であり、花を振動させることでより遠くへ花粉を飛ばすことが可能です。植物により適しているポリネーターが変わってくるため、農園ごとにどのポリネーターが一番多いのかという情報は、収穫量に大きな影響を与えると推測します。

実行するコード

# 複数のポリネーターの中で一番密度が高いポリネーターを記録した特徴量を作成する
def get_dominant_pollinator_corr_priority(row):
    pollinators = {
        "honeybee_log": row["honeybee_log"],
        "bumbles": row["bumbles"],
        "osmia": row["osmia"],
        "andrena": row["andrena"]
    }
    # 最大値を見つける
    max_count = max(pollinators.values())
    # 最大値を持つ蜂の種類をリストアップする
    dominant_types = [bee_type for bee_type, count in pollinators.items() if count == max_count]
    # 優先順位を設定（yield列への正の相関が大きい順）
    if "osmia" in dominant_types:
        return "osmia"
    elif "bumbles" in dominant_types:
        return "bumbles"
    elif "andrena" in dominant_types:
        return "andrena"
    elif "honeybee_log" in dominant_types:
        return "honeybee_log"
    else:
        return "none"

# 学習データに適用
train["DominantPollinator"] = train.apply(get_dominant_pollinator_corr_priority, axis=1)

# このカテゴリ列をワンホットエンコーディング化する
train = pd.get_dummies(train, drop_first=False)

# clonesize列とhoneybee列をまとめたinteraction列を作成
train["interaction"] = train["clonesize"] * train["honeybee_log"]

# 新しい気温の特徴量である平均気温の範囲("AvgTempRange"）という列を作成
train["AvgTempRange"] = train["AverageOfUpperTRange"] - train["AverageOfLowerTRange"]

# 必要のなくなった列を削除する
train = train.drop(["clonesize", "honeybee", "honeybee_log", "MaxOfUpperTRange", "MinOfUpperTRange", "AverageOfUpperTRange", "MaxOfLowerTRange", "MinOfLowerTRange", "AverageOfLowerTRange", "AverageRainingDays"], axis=1, errors="ignore")

# データ前処理後のヒートマップを作成
corr_matrix = train.corr()
plt.figure(figsize=(12, 10))
sns.heatmap(corr_matrix, annot=True, cmap="coolwarm", fmt=".2f", linewidth=.5)
plt.title("整理後の相関関係ヒートマップ")
plt.show()

実行結果

多重共線性はかなり解消されました。
残りのfruitmass, seeds, fruitsetに関しては、無理に統合するとそれぞれが持つ特徴を捉えられなくなる可能性が高いと判断したため,前処理はせずにモデルへ学習させます。
そのため、今回のモデルには多重共線性の影響を受けにくいXGBoost回帰モデルを採用します。

予測モデルを作成し、平均絶対誤差で評価

ハイパーパラメータチューニングで最適なパラメータを探索し、平均絶対誤差MAEで評価します。
平均絶対誤差は予測値と正解値の誤差の絶対値を計算し、その総和の平均を出力します。値が小さければ小さいほどモデルとして優れていると判断できる指標です。
過学習を起こしていないか判断するためにtrain_test_split関数も使用します。

実行するコード

# ライブラリのインポート
from sklearn.model_selection import train_test_split, GridSearchCV
import xgboost as xgb
from sklearn.metrics import mean_absolute_error

# 特徴量と目的変数に分割
X = train.drop("yield", axis=1, errors="ignore")
y = train["yield"]

# 学習データを訓練データと検証データに分割する
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.3, random_state=0)

# XGBモデルのインスタンス化を実行
model = xgb.XGBRegressor(objective="reg:absoluteerror", random_state=0, n_jobs=-1)

# 探索するハイパーパラメータを定義する
param_grid = {"n_estimators": [800, 850, 900],
              "learning_rate": [0.01, 0.03, 0.05],
              "max_depth": [5, 6, 7],
              "subsample": [0.8, 1.0],
              "colsample_bytree": [0.8, 1.0]}

# グリッドサーチをインスタンス化する
# scoring="neg_mean_absolute_error"を指定し、MAEが最小になるパラメータを探索する
grid_search = GridSearchCV(estimator=model,
                           param_grid=param_grid,
                          scoring="neg_mean_absolute_error",
                          cv=5, n_jobs=-1)

# グリッドサーチを実行
grid_search.fit(X_train, y_train)

print(grid_search.best_params_, grid_search.best_score_)

実行結果

・ベストパラメータ
{'colsample_bytree': 1.0,
'learning_rate': 0.03,
'max_depth': 5,
'n_estimators': 850,
'subsample': 0.8}

・ベストスコア
-344.5554691879012

グリッドサーチはベストスコアが最も大きい物を評価するため、neg_で逆転させマイナスから0へ近づけさせる処理にしています。

そのためMAEは344.5554691879012となります。

バリデーションデータで評価

# 最適なモデルを取得
best_model = grid_search.best_estimator_

# テストデータで予測を実行
y_pred = best_model.predict(X_val)

# 予測結果の平均絶対誤差を計算する
mae = mean_absolute_error(y_val, y_pred)

print(mae)

実行結果

344.4279571434898

検証スコアとベストスコアの差がほとんどないため、過学習は起きていないと判断できます。

kaggleコンペティションへ投稿し最終評価を確認する

# 学習データに行った前処理をテストデータにも実行する。
test["honeybee_log"] = np.log1p(test["honeybee"])
test["interaction"] = test["clonesize"] * test["honeybee_log"]
test["AvgTempRange"] = test["AverageOfUpperTRange"] - test["AverageOfLowerTRange"]
test["DominantPollinator"] = test.apply(get_dominant_pollinator_corr_priority, axis=1)
test = pd.get_dummies(test, drop_first=False)

# 予測に必要な特徴量のみを抽出する
test_features = test[X.columns]

# 予測を実行
test_pred = best_model.predict(test_features)

# 念のため負の値が発生してしまった場合に保険をかけておく
test_pred = np.clip(test_pred, 0, None)

# 提出ファイルの作成
submission = pd.DataFrame({"id": test["id"], "yield": test_pred})
submission.to_csv("submisson.csv", index=False)

実行結果

kaggleスコアは上記の通りになりました。
Public ScoreとPrivate Scoreに差がありますが、他のノート提出者の記録も同程度の差があるため問題なしと判断します。

作成した予測モデルの特徴量重要度の可視化

importance_typeはgain（予測性能改善にどれだけ貢献したかの指標）で可視化します。

実行するコード

# 特徴量重要度を取得
feature_importances = best_model.get_booster().get_score(importance_type="gain")

# 特徴量重要度をデータフレームに変換し、ソートする
importance_df = pd.DataFrame({"Feature": list(feature_importances.keys()), "Importance": list(feature_importances.values())})
importance_df = importance_df.sort_values(by="Importance", ascending=False)

# 特徴量重要度を可視化
plt.figure(figsize=(10, 8))
sns.barplot(x="Importance", y="Feature", data=importance_df)
plt.title("特徴量重要度")
plt.xlabel("重要度（gain）")
plt.ylabel("特徴量")
plt.tight_layout()
plt.show()

実行結果

分析したデータから得られた情報

傾向

予測モデルの重要度は7割以上fruitset列に重みがあようです。
次点でseeds列、AvgTempRange列、fruitmass列が予測に貢献しているる傾向がみられます。
ポリネーター系の特徴量はosmia列とandrena列がスコア予測に貢献しているようです。

課題

fruitset、seeds、fruitmassの各列は、本来収穫してからでないと記録できない情報です。今回のデータセットでは収穫量に直結するデータが記録されていたため、スコア予測の重要度の比重がその3種の列に寄ってしまいました。
これらのデータが無い前提で考えた際、土のph度数や野鳥、害虫対策などのさらなる環境的要因の情報が必要不可欠であると感じました。

考察

課題の項目で記述した通り、fruitset、seeds、fruitmassの各列は収穫前の生育段階の時点では知る由もない情報のため、その3つを除外して重要度を確認します。
最もスコア改善に影響しているのはAvgTempRange列（平均気温差）。寒暖差が重要だということが分かりました。寒暖差が大きいほど果実の糖度は高くなる傾向があります。
次点で重要度が高いのはRainingDays列（降雨日数）です。雨が少ない乾燥地帯の方が収穫量は増える傾向があります。
例として、フルーツトマトという野菜は果物に匹敵する糖度を持っているが、果皮が硬いという特徴をもっています。これはわざと水分を与えず、寒暖差をあたえることで優先的に果実に糖分を蓄積させているのです（糖分は水分よりも凝固点が低く凍結しづらい特徴をもつ）。その代わりに果実の水分が減るため果皮は硬くなります。
果皮が硬いという事は裂果が起きづらくなるため、ダメになる果実は減り収穫量が上昇する。平均気温差と降雨日数は収穫量を上昇させるための重要な情報であると考えられます。
osmia（ツツハナバチ）が農園に多いことも重要であるようです。ツツハナバチ（11mm）はミツバチ（13mm）よりも小さいのですが比較的毛深い種類らしく、その体毛に花粉をまとわせて受粉媒介をします。その受粉活動が最も収穫量に影響を与えていると考えられます。
最後に、interactionも影響があります。同じ種はできる限り密集させず、別種を近くに配置してミツバチの巣の数も程ほどにすることが重要だと考えられます。

まとめ

分析結果から、温室など一定の気温を保つ場所で育てず、寒暖差のある場所で生育し、程よく乾燥した環境で育てる。そのうえで、プランターなどで生育する場合は隣に別種のブルーベリーを配置することが重要であるとわかりました。ポリネーターはツツハナバチが最も優れているとわかりましたが、日本で育てることを考えるとミツバチでも十分な結果を得ることができると思います。
今回の経験でデータから背景知識を読み取ることの重要性を学びました。次回はこの経験を生かし、農学系とはまた別のデータ分析に挑戦し知見を深め、一人前のデータサイエンティストを目指したいと思います。