回帰

glmnetの引数は
x:数値からなる行列
y:応答を表すベクトル
family:二値判別分析、多値判別分析も可能。ポアソンとかcoxとかもある。
となっているので、これに従ってやる。
 
とりあえず、caretパッケージのcarsを使う。
このデータは、

data frame of the suggested retail price (column Price) and various characteristics of each car (columns Mileage, Cylinder, Doors, Cruise, Sound, Leather, Buick, Cadillac, Chevy, Pontiac, Saab, Saturn, convertible, coupe, hatchback, sedan and wagon)

t(head(cars))
                   1        2        3        4        5        6
Price       22661.05 21725.01 29142.71 30731.94 33358.77 30315.17
Mileage     20105.00 13457.00 31655.00 22479.00 17590.00 23635.00
Cylinder        6.00     6.00     4.00     4.00     4.00     4.00
Doors           4.00     2.00     2.00     2.00     2.00     2.00
Cruise          1.00     1.00     1.00     1.00     1.00     1.00
Sound           0.00     1.00     1.00     0.00     1.00     0.00
Leather         0.00     0.00     1.00     0.00     1.00     0.00
Buick           1.00     0.00     0.00     0.00     0.00     0.00
Cadillac        0.00     0.00     0.00     0.00     0.00     0.00
Chevy           0.00     1.00     0.00     0.00     0.00     0.00
Pontiac         0.00     0.00     0.00     0.00     0.00     0.00
Saab            0.00     0.00     1.00     1.00     1.00     1.00
Saturn          0.00     0.00     0.00     0.00     0.00     0.00
convertible     0.00     0.00     1.00     1.00     1.00     1.00
coupe           0.00     1.00     0.00     0.00     0.00     0.00
hatchback       0.00     0.00     0.00     0.00     0.00     0.00
sedan           1.00     0.00     0.00     0.00     0.00     0.00
wagon           0.00     0.00     0.00     0.00     0.00     0.00

となっていて、価格(Price)と走行距離(Mileage)は連続値、エンジン?(Cylinder)とDoorは自然数、それ以外は0or1っぽい。

cars.glmnet <- glmnet(as.matrix(cars[, -1]), cars[, 1])
plot(cars.glmnet)
Call:  glmnet(x = as.matrix(cars[, -1]), y = cars[, 1]) 

      Df   %Dev  Lambda
 [1,]  0 0.0000 6513.00
 [2,]  1 0.0738 5934.00
 [3,]  1 0.1351 5407.00
.
.
.


この結果のグラフで、横軸がペナルティの強さ、縦軸が回帰係数の推定結果を示している、らしい。
それぞれの折れ線が説明変数xに相当する、らしい。
ペナルティは左が強く、右に行くに従って弱くなっている、らしい。
 
次に、cv.glmnet関数で最適な\lambdaの値を求める。これを用いると、predict.glmnet関数で未知データの回帰が可能になる。

cars.cv.glmnet <- cv.glmnet(as.matrix(cars[, -1]), cars[, 1])
prd.cars.glmnet <- predict(cars.glmnet, as.matrix(cars[, -1]), s=cars.cv.glmnet$lambda.min)
plot(cars[, 1], prd.cars.glmnet, xlim=range(cars[, 1], prd.cars.glmnet), ylim=range(cars[, 1], prd.cars.glmnet), xlab="cars price", ylab="predict")


 

交差検証で最適な\lambdaが得られましたが、それは65ステップ目であることがわかります。そこで、65ステップ目における各係数を出力させてみたところ以下のようになりました。これをみると、やはり車種がキェでラック(Cadillac)の場合は価格に反映されますが、革張り(Leather)の効果は低く、クーペ(Coupe)かセダン(Sedan)かどうかは価格には影響しないことがわかります。

 match(cars.cv.glmnet$lambda.min, cars.glmnet$lambda)
[1] 65
cars.glmnet$beta[, match(cars.cv.glmnet$lambda.min, cars.glmnet$lambda)]
                     [,1]
Mileage        -0.1818766
Cylinder     3636.7109817
Doors        -631.4904976
Cruise        340.5460071
Sound         387.1736172
Leather       756.9417865
Buick         928.2740498
Cadillac    13399.6038147
Chevy        -487.6230684
Pontiac     -1302.5800365
Saab        12275.1021808
Saturn          0.0000000
convertible 11016.5702170
coupe           0.0000000
hatchback   -1884.5001090
sedan           0.0000000
wagon        4328.7502752

Cadillacは高いけど、coupeやsedanは0なのでそういうことなんだろう。
 
ここで得られた係数とそれぞれの値の内積を算出すると、回帰式のY切片が得られるらしい。

summary(as.matrix(cars[, -1]) %*% cars.glmnet$beta[, match(cars.cv.glmnet$lambda.min, cars.glmnet$lambda)])
       V1       
 Min.   : 4573  
 1st Qu.:10644  
 Median :16141  
 Mean   :18011  
 3rd Qu.:21721  
 Max.   :53626