データの分割

　最初にinsuranceデータフレームをトレーニング用とテスト用に分割します。線形回帰用に追加した非線形項（age2とbmi30）は外しておきます。なお、分割比率はデフォルトの\(3:1\)です。

ins_train

## # A tibble: 1,004 x 7
##      age sex      bmi children smoker region    expenses
##    <int> <fct>  <dbl>    <int> <fct>  <fct>        <dbl>
##  1    19 female  27.9        0 yes    southwest   16885.
##  2    18 male    33.8        1 no     southeast    1726.
##  3    28 male    33          3 no     southeast    4449.
##  4    32 male    28.9        0 no     northwest    3867.
##  5    46 female  33.4        1 no     southeast    8241.
##  6    37 female  27.7        3 no     northwest    7282.
##  7    37 male    29.8        2 no     northeast    6406.
##  8    60 female  25.8        0 no     northwest   28923.
##  9    62 female  26.3        0 yes    southeast   27809.
## 10    23 male    34.4        0 no     southwest    1827.
## # ... with 994 more rows

ins_test

## # A tibble: 334 x 7
##      age sex      bmi children smoker region    expenses
##    <int> <fct>  <dbl>    <int> <fct>  <fct>        <dbl>
##  1    33 male    22.7        0 no     northwest   21984.
##  2    31 female  25.7        0 no     southeast    3757.
##  3    25 male    26.2        0 no     northeast    2721.
##  4    19 male    24.6        1 no     southwest    1837.
##  5    52 female  30.8        1 no     northeast   10797.
##  6    18 male    34.1        0 no     southeast    1137.
##  7    55 female  32.8        2 no     northwest   12269.
##  8    63 male    28.3        0 no     northwest   13770.
##  9    19 male    20.4        0 no     northwest    1625.
## 10    26 male    20.8        0 no     southwest    2302.
## # ... with 324 more rows

モデルの作成

トレーニング用データを用いてモデルを作成します。

m_ins <- Cubist::cubist(ins_train[, -7], ins_train[, 7])
m_ins

## 
## Call:
## cubist.default(x = ins_train[, -7], y = ins_train[, 7])
## 
## Number of samples: 1004 
## Number of predictors: 6 
## 
## Number of committees: 1 
## Number of rules: 4

　モデルのルールは以下の通りです。

m_ins %>% 
  summary()

## 
## Call:
## cubist.default(x = ins_train[, -7], y = ins_train[, 7])
## 
## 
## Cubist [Release 2.07 GPL Edition]  Wed Nov 20 13:03:03 2019
## ---------------------------------
## 
##     Target attribute `outcome'
## 
## Read 1004 cases (7 attributes) from undefined.data
## 
## Model:
## 
##   Rule 1: [253 cases, mean 4366.201, range 1135.94 to 27724.29, est err 2023.477]
## 
##     if
##  age <= 29
##  smoker = no
##     then
##  outcome = -2678.36 + 214 age + 820 children
## 
##   Rule 2: [549 cases, mean 10118.000, range 3260.2 to 36910.61, est err 1510.557]
## 
##     if
##  age > 29
##  smoker = no
##     then
##  outcome = -6384.48 + 315 age + 506 children
## 
##   Rule 3: [97 cases, mean 21791.002, range 12829.46 to 38245.59, est err 1991.133]
## 
##     if
##  bmi <= 30.1
##  smoker = yes
##     then
##  outcome = 1722.86 + 247 age + 374 bmi + 125 children
## 
##   Rule 4: [105 cases, mean 41641.012, range 32548.34 to 62592.87, est err 1462.996]
## 
##     if
##  bmi > 30.1
##  smoker = yes
##     then
##  outcome = 14924.378 + 271 age + 425 bmi + 210 children
## 
## 
## Evaluation on training data (1004 cases):
## 
##     Average  |error|           2570.024
##     Relative |error|               0.28
##     Correlation coefficient        0.88
## 
## 
##  Attribute usage:
##    Conds  Model
## 
##    100%           smoker
##     80%   100%    age
##     20%    20%    bmi
##           100%    children
## 
## 
## Time: 0.0 secs

Model:

  Rule 1: [253 cases, mean 4366.201, range 1135.94 to 27724.29, est err 2023.477]

    if
    age <= 29
    smoker = no
    then
    outcome = -2678.36 + 214 age + 820 children

Evaluation on training data (1004 cases):

    Average  |error|           2570.024
    Relative |error|               0.28
    Correlation coefficient        0.88

    Attribute usage:
      Conds  Model

      100%           smoker
       80%   100%    age
       20%    20%    bmi
             100%    children

モデルを用いた予測

　作成したモデルにテスト用データを適用してモデル精度を確認します。線形回帰モデルと比較するために線形回帰モデルのときと同じプロットを行います。

p_ins <- m_ins %>% 
  predict(ins_test)

ins_test %>% 
  dplyr::bind_cols(pred = p_ins) %>% 
  ggplot2::ggplot(ggplot2::aes(x = expenses, y = pred)) + 
    ggplot2::geom_abline(slope = 1, colour = "red", linetype = "dotted") + 
    ggplot2::geom_point() + 
    ggplot2::labs(title = "予測値-観測値プロット（モデル木）",
                  x = "観測値", y = "予測値")

ins_test %>% 
  dplyr::bind_cols(pred = p_ins) %>% 
  dplyr::mutate(resid = expenses - pred) %>% 
  ggplot2::ggplot(ggplot2::aes(x = pred, y = resid)) + 
    ggplot2::geom_hline(yintercept = 0, colour = "red", linetype = "dotted") +
    ggplot2::geom_point() + ggplot2::ylim(-20000, 30000) +
    ggplot2::labs(title = "残差プロット（モデル木）", x = "予測値", y = "残差")

ins_test %>% 
  dplyr::bind_cols(pred = p_ins) %>% 
  with(cor(expenses, pred))

## [1] 0.9296292

　相関係数もかなり高く、まずまずのモデルになっていると言えます。

全データに対する予測

　線形回帰モデルと比較を行うために全データに対する予測を行います。こちらでも線形回帰用に追加した非線形項（age2とbmi30）を外します。

p_ins <- insurance %>% 
  dplyr::select(-age2, -bmi30) %>% 
  predict(m_ins, .)

gg_cubist <- insurance %>% 
  dplyr::select(-age2, -bmi30) %>% 
  dplyr::bind_cols(pred = p_ins) %>% 
  ggplot2::ggplot(ggplot2::aes(x = expenses, y = pred)) + 
    ggplot2::geom_abline(slope = 1, colour = "red", linetype = "dotted") + 
    ggplot2::geom_point() + 
    ggplot2::labs(title = "予測値-観測値プロット（モデル木）",
                  x = "観測値", y = "予測値")
gg_cubist

gg_cubist_resid <- insurance %>% 
  dplyr::select(-age2, -bmi30) %>% 
  dplyr::bind_cols(pred = p_ins) %>% 
  dplyr::mutate(resid = expenses - pred) %>% 
  ggplot2::ggplot(ggplot2::aes(x = pred, y = resid)) + 
    ggplot2::geom_hline(yintercept = 0, colour = "red", linetype = "dotted") +
    ggplot2::geom_point() + ggplot2::ylim(-20000, 30000) +
    ggplot2::labs(title = "残差プロット（モデル木）", x = "予測値", y = "残差")
gg_cubist_resid