演習の概要

注意事項

  • 慶應義塾大学SFCで開講している「統計解析」の授業履修者向けの演習用ページです。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.6.0を使っています。
  • Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

データ分析の準備

パッケージのインストール

  • 今回の演習では、以下のパッケージを使います。
install.packages("ggplot2")
install.packages("glmnet")
install.packages("lars")
library(ggplot2)
library(glmnet)
library(lars)

モデル推定

データの読み込み

  • 今回の演習では、larsパッケージに含まれるdiabetesデータを用います。
data("diabetes")
head(diabetes)

線形回帰モデル

  • まず重回帰モデルを推定します。
dia.lm <- lm(y~x, data=diabetes)
summary(dia.lm)

Call:
lm(formula = y ~ x, data = diabetes)

Residuals:
     Min       1Q   Median       3Q      Max 
-155.829  -38.534   -0.227   37.806  151.355 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  152.133      2.576  59.061  < 2e-16 ***
xage         -10.012     59.749  -0.168 0.867000    
xsex        -239.819     61.222  -3.917 0.000104 ***
xbmi         519.840     66.534   7.813 4.30e-14 ***
xmap         324.390     65.422   4.958 1.02e-06 ***
xtc         -792.184    416.684  -1.901 0.057947 .  
xldl         476.746    339.035   1.406 0.160389    
xhdl         101.045    212.533   0.475 0.634721    
xtch         177.064    161.476   1.097 0.273456    
xltg         751.279    171.902   4.370 1.56e-05 ***
xglu          67.625     65.984   1.025 0.305998    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 54.15 on 431 degrees of freedom
Multiple R-squared:  0.5177,    Adjusted R-squared:  0.5066 
F-statistic: 46.27 on 10 and 431 DF,  p-value: < 2.2e-16

Ridge回帰モデル

  • Ridge回帰モデルとLasso回帰モデルは、glmnetパッケージのglmnet()関数を用いて推定します。
  • Ridge回帰モデルは、glmnet()関数の引数alpha=0を指定して推定します。
dia.ridge <- glmnet(diabetes$x, diabetes$y, 
                    # family="gaussian", 
                    alpha=0)
  • 回帰係数と\(log(\lambda)\)との関係を可視化します。
plot(dia.ridge, label=TRUE)

  • 回帰係数と\(\lambda\)との関係を可視化します。
plot(dia.ridge, xvar="lambda", label=TRUE)

  • クロスバリデーションによって、最適な\(\lambda\)を決定します。
dia.ridge.cv <- cv.glmnet(diabetes$x, diabetes$y, 
                    # family="gaussian", 
                    alpha=0)
  • クロスバリデーションの結果を可視化します。
plot(dia.ridge.cv)