boston

Bostonデータの分析

全変数を用いた線形回帰

# 必要なライブラリを読み込む
library(tidyverse)

Warning: package 'ggplot2' was built under R version 4.3.1

Warning: package 'tidyr' was built under R version 4.3.1

Warning: package 'readr' was built under R version 4.3.1

Warning: package 'dplyr' was built under R version 4.3.1

Warning: package 'stringr' was built under R version 4.3.1

Warning: package 'lubridate' was built under R version 4.3.1

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.0     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# CSVファイルを読み込む
boston_data <- read.csv("boston.csv")

# 線形回帰モデルを作成
# ここでは、MEDV（住宅価格の中央値）を目的変数、他のすべての変数を説明変数として使用
lm_model <- lm(MEDV ~ ., data = boston_data)

# モデルの要約を表示
summary(lm_model)

Call:
lm(formula = MEDV ~ ., data = boston_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.595  -2.730  -0.518   1.777  26.199 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
CRIM        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
ZN           4.642e-02  1.373e-02   3.382 0.000778 ***
INDUS        2.056e-02  6.150e-02   0.334 0.738288    
CHAS         2.687e+00  8.616e-01   3.118 0.001925 ** 
NOX         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
RM           3.810e+00  4.179e-01   9.116  < 2e-16 ***
AGE          6.922e-04  1.321e-02   0.052 0.958229    
DIS         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
RAD          3.060e-01  6.635e-02   4.613 5.07e-06 ***
TAX         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
PTRATIO     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
B            9.312e-03  2.686e-03   3.467 0.000573 ***
LSTAT       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406,    Adjusted R-squared:  0.7338 
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

1. ライブラリの読み込み: tidyverseパッケージを使用してデータ操作と解析を行います。

2. データの読み込み: read_csv関数を使用して、CSVファイルからデータを読み込みます。ファイルパスは適宜調整してください。

3. 線形回帰モデルの作成: lm関数を使用して、線形回帰モデルを作成します。MEDV（住宅価格の中央値）を目的変数とし、他のすべての変数を説明変数として指定します。

4. モデルの要約表示: summary関数を使用して、線形回帰モデルの結果を表示します。これには、係数、標準誤差、t値、p値、決定係数などの統計情報が含まれます。

ステップワイズ法による変数選択

# 必要なライブラリを読み込む
library(tidyverse)
library(MASS)  # ステップワイズ法に使用

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select

# CSVファイルを読み込む
boston_data <- read_csv("boston.csv")

Rows: 506 Columns: 14

── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (14): CRIM, ZN, INDUS, CHAS, NOX, RM, AGE, DIS, RAD, TAX, PTRATIO, B, LS...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# フルモデルを作成
full_model <- lm(MEDV ~ ., data = boston_data)

# ステップワイズ法による変数選択を実行
stepwise_model <- stepAIC(full_model, direction = "both")

Start:  AIC=1589.64
MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + 
    TAX + PTRATIO + B + LSTAT

          Df Sum of Sq   RSS    AIC
- AGE      1      0.06 11079 1587.7
- INDUS    1      2.52 11081 1587.8
<none>                 11079 1589.6
- CHAS     1    218.97 11298 1597.5
- TAX      1    242.26 11321 1598.6
- CRIM     1    243.22 11322 1598.6
- ZN       1    257.49 11336 1599.3
- B        1    270.63 11349 1599.8
- RAD      1    479.15 11558 1609.1
- NOX      1    487.16 11566 1609.4
- PTRATIO  1   1194.23 12273 1639.4
- DIS      1   1232.41 12311 1641.0
- RM       1   1871.32 12950 1666.6
- LSTAT    1   2410.84 13490 1687.3

Step:  AIC=1587.65
MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + DIS + RAD + TAX + 
    PTRATIO + B + LSTAT

          Df Sum of Sq   RSS    AIC
- INDUS    1      2.52 11081 1585.8
<none>                 11079 1587.7
+ AGE      1      0.06 11079 1589.6
- CHAS     1    219.91 11299 1595.6
- TAX      1    242.24 11321 1596.6
- CRIM     1    243.20 11322 1596.6
- ZN       1    260.32 11339 1597.4
- B        1    272.26 11351 1597.9
- RAD      1    481.09 11560 1607.2
- NOX      1    520.87 11600 1608.9
- PTRATIO  1   1200.23 12279 1637.7
- DIS      1   1352.26 12431 1643.9
- RM       1   1959.55 13038 1668.0
- LSTAT    1   2718.88 13798 1696.7

Step:  AIC=1585.76
MEDV ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + 
    B + LSTAT

          Df Sum of Sq   RSS    AIC
<none>                 11081 1585.8
+ INDUS    1      2.52 11079 1587.7
+ AGE      1      0.06 11081 1587.8
- CHAS     1    227.21 11309 1594.0
- CRIM     1    245.37 11327 1594.8
- ZN       1    257.82 11339 1595.4
- B        1    270.82 11352 1596.0
- TAX      1    273.62 11355 1596.1
- RAD      1    500.92 11582 1606.1
- NOX      1    541.91 11623 1607.9
- PTRATIO  1   1206.45 12288 1636.0
- DIS      1   1448.94 12530 1645.9
- RM       1   1963.66 13045 1666.3
- LSTAT    1   2723.48 13805 1695.0

# 選択されたモデルの要約を表示
summary(stepwise_model)

Call:
lm(formula = MEDV ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + 
    TAX + PTRATIO + B + LSTAT, data = boston_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.5984  -2.7386  -0.5046   1.7273  26.2373 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
CRIM         -0.108413   0.032779  -3.307 0.001010 ** 
ZN            0.045845   0.013523   3.390 0.000754 ***
CHAS          2.718716   0.854240   3.183 0.001551 ** 
NOX         -17.376023   3.535243  -4.915 1.21e-06 ***
RM            3.801579   0.406316   9.356  < 2e-16 ***
DIS          -1.492711   0.185731  -8.037 6.84e-15 ***
RAD           0.299608   0.063402   4.726 3.00e-06 ***
TAX          -0.011778   0.003372  -3.493 0.000521 ***
PTRATIO      -0.946525   0.129066  -7.334 9.24e-13 ***
B             0.009291   0.002674   3.475 0.000557 ***
LSTAT        -0.522553   0.047424 -11.019  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.736 on 494 degrees of freedom
Multiple R-squared:  0.7406,    Adjusted R-squared:  0.7348 
F-statistic: 128.2 on 11 and 494 DF,  p-value: < 2.2e-16

1. ライブラリの読み込み:

   • tidyverse: データ操作と解析を行うためのパッケージ。

   • MASS: ステップワイズ法を実行するための関数 stepAIC が含まれるパッケージ。

2. データの読み込み:

   • read_csv 関数を使用して、CSVファイルからデータを読み込みます。ファイルパスは適宜調整してください。

3. フルモデルの作成:

   • lm 関数を使用して、すべての説明変数を含む線形回帰モデル（フルモデル）を作成します。

4. ステップワイズ法による変数選択の実行:

   • stepAIC 関数を使用して、ステップワイズ法を実行します。direction = "both" は、前進選択と後退削除の両方を行うことを意味します。

   • stepAIC は、Akaike情報量基準 (AIC) を使用してモデルの適合度を評価し、最適な変数の組み合わせを選択します。

5. 選択されたモデルの要約を表示:

   • summary 関数を使用して、選択されたモデルの結果を表示します。これには、選択された変数、係数、標準誤差、t値、p値、決定係数などの統計情報が含まれます。

Lasso回帰

# 必要なライブラリを読み込む
library(tidyverse)
library(glmnet)

Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

Loaded glmnet 4.1-8

# CSVファイルを読み込む
boston_data <- read_csv("boston.csv")

Rows: 506 Columns: 14

── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (14): CRIM, ZN, INDUS, CHAS, NOX, RM, AGE, DIS, RAD, TAX, PTRATIO, B, LS...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# 説明変数（X）と目的変数（y）を設定
X <- as.matrix(boston_data %>% dplyr::select(-MEDV))
y <- boston_data$MEDV

# Lasso回帰のモデルを構築
lasso_model <- cv.glmnet(X, y, alpha = 1)

# 最適なラムダ値（正則化パラメータ）を取得
best_lambda <- lasso_model$lambda.min

# 最適なラムダ値で再度モデルを構築
final_model <- glmnet(X, y, alpha = 1, lambda = best_lambda)

# モデルの係数を表示
coef(final_model)

14 x 1 sparse Matrix of class "dgCMatrix"
                       s0
(Intercept)  34.767597483
CRIM         -0.100135466
ZN            0.042313328
INDUS         .          
CHAS          2.689677181
NOX         -16.487728982
RM            3.854232596
AGE           .          
DIS          -1.413899564
RAD           0.261086342
TAX          -0.010183187
PTRATIO      -0.932475665
B             0.009069592
LSTAT        -0.522594035

1. 必要なライブラリを読み込む

まず、tidyverseとglmnetの2つのライブラリを読み込みます。tidyverseはデータ操作や解析を行うためのパッケージ群で、特にデータの読み込みと前処理に使用されます。一方、glmnetはLasso回帰やRidge回帰を実行するためのパッケージで、正則化を用いた回帰分析を行うために使用されます。

2. CSVファイルを読み込む

次に、read_csv関数を使ってCSVファイルからデータを読み込みます。この関数は指定されたファイルパスにあるCSVファイルを読み込み、データフレーム形式で格納します。このデータフレームをboston_dataという変数に保存します。

3. 説明変数（X）と目的変数（y）を設定

データフレームから目的変数（住宅価格の中央値、MEDV列）を除いた全ての列を説明変数として選択し、行列形式に変換します。これをXという変数に格納します。また、MEDV列自体を目的変数としてベクトル形式でyという変数に保存します。

4. Lasso回帰のモデルを構築

cv.glmnet関数を使用して、クロスバリデーションを行いながらLasso回帰モデルを構築します。ここでalpha = 1を指定することで、Lasso回帰が実行されます。クロスバリデーションにより、データセットを複数の部分に分割してモデルを評価し、最適な正則化パラメータ（ラムダ）を選択します。

5. 最適なラムダ値（正則化パラメータ）を取得

cv.glmnetの結果から、最適なラムダ値を取得します。最適なラムダ値は、クロスバリデーションの結果に基づいて決定され、モデルの過剰適合を防ぎつつ予測性能を最適化するために使用されます。

6. 最適なラムダ値で再度モデルを構築

先ほど取得した最適なラムダ値を用いて、再度Lasso回帰モデルを構築します。ここでもalpha = 1を指定してLasso回帰を行います。このモデルは、最適な正則化パラメータを適用することで、不要な変数の係数をゼロにし、重要な変数だけを選択します。

7. モデルの係数を表示

最終的に構築されたモデルの係数を表示します。Lasso回帰では、不要な変数の係数がゼロになるため、どの変数が選択されたか、どの変数が選択されなかったかを確認することができます。このステップにより、モデルに含まれる重要な変数とその影響度を評価することができます。

以上が、Rを使ったLasso回帰による変数選択のプロセスの詳細な説明です。これにより、データの中から重要な説明変数を特定し、過剰適合を防ぎつつ、予測性能の高いモデルを構築することができます。

StepwiseとLasso回帰の比較

# 必要なライブラリを読み込む
library(tidyverse)
library(MASS)
library(glmnet)
library(caret)

Loading required package: lattice

Attaching package: 'caret'

The following object is masked from 'package:purrr':

    lift

# データを読み込む
boston_data <- read_csv("boston.csv")

Rows: 506 Columns: 14

── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (14): CRIM, ZN, INDUS, CHAS, NOX, RM, AGE, DIS, RAD, TAX, PTRATIO, B, LS...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# ステップワイズ法のモデルを作成
full_model <- lm(MEDV ~ ., data = boston_data)
stepwise_model <- stepAIC(full_model, direction = "both")

Start:  AIC=1589.64
MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + 
    TAX + PTRATIO + B + LSTAT

          Df Sum of Sq   RSS    AIC
- AGE      1      0.06 11079 1587.7
- INDUS    1      2.52 11081 1587.8
<none>                 11079 1589.6
- CHAS     1    218.97 11298 1597.5
- TAX      1    242.26 11321 1598.6
- CRIM     1    243.22 11322 1598.6
- ZN       1    257.49 11336 1599.3
- B        1    270.63 11349 1599.8
- RAD      1    479.15 11558 1609.1
- NOX      1    487.16 11566 1609.4
- PTRATIO  1   1194.23 12273 1639.4
- DIS      1   1232.41 12311 1641.0
- RM       1   1871.32 12950 1666.6
- LSTAT    1   2410.84 13490 1687.3

Step:  AIC=1587.65
MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + DIS + RAD + TAX + 
    PTRATIO + B + LSTAT

          Df Sum of Sq   RSS    AIC
- INDUS    1      2.52 11081 1585.8
<none>                 11079 1587.7
+ AGE      1      0.06 11079 1589.6
- CHAS     1    219.91 11299 1595.6
- TAX      1    242.24 11321 1596.6
- CRIM     1    243.20 11322 1596.6
- ZN       1    260.32 11339 1597.4
- B        1    272.26 11351 1597.9
- RAD      1    481.09 11560 1607.2
- NOX      1    520.87 11600 1608.9
- PTRATIO  1   1200.23 12279 1637.7
- DIS      1   1352.26 12431 1643.9
- RM       1   1959.55 13038 1668.0
- LSTAT    1   2718.88 13798 1696.7

Step:  AIC=1585.76
MEDV ~ CRIM + ZN + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + 
    B + LSTAT

          Df Sum of Sq   RSS    AIC
<none>                 11081 1585.8
+ INDUS    1      2.52 11079 1587.7
+ AGE      1      0.06 11081 1587.8
- CHAS     1    227.21 11309 1594.0
- CRIM     1    245.37 11327 1594.8
- ZN       1    257.82 11339 1595.4
- B        1    270.82 11352 1596.0
- TAX      1    273.62 11355 1596.1
- RAD      1    500.92 11582 1606.1
- NOX      1    541.91 11623 1607.9
- PTRATIO  1   1206.45 12288 1636.0
- DIS      1   1448.94 12530 1645.9
- RM       1   1963.66 13045 1666.3
- LSTAT    1   2723.48 13805 1695.0

# Lasso回帰のモデルを作成
X <- as.matrix(boston_data %>% dplyr::select(-MEDV))
y <- boston_data$MEDV
lasso_model <- cv.glmnet(X, y, alpha = 1)
best_lambda <- lasso_model$lambda.min
final_lasso_model <- glmnet(X, y, alpha = 1, lambda = best_lambda)

# ステップワイズ法で選択された変数
stepwise_selected_vars <- names(coef(stepwise_model))[-1]

# Lasso回帰で選択された変数
lasso_coef <- coef(final_lasso_model)
lasso_selected_vars <- rownames(lasso_coef)[lasso_coef[,1] != 0]

# 選択された変数を表示
cat("ステップワイズ法で選択された変数:\n")

ステップワイズ法で選択された変数:

print(stepwise_selected_vars)

 [1] "CRIM"    "ZN"      "CHAS"    "NOX"     "RM"      "DIS"     "RAD"    
 [8] "TAX"     "PTRATIO" "B"       "LSTAT"  

cat("Lasso回帰で選択された変数:\n")

Lasso回帰で選択された変数:

print(lasso_selected_vars)

 [1] "(Intercept)" "CRIM"        "ZN"          "CHAS"        "NOX"        
 [6] "RM"          "DIS"         "RAD"         "TAX"         "PTRATIO"    
[11] "B"           "LSTAT"      

# モデル性能の比較
# データを訓練セットとテストセットに分割
set.seed(123)
trainIndex <- createDataPartition(boston_data$MEDV, p = 0.8, list = FALSE)
train_data <- boston_data[trainIndex, ]
test_data <- boston_data[-trainIndex, ]

# ステップワイズ法のモデルで予測
stepwise_model_final <- lm(MEDV ~ ., data = train_data)
stepwise_predictions <- predict(stepwise_model_final, newdata = test_data)

# Lasso回帰のモデルで予測
X_train <- as.matrix(train_data %>% dplyr::select(-MEDV))
y_train <- train_data$MEDV
X_test <- as.matrix(test_data %>% dplyr::select(-MEDV))
lasso_model_final <- glmnet(X_train, y_train, alpha = 1, lambda = best_lambda)
lasso_predictions <- predict(lasso_model_final, newx = X_test)

# MSEの計算
stepwise_mse <- mean((test_data$MEDV - stepwise_predictions)^2)
lasso_mse <- mean((test_data$MEDV - lasso_predictions)^2)

# 決定係数（R²）の計算
stepwise_r2 <- 1 - sum((test_data$MEDV - stepwise_predictions)^2) / sum((test_data$MEDV - mean(test_data$MEDV))^2)
lasso_r2 <- 1 - sum((test_data$MEDV - lasso_predictions)^2) / sum((test_data$MEDV - mean(test_data$MEDV))^2)

# 結果を表示
cat("ステップワイズ法のMSE:", stepwise_mse, "\n")

ステップワイズ法のMSE: 21.05844 

cat("Lasso回帰のMSE:", lasso_mse, "\n")

Lasso回帰のMSE: 21.03565 

cat("ステップワイズ法のR²:", stepwise_r2, "\n")

ステップワイズ法のR²: 0.7572919 

cat("Lasso回帰のR²:", lasso_r2, "\n")

Lasso回帰のR²: 0.7575546 

1. 必要なライブラリを読み込む

まず、以下のライブラリを読み込みます：

• tidyverse: データ操作と解析を行うためのパッケージ群。

• MASS: ステップワイズ回帰を行うためのstepAIC関数を含むパッケージ。

• glmnet: Lasso回帰やRidge回帰を実行するためのパッケージ。

• caret: データの訓練セットとテストセットへの分割など、機械学習用のユーティリティ関数を提供。

2. データを読み込む

read_csv関数を使ってCSVファイルからデータを読み込み、boston_dataというデータフレームに格納します。

3. ステップワイズ法のモデルを作成

まず、全ての変数を用いたフルモデルをlm関数で作成します。その後、stepAIC関数を使用してステップワイズ法による変数選択を行い、最適なモデルをstepwise_modelに保存します。この過程では前進選択と後退削除の両方を行います。

4. Lasso回帰のモデルを作成

説明変数を行列形式に変換し、目的変数をベクトル形式に設定します。cv.glmnet関数を用いてクロスバリデーションを行いながらLasso回帰モデルを構築し、最適なラムダ値を選択します。その後、最適なラムダ値を用いてLasso回帰モデルを再構築します。

5. ステップワイズ法で選択された変数

coef関数を使用して、ステップワイズ法で選択された変数のリストを取得します。names関数を使って変数名を取得し、先頭のインターセプトを除いたリストを表示します。

6. Lasso回帰で選択された変数

coef関数を使ってLasso回帰モデルの係数を取得し、非ゼロの係数を持つ変数を選択します。rownames関数を使用して変数名を取得し、非ゼロの係数を持つ変数のリストを表示します。

7. 選択された変数を表示

catとprint関数を使って、ステップワイズ法とLasso回帰で選択された変数をそれぞれ表示します。

8. モデル性能の比較

データを訓練セットとテストセットに分割します。set.seed関数でランダムシードを設定し、createDataPartition関数を使ってデータを分割します。

9. ステップワイズ法のモデルで予測

訓練セットを用いてステップワイズ法の最終モデルを構築し、テストセットで予測を行います。

10. Lasso回帰のモデルで予測

訓練セットを行列形式に変換し、再度Lasso回帰モデルを構築します。その後、テストセットで予測を行います。

11. MSEの計算

ステップワイズ法とLasso回帰のそれぞれについて、予測値と実際の値の二乗誤差の平均（MSE）を計算します。

12. 決定係数（R²）の計算

ステップワイズ法とLasso回帰のそれぞれについて、決定係数（R²）を計算します。これはモデルの予測精度を評価するための指標です。

13. 結果を表示

cat関数を使って、ステップワイズ法とLasso回帰のそれぞれのMSEとR²の結果を表示します。

この一連の手順により、ステップワイズ法とLasso回帰のモデル選択結果と予測性能を比較することができます。