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関数で最適なの値を求める。これを用いると、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")
交差検証で最適なが得られましたが、それは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