# Lasso(ラッソ)回帰・Ridge(リッジ)回帰を線形回帰と比較して実装してみよう！

> 当サイト【スタビジ】の本記事では、線形回帰・Lasso(ラッソ)回帰・Ridge(リッジ)回帰についてRでの実装も行いながらまとめていきます。高次元データを分析する際はぜひLasso(ラッソ)回帰・Ridge(リッジ)回帰を利用してみてくださいね！

- URL: https://toukei-lab.com/%e7%b7%9a%e5%bd%a2%e5%9b%9e%e5%b8%b0%e3%83%bblasso%e3%83%bbridge%e3%82%92%e5%ae%9f%e8%a3%85%e3%81%97%e3%81%a6%e3%81%bf%e3%82%88%e3%81%86%ef%bc%81
- 公開: 2022.11.29 / 更新: 2024.10.31
- サイト: スタビジ（https://toukei-lab.com）

:::talk{who="umatan"}
本記事では、線形回帰・Lasso(ラッソ)回帰・Ridge(リッジ)回帰についてRでの実装も行いながらまとめていきます。高次元データを分析する際はぜひラッソ回帰・リッジ回帰を利用してみてくださいね！
:::

変数が多い高次元データでは、通常の[回帰分析](/%e5%9b%9e%e5%b8%b0%e5%88%86%e6%9e%90)が上手くいかない/そもそも解析不能なんてことが往々に起こります。

そんな高次元データの不都合を解決するのが「**Lasso(ラッソ)回帰**」や「**Ridge(リッジ)回帰**」！！

Lasso回帰・Ridge回帰は[**正則化**](/%e6%ad%a3%e5%89%87%e5%8c%96)や[**スパース推定**](/%e3%82%b9%e3%83%91%e3%83%bc%e3%82%b9%e6%8e%a8%e5%ae%9a)と関連性が深いので、そちらを先にご覧いただけると理解が深まると思います。

<https://toukei-lab.com/%e6%ad%a3%e5%89%87%e5%8c%96>

<https://toukei-lab.com/%e3%82%b9%e3%83%91%e3%83%bc%e3%82%b9%e6%8e%a8%e5%ae%9a>

:::talk{who="umatan"}
回帰分析にも色んな種類があるんだよねー！
:::

この記事では、高次元データに対して効果を発揮する「Lasso(ラッソ)回帰」と「Ridge(リッジ)回帰」について見ていきましょう！

式や概念の説明だといまいちどんなことが起きているのかわかりにくいので<mark>実際にR言語を用いて実装していきます！</mark>

## 線形回帰・Lasso回帰・Ridge回帰とは

![](/wp-content/uploads/2019/04/52ceff2ccf51b6465679ab4c976fc8b0_s.jpg)

今回比較するのは、**通常の線形回帰・Lasso回帰・Ridge回帰**です。\
簡単に、三つの手法について見ていきましょう！

### 通常の線形回帰

通常の線形回帰は残差${\bf y}-{\bf X}{\bf \beta}$が最小になるように$\beta$を推定します。\
つまり、次の式を最小化するように$\beta$を推定します。\
\begin{eqnarray\*}\
\\|{\bf y}-{\bf X}{\bf \beta}\\|\_2^2\\\\\
\end{eqnarray\*}

:::talk{who="robotan"}
線形回帰は、「データの間にある関係性を見つけて、**まっすぐな線**で予測しちゃおう！」という感じの手法だよ。データの変数（例えば身長と体重とか）の関係を直線で表して、この線を使って新しいデータの予測をするんだ。
:::

### Ridge回帰

![](/wp-content/uploads/2017/11/91dc209de29f987fd2a4d32f9a55cfeb-1024x443.png)

Ridge回帰は正則化手法であり、<mark>多重共線性などが起きていて通常の線形回帰ではパラメータが推定できないような場面でも罰則を与えることで推定できるようにしたものです。</mark>

:::talk{who="robotan"}
「多重共線性」についても確認しておこう！
:::

次の式を最小化するように$\beta$を推定します。\
\begin{eqnarray\*}\
\\|{\bf y}-{\bf X}{\bf \beta}\\|\_2^2+\lambda\\|{\bf \beta}\\|\_2^2\\\\\
\end{eqnarray\*}\
第２項は**正則化項**と呼ばれ、なるべくパラメータ$\beta$の絶対値を小さくするような働きをします。

ここで、$\lambda$は正則化パラメータであり、解析者が設定するかクロスバリデーションなどで決定する必要があります。

詳しく知りたい場合、「Ridge回帰」の論文を見るとよいでしょう！

> Proposed is an estimation procedure based on adding small positive quantities to the diagonal of **X**′**X**. Introduced is the ridge trace, a method for showing in two dimensions the effects of nonorthogonality.
>
> 引用：[Google-“Ridge Regression: Biased Estimation for Nonorthogonal Problems”](https://www.tandfonline.com/doi/abs/10.1080/00401706.2000.10485983?role=button\&needAccess=true\&journalCode=utch20)

### Lasso回帰

Lasso回帰も**正則化手法**です。Ridge回帰とは異なり<mark>パラメータのいくつかを0に推定できる</mark>という特徴があります。

これを<mark>**スパース推定**</mark>と呼びます。次の式を最小化するように$\beta$を推定します。\
\begin{eqnarray\*}\
\\|{\bf y}-{\bf X}{\bf \beta}\\|\_2^2+\lambda\\|{\bf \beta}\\|\_1\\\\\
\end{eqnarray\*}

今回は、線形回帰はglm()関数、LassoとRidgeはglmnetパッケージのglmnet()関数を用います。

詳しく知りたい場合、「Lasso回帰」の論文を見るとよいでしょう！

> We propose a new method for estimation in linear models. the ‘lasso’ minimizes the residual sum of squares subject to the sum of the absolute value of the coefficients being less than a constant. Because of the nature of this constraint it tends to produce some coefficients that are exactly 0 and hence gives interpretable models.
>
> 引用：[Google-“Regression Shrinkage and Selection via the Lasso”](https://www.ccs.neu.edu/home/eelhami/courses/EE290A/LASSO_Tibshirani.pdf)

## Lasso回帰・Ridge回帰用の人工データ生成

![](/wp-content/uploads/2019/06/PAK12_10many_TP_V-min-e1560939461601-1024x576.jpg)

今回は、実データではなく自分で発生させる人工データを扱います！\
今回のデータの生成モデルを次に示します。

\begin{eqnarray\*}\
y_i&=&1+{\bf x}\_i^T{\bf \beta}+\epsilon_i\\\\\
{\bf \beta}&=&(1,0,-1,0,0)^T\\\\\
{\bf x_i}&\sim& N({\bf \mu},{\bf \Sigma})\\\\\
{\bf \mu}&=&{\bf 0}\_5\\\\\
{\bf \Sigma}&=&{\bf I}\_5\\\\\
\epsilon_i&\sim& N(0,1^2)\\\\\
i&=&1,2,\cdots,n\\\\\
\end{eqnarray\*}

説明変数は、平均0・独立で発生させています。\
また、$\beta$の設定から変数2・変数4・変数5は目的変数に寄与していません。

## 線形回帰・Lasso回帰・Ridge回帰をRで実装

![](/wp-content/uploads/2018/09/goriIMGL9703_TP_V1-e1555747456954.jpg)

先ほど説明したデータで解析をしてみましょう！

サンプルサイズ$n$は1100とします。\
<mark>100個のデータを学習データ、残りの1000個のデータをテストデータ</mark>として実験します。

それぞれ３手法で推定された回帰係数と予測誤差(MSPE;Mean Squared Prediction Error)を計算してみましょう。

```text title="LinearRegression.txt"
library(MASS)
library(mvtnorm)
library(glmnet)

n_train <- 100
n_test <- 1000
p <- 5

β0 <- 1 ##定数項
β <- c(1,0,-1,0,0)
β <- c(β0,β)

μ <- numeric(p)
Σ <- diag(p)

σ <- 1 ##誤差標準偏差

X <- rmvnorm((n_train+n_test),μ,Σ)
X <- cbind(rep(1,(n_train+n_test)),X)
ε <- rnorm((n_train+n_test),0,σ)
y <- X%*%β+ε

data <- cbind(X,y)
data <- as.data.frame(data)
colnames(data) <- c("Int","x1","x2","x3","x4","x5","y")

##学習データとテストデータを作成
data_train <- data[1:n_train,]
data_test <- data[-(1:n_train),]

##glm(普通の線形回帰)
result_glm <- glm(formula=y~x1+x2+x3+x4+x5, data=data_train[,-1], family=gaussian(link="identity"))
##推定された回帰係数
as.matrix(result_glm$coefficients)

##glmnet(lasso)
result_lasso <- cv.glmnet(x=as.matrix(data_train[,c("x1","x2","x3","x4","x5")]), y=data_train$y, family="gaussian",alpha=1)
##推定された回帰係数
coefficients(result_lasso,s="lambda.min")

##glmnet(ridge)
result_ridge <- cv.glmnet(x=as.matrix(data_train[,c("x1","x2","x3","x4","x5")]), y=data_train$y, family="gaussian",alpha=0)
##推定された回帰係数
coefficients(result_ridge,s="lambda.min")

##予測誤差を算出
yhat_glm <- as.matrix(data_test[,colnames(data_test)!="y"])%*%as.matrix(result_glm$coefficients)
mean((data_test$y-yhat_glm)^2)

yhat_lasso <- as.matrix(data_test[,colnames(data_test)!="y"])%*%coefficients(result_lasso,s="lambda.min")
mean((data_test$y-yhat_lasso)^2)

yhat_ridge <- as.matrix(data_test[,colnames(data_test)!="y"])%*%coefficients(result_ridge,s="lambda.min")
mean((data_test$y-yhat_ridge)^2)
```

### 線形回帰・Lasso回帰・Ridge回帰の回帰係数

**■通常の線形回帰で推定された回帰係数**

(Intercept) 1.110727505\
x1 0.862992782\
x2 -0.024437489\
x3 -0.757605621\
x4 -0.098294387\
x5 -0.008577822

**■Lasso回帰で推定された回帰係数**

(Intercept) 1.0713394\
x1 0.7781268\
x2 .\
x3 -0.6778876\
x4 .\
x5 .

**■Ridge回帰で推定された回帰係数**

(Intercept) 1.08903319\
x1 0.81095492\
x2 -0.02288852\
x3 -0.71278554\
x4 -0.09335028\
x5 -0.01748969

(Intercept)は定数項を表しています。\
どうでしょうか、Lasso回帰は変数1と変数3以外は完全に0になっていることがわかります。\
つまり、<mark>自動的に変数選択ができました！これは嬉しいですね！</mark>

:::talk{who="umatan"}
「.」と書いてあるところは０と推定されているよ！
:::

### 回帰係数の2乗和を比較

また、推定された回帰係数の２乗和を計算すると次のようになりました。

**■通常の線形回帰**\
\> sum(as.matrix(result_glm\$coefficients)^2)\
2.562771

**■Lasso**\
\> sum(coefficients(result_lasso,s=”lambda.min”)^2)\
**2.212781**

**■Ridge**\
\> sum(coefficients(result_ridge,s=”lambda.min”)^2)\
**2.361248**

やはり、通常の線形回帰と比較して小さくなっています。

### MSPEを比較

最後に、MSPEを比較します。

※MSPEは、**Mean square prediction error**で予測値と観測値の差の二乗平均を取った値です。

**■通常の線形回帰**\
\> mean((data_test\$y-yhat_glm)^2)\
1.119285

**■Lasso**\
\> mean((data_test\$y-yhat_lasso)^2)\
1.171355

**■Ridge**\
\> mean((data_test\$y-yhat_ridge)^2)\
1.153019

最も良いのは通常の線形回帰になりました。

これはサンプルをたくさんあって、<mark>説明変数に相関もなく多重共線性が起きていないから当たり前ですね！</mark>

もっと、変数を増やしてみたりするとLassoやRidgeがMSPEで勝つこともあるでしょう。

ちなみに、このプログラムでは一回の計算なので結果の信頼性には欠けます。

そのため、<mark>続いてシミュレーションを何回か行って精度を検証していきましょう！</mark>

## 線形回帰・Lasso回帰・Ridge回帰の精度をシミュレーションで比較

![](/wp-content/uploads/2018/08/PAK85_MBAdesagyou20140312_TP_V-e1557018085482-1024x575.jpg)

どの手法が良いか・どのモデルが良いかを比較する為にはシミュレーションを繰り返して、<mark>指標の平均やばらつきをみて判断することが大切です。</mark>

たった一回の結果では本当に差があるのか偶然なのか信頼性に欠けるからです。

今回は<mark>線形回帰モデルの人工データを用いて線形回帰・Lasso・Ridgeの精度比較を行います。</mark>

### 実験の設定と手順

:::box
**・比較手法：通常の線形回帰・Lasso回帰・Ridge回帰**\
**・評価指標：MSPE(平均二乗予測誤差)**\
**・シミュレーション回数：100回**
:::

実験の流れは簡単にこんな感じです。

1. 人工データを生成
2. データを解析
3. 指標を算出して結果を格納
4. 1\~3を100回繰り返す
5. 結果の平均とばらつきで比較

### 人工データ生成

今回のデータの生成モデルを次に示します。

\begin{eqnarray\*}\
y_i&=&1+{\bf x}\_i^T{\bf \beta}+\epsilon_i\\\\\
{\bf \beta}&=&(1,-1,0,\cdots,0)^T\in {\bf R}^p\\\\\
{\bf x_i}&\sim& N({\bf \mu},{\bf \Sigma})\\\\\
{\bf \mu}&=&{\bf 0}\_p\\\\\
{\bf \Sigma}&=&{\bf I}\_p\\\\\
\epsilon_i&\sim& N(0,1^2)\\\\\
i&=&1,2,\cdots,n\\\\\
\end{eqnarray\*}

説明変数は$p$次元であり、平均0・独立で発生させています。

また、$\beta$はp次元ベクトルであり、変数1と変数2のみが目的変数へ寄与があります。

つまり、かなりスパースな設定になっています。

今回は、$p$を5,10,50と変化させて実験してみます。

### シミュレーション実験

シミュレーション実験をしてみましょう。

サンプルサイズ$n$は1100とします。

100個のデータを学習データ、残りの1000個のデータをテストデータとして実験します。

今回用いるコード（$p=50$の場合）を載せておきます。\
実行する前に入っていないパッケージはインストールしておいてください。

```text title="LinearRegressionSimulaion"
library(MASS)
library(mvtnorm)
library(glmnet)
library(ggplot2)
library(tidyr)

times <- 100 ##シミュレーション回数

n_train <- 100
n_test <- 1000
p <- 50 #変化させて実験

β0 <- 1 ##定数項
β <- c(1,-1,numeric(p-2))
β <- c(β0,β)

μ <- numeric(p)
Σ <- diag(p)

σ <- 1 ##誤差標準偏差

MSPE <- data.frame(glm=numeric(times)+NaN, lasso=numeric(times)+NaN, ridge=numeric(times)+NaN)

for(t in 1:times){
  X <- rmvnorm((n_train+n_test),μ,Σ)
  X <- cbind(rep(1,(n_train+n_test)),X)
  ε <- rnorm((n_train+n_test),0,σ)
  y <- X%*%β+ε
  
  data <- cbind(X,y)
  data <- as.data.frame(data)
  colnames(data) <- c("Int",paste0('x', 1:p),"y")

  ##学習データとテストデータを作成
  data_train <- data[1:n_train,]
  data_test <- data[-(1:n_train),]
  
  ##glm(普通の線形回帰)
  result_glm <- glm(formula=y~., data=data_train[,colnames(data_train)!="Int"], family=gaussian(link="identity"))
  
  ##glmnet(lasso)
  result_lasso <- cv.glmnet(x=as.matrix(data_train[,!(colnames(data_train) %in% c("Int","y"))]), 
                            y=data_train$y, family="gaussian",alpha=1)
  
  ##glmnet(ridge)
  result_ridge <- cv.glmnet(x=as.matrix(data_train[,!(colnames(data_train) %in% c("Int","y"))]), 
                            y=data_train$y, family="gaussian",alpha=0)
  
  ##予測誤差を算出
  yhat_glm <- as.matrix(data_test[,colnames(data_test)!="y"])%*%as.matrix(result_glm$coefficients)
  MSPE$glm[t] <- mean((data_test$y-yhat_glm)^2)
  
  yhat_lasso <- as.matrix(data_test[,colnames(data_test)!="y"])%*%coefficients(result_lasso,s="lambda.min")
  MSPE$lasso[t] <- mean((data_test$y-yhat_lasso)^2)
  
  yhat_ridge <- as.matrix(data_test[,colnames(data_test)!="y"])%*%coefficients(result_ridge,s="lambda.min")
  MSPE$ridge[t] <- mean((data_test$y-yhat_ridge)^2)
  cat("t=",t," ")
}

##ggplotで箱ひげ図を描く
g <- ggplot(gather(MSPE), aes(y=value,x=key))
g <- g + geom_boxplot()
g <- g + xlab("Method")
g <- g + ylab("MSPE")

plot(g)
```

ばらつきも含めて見たいので結果を[箱ひげ図](/python-boxplot)で表します。\
通常の線形回帰はglmと表記してあります。

・$p=5$の結果

![](/wp-content/uploads/2018/02/3e65cfe65c91ce8abec5ba00e6b79e5e.png)

説明変数の数も多くないのであまり手法の差は見えません。

・$p=10$の結果

![](/wp-content/uploads/2018/02/34cbcc674ea1b1a90be689239388ac0c.png)

通常の線形回帰とRidgeは同じくらいですが、Lassoが良い精度を示しています。\
説明変数が多くなってスパース性が上がったことが原因でしょう。

・$p=50$の結果

![](/wp-content/uploads/2018/02/19b6ed49d89b805d0ecf35d2764ce638.png)

通常の線形回帰が最も悪くなり、Lassoが最も良くなりました。

差もかなり大きいです。

:::point
**通常の線形回帰が悪い理由**：サンプルサイズが100なのに対して説明変数が50次元なので推定が不安定になった。\
**Ridgeが２番目の理由**：正則化しているので推定は安定しているがスパース性を考慮していないのでLassoより劣った。
:::

よって、<mark>説明変数が多くて回帰係数がスパースのとき、Lassoは予測精度が通常の線形回帰・Ridgeと比べて良い</mark>ことがわかりました！

こんな感じで手法の精度を比べたいときは、シミュレーションを繰り返して比較しましょう！

しかし、「**今回は人工データだからたくさんデータを生成できるけど実データで精度比較したいときはどうするの？**」という疑問が出てくると思います。その場合には<mark>交差検証法</mark>を使って比較します。

<https://toukei-lab.com/cross-validation>

Ridge・Lassoについて理論から知りたい方はこちらの記事をご覧ください！

<https://toukei-lab.com/%e3%82%b9%e3%83%91%e3%83%bc%e3%82%b9%e6%8e%a8%e5%ae%9a>

また、機械学習手法やデータサイエンスの広範な範囲について知りたい方はぜひ以下の記事をチェックしてみてくださいね！

<https://toukei-lab.com/machine-learning-study>

<https://toukei-lab.com/data-scientist>

さらに当メディアでは機械学習やデータサイエンスについて幅広く学べるスクール「スタアカ」を展開しています！

:::box
**[業界最安値のAIデータサイエンス特化スクール「スタアカ」](https://toukei-lab.com/achademy/)**
:::

ぜひチェックしてみてください！