Main Content

回帰用の一般化加法モデルの学習

この例では、最適なパラメーターで回帰用の一般化加法モデル (GAM)に学習させる方法と、学習済みモデルの予測性能を評価する方法を示します。この例では、最初に一変量の GAM に最適なパラメーター値 (線形項のパラメーター) を特定し、次に二変量の GAM の値 (交互作用項のパラメーター) を特定します。また、この例では、特定の予測に対する項の局所的効果を調べて、予測子に対する予測の部分依存を計算することで、学習済みモデルを解釈する方法についても説明します。

標本データの読み込み

標本データ セット NYCHousing2015 を読み込みます。

load NYCHousing2015

データ セットには、2015 年のニューヨーク市における不動産の売上に関する情報を持つ 10 の変数が含まれます。この例では、これらの変数を使用して売価 (SALEPRICE) を解析します。

データ セットを前処理します。1000 ドル以下の SALEPRICE は、現金対価なしの所有権移転を示すものと仮定します。このような SALEPRICE をもつ標本を削除します。さらに、関数 isoutlier で識別される外れ値も削除します。その後、datetime 配列 (SALEDATE) を月番号に変換して、応答変数 (SALEPRICE) を最後の列に移動します。LANDSQUAREFEETGROSSSQUAREFEET、および YEARBUILT に含まれるゼロを NaN に変更します。

idx1 = NYCHousing2015.SALEPRICE <= 1000;
idx2 = isoutlier(NYCHousing2015.SALEPRICE);
NYCHousing2015(idx1|idx2,:) = [];
NYCHousing2015.SALEDATE = month(NYCHousing2015.SALEDATE);
NYCHousing2015 = movevars(NYCHousing2015,'SALEPRICE','After','SALEDATE');
NYCHousing2015.LANDSQUAREFEET(NYCHousing2015.LANDSQUAREFEET == 0) = NaN; 
NYCHousing2015.GROSSSQUAREFEET(NYCHousing2015.GROSSSQUAREFEET == 0) = NaN; 
NYCHousing2015.YEARBUILT(NYCHousing2015.YEARBUILT == 0) = NaN; 

テーブルの最初の 3 行を表示します。

head(NYCHousing2015,3)
    BOROUGH    NEIGHBORHOOD       BUILDINGCLASSCATEGORY        RESIDENTIALUNITS    COMMERCIALUNITS    LANDSQUAREFEET    GROSSSQUAREFEET    YEARBUILT    SALEDATE    SALEPRICE
    _______    ____________    ____________________________    ________________    _______________    ______________    _______________    _________    ________    _________

       2       {'BATHGATE'}    {'01  ONE FAMILY DWELLINGS'}           1                   0                1103              1290            1910          2           3e+05 
       2       {'BATHGATE'}    {'01  ONE FAMILY DWELLINGS'}           1                   1                2500              2452            1910          7           4e+05 
       2       {'BATHGATE'}    {'01  ONE FAMILY DWELLINGS'}           1                   2                1911              4080            1931          1         5.1e+05 

関数datasampleを使用して 1000 個の標本を無作為に選択し、関数cvpartitionを使用して観測値を学習セットとテスト セットに分割します。テスト用に 10% のホールドアウト標本を指定します。

rng('default') % For reproducibility
NumSamples = 1e3;
NYCHousing2015 = datasample(NYCHousing2015,NumSamples,'Replace',false);
cv = cvpartition(size(NYCHousing2015,1),'HoldOut',0.10);

学習インデックスとテスト インデックスを抽出し、学習データ セット用とテスト データ セット用の table を作成します。

tbl_training = NYCHousing2015(training(cv),:);
tbl_test = NYCHousing2015(test(cv),:);

最適なハイパーパラメーターでの GAM の学習

名前と値の引数 OptimizeHyperparameters を使用して、交差検証損失を最小化するハイパーパラメーターで GAM に学習させます。

OptimizeHyperparameters'auto' または 'all' として指定し、一変量パラメーターおよび二変量パラメーター両方に最適なハイパーパラメーター値を計算できます。あるいは、'auto-univariate' または 'all-univariate' オプションを使用して一変量パラメーターに最適な値を求めてから、'auto-bivariate' または 'all-bivariate' オプションを使用して二変量パラメーターに最適な値を求めることもできます。この例では、'all-univariate' および 'all-bivariate' を使用します。

一変量の GAM に学習させます。FitStandardDeviationtrue として指定し、応答変数の標準偏差にもモデルを当てはめます。標準偏差モデルを標準偏差推定の精度に当てはめる場合、最適なハイパーパラメーターの使用が推奨されます。名前と値の引数 InitialLearnRateForPredictorsMaxNumSplitsPerPredictor および NumTreesPerPredictor の最適な値を fitrgam で求めるため、OptimizeHyperparameters として 'all-univariate' を指定します。再現性を得るために、'expected-improvement-plus' の獲得関数を使用します。ShowPlotsfalse として、Verbose を 0 として指定し、プロットの表示とメッセージの表示をそれぞれ無効にします。

Mdl_univariate = fitrgam(tbl_training,'SALEPRICE','FitStandardDeviation',true, ...
    'OptimizeHyperparameters','all-univariate', ...
    'HyperparameterOptimizationOptions',struct('AcquisitionFunctionName','expected-improvement-plus', ...
    'ShowPlots',false,'Verbose',0))
Mdl_univariate = 
  RegressionGAM
                       PredictorNames: {'BOROUGH'  'NEIGHBORHOOD'  'BUILDINGCLASSCATEGORY'  'RESIDENTIALUNITS'  'COMMERCIALUNITS'  'LANDSQUAREFEET'  'GROSSSQUAREFEET'  'YEARBUILT'  'SALEDATE'}
                         ResponseName: 'SALEPRICE'
                CategoricalPredictors: [2 3]
                    ResponseTransform: 'none'
                            Intercept: 5.1868e+05
               IsStandardDeviationFit: 1
                      NumObservations: 900
    HyperparameterOptimizationResults: [1×1 BayesianOptimization]


  Properties, Methods

fitrgam は、最適な推定実行可能点を使用する RegressionGAM モデル オブジェクトを返します。最適な推定実行可能点は、ベイズ最適化プロセスの基となる目的関数モデルに基づいて目的関数値の信頼限界の上限を最小化するハイパーパラメーターのセットを示します。HyperparameterOptimizationResults プロパティから、または関数 bestPoint を使用して、最適な点を取得できます。

x = Mdl_univariate.HyperparameterOptimizationResults.XAtMinEstimatedObjective
x=1×3 table
    InitialLearnRateForPredictors    MaxNumSplitsPerPredictor    NumTreesPerPredictor
    _____________________________    ________________________    ____________________

              0.063687                          1                         61         

bestPoint(Mdl_univariate.HyperparameterOptimizationResults)
ans=1×3 table
    InitialLearnRateForPredictors    MaxNumSplitsPerPredictor    NumTreesPerPredictor
    _____________________________    ________________________    ____________________

              0.063687                          1                         61         

最適化プロセスの詳細については、OptimizeHyperparameters を使用した GAM の最適化を参照してください。

二変量の GAM に学習させます。名前と値の引数 InteractionsInitialLearnRateForInteractionsMaxNumSplitsPerInteraction および NumTreesPerInteraction の最適な値を fitrgam で求めるため、OptimizeHyperparameters として 'all-bivariate' を指定します。x の一変量パラメーターの値を使用して、交互作用項に最適なパラメーター値が x の値に基づいて特定されるようにします。

Mdl = fitrgam(tbl_training,'SALEPRICE','FitStandardDeviation',true, ...
    'InitialLearnRateForPredictors',x.InitialLearnRateForPredictors, ...
    'MaxNumSplitsPerPredictor',x.MaxNumSplitsPerPredictor, ...
    'NumTreesPerPredictor',x.NumTreesPerPredictor, ...
    'OptimizeHyperparameters','all-bivariate', ...
    'HyperparameterOptimizationOptions',struct('AcquisitionFunctionName','expected-improvement-plus', ...
    'ShowPlots',false,'Verbose',0))
Mdl = 
  RegressionGAM
                       PredictorNames: {'BOROUGH'  'NEIGHBORHOOD'  'BUILDINGCLASSCATEGORY'  'RESIDENTIALUNITS'  'COMMERCIALUNITS'  'LANDSQUAREFEET'  'GROSSSQUAREFEET'  'YEARBUILT'  'SALEDATE'}
                         ResponseName: 'SALEPRICE'
                CategoricalPredictors: [2 3]
                    ResponseTransform: 'none'
                            Intercept: 5.1679e+05
                         Interactions: [3×2 double]
               IsStandardDeviationFit: 1
                      NumObservations: 900
    HyperparameterOptimizationResults: [1×1 BayesianOptimization]


  Properties, Methods

最適な二変量のハイパーパラメーターを表示します。

Mdl.HyperparameterOptimizationResults.XAtMinEstimatedObjective
ans=1×4 table
    Interactions    InitialLearnRateForInteractions    MaxNumSplitsPerInteraction    NumTreesPerInteraction
    ____________    _______________________________    __________________________    ______________________

         3                     0.0010182                           21                         302          

Mdl のモデル表示には、モデルのプロパティの一部のみが表示されます。モデルのプロパティの完全な一覧を表示するには、ワークスペースで変数名 Mdl をダブルクリックします。Mdl の変数エディターが開きます。あるいは、コマンド ウィンドウでドット表記を使用してプロパティを表示できます。たとえば、ReasonForTermination プロパティを表示します。

Mdl.ReasonForTermination
ans = struct with fields:
      PredictorTrees: 'Terminated after training the requested number of trees.'
    InteractionTrees: 'Terminated after training the requested number of trees.'

ReasonForTermination プロパティを使用して、各線形項および各交互作用項に対して、学習済みモデルに指定した数の木が含まれているかどうかを確認できます。

Mdl の交互作用項を表示します。

Mdl.Interactions
ans = 3×2

     3     6
     4     6
     5     8

Interactions の各行は 1 つの交互作用項を表し、交互作用項の予測子変数の列インデックスを格納します。Interactions プロパティを使用して、モデル内の交互作用項とそれらが fitrgam でモデルに追加された順序を確認できます。

予測子名を使用して Mdl の交互作用項を表示します。

Mdl.PredictorNames(Mdl.Interactions)
ans = 3×2 cell
    {'BUILDINGCLASSCATEGORY'}    {'LANDSQUAREFEET'}
    {'RESIDENTIALUNITS'     }    {'LANDSQUAREFEET'}
    {'COMMERCIALUNITS'      }    {'YEARBUILT'     }

新しい観測値での予測性能の評価

テスト標本 tbl_test とオブジェクト関数 predict および loss を使用して、学習済みモデルの性能を評価します。これらの関数をもつ完全またはコンパクトなモデルを使用できます。

  • predict — 応答を予測

  • loss — 回帰損失 (既定では平均二乗誤差) を計算

学習データ セットの性能を評価するには、再代入オブジェクト関数resubPredictおよびresubLossを使用します。これらの関数を使用するには、学習データを含む完全なモデルを使用しなければなりません。

コンパクトなモデルを作成して、学習済みモデルのサイズを縮小します。

CMdl = compact(Mdl);
whos('Mdl','CMdl')
  Name      Size               Bytes  Class                                          Attributes

  CMdl      1x1             11975596  classreg.learning.regr.CompactRegressionGAM              
  Mdl       1x1             12170960  RegressionGAM                                            

線形項と交互作用項の両方を含めることにより得た結果と線形項のみを含めることにより得た結果を比較します。

テスト データ セット tbl_test について、応答を予測し、平均二乗誤差を計算します。

[yFit,ySD,yInt] = predict(CMdl,tbl_test);
L = loss(CMdl,tbl_test)
L = 1.2746e+11

学習済みモデルに交互作用項を含めずに予測された応答と誤差を調べます。

[yFit_nointeraction,ySD_nointeraction,yInt__nointeraction] = predict(CMdl,tbl_test,'IncludeInteractions',false);
L_nointeractions = loss(CMdl,tbl_test,'IncludeInteractions',false)
L_nointeractions = 1.2531e+11

交互作用項が含まれていない場合の方がテスト データ セットの誤差が小さくなっています。

並べ替えられた真の応答を予測応答および予測区間と一緒にプロットします。

yTrue = tbl_test.SALEPRICE;
[sortedYTrue,I] = sort(yTrue);

figure
ax = nexttile;
plot(sortedYTrue,'o')
hold on
plot(yFit(I))
plot(yInt(I,1),'k:')
plot(yInt(I,2),'k:')
legend('True responses','Predicted responses', ...
    '95% Prediction interval limits','Location','best')
title('Linear and interaction terms')
hold off

nexttile
plot(sortedYTrue,'o')
hold on
plot(yFit_nointeraction(I))
plot(yInt__nointeraction(I,1),'k:')
plot(yInt__nointeraction(I,2),'k:')
ylim(ax.YLim)
title('Linear terms only')
hold off

2 つのプロットの予測区間は同じような幅になっています。

予測の解釈

関数plotLocalEffectsを使用して、最初のテスト観測値についての予測を解釈します。また、関数plotPartialDependenceを使用して、モデル内のいくつかの重要な項の部分依存プロットを作成します。

テスト データの最初の観測値に対する応答値を予測し、予測に対する CMdl 内の項の局所的効果をプロットします。プロットに切片項を含めるには、'IncludeIntercept',true を指定します。

yFit = predict(CMdl,tbl_test(1,:))
yFit = 5.3526e+05
figure
plotLocalEffects(CMdl,tbl_test(1,:),'IncludeIntercept',true)

関数 predict で、最初の観測値 tbl_test(1,:) の売価を返します。関数 plotLocalEffects で、予測に対する CMdl 内の項の局所的効果を示す横棒グラフを作成します。局所的効果の各値は、tbl_test(1,:) の予測売価への各項の寄与を示します。

BUILDINGCLASSCATEGORY の部分従属の値を計算し、並べ替えられた値をプロットします。学習データ セットとテスト データ セットの両方を指定して、両方のセットの部分依存の値を計算します。

[pd,x] = partialDependence(CMdl,'BUILDINGCLASSCATEGORY',[tbl_training; tbl_test]);
[pd_sorted,I] = sort(pd);
x_sorted = x(I);
x_sorted = reordercats(x_sorted,I);
figure
plot(x_sorted,pd_sorted,'o:')
xlabel('BUILDINGCLASSCATEGORY')
ylabel('SALEPRICE')
title('Patial Dependence Plot')

プロットされたラインは、学習済みモデルにおける予測子 BUILDINGCLASSCATEGORY と応答 SALEPRICE の間の平均化された部分関係性を表します。

テスト データ セットを使用して、項 RESIDENTIALUNITSLANDSQUAREFEET の部分依存プロットを作成します。

figure
plotPartialDependence(CMdl,["RESIDENTIALUNITS","LANDSQUAREFEET"],tbl_test)

x 軸 (RESIDENTIALUNITS) と y 軸 (LANDSQUAREFEET) 小目盛りは、指定されたデータ内の一意の予測子値を表します。予測子値にはいくつかの外れ値があり、RESIDENTIALUNITSLANDSQUAREFEET の値のほとんどはそれぞれ 5 未満と 5000 未満です。プロットは、RESIDENTIALUNITS の値が 5 を超えると、SALEPRICE の値が大きくは変わらないことを示しています。

参考

| | | | | |

関連するトピック