Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2012-12-26

MikuHatsune2012-12-26

LASSOをやってみる

この1年の最後のほうで、LASSOというものがやたらと出てきたのでやってみる。

数式や意味などはこちらで解説されているので、実際にデータがあったらどうするかということでいつもどおりやる。

今回もいつもどおりぐぐると先生にあたるので、

からglmnetパッケージの例をそのままやってみる。

回帰

二値判別

多値判別

 

回帰

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
.
.
.

f:id:MikuHatsune:20121225161015p:image

この結果のグラフで、横軸がペナルティの強さ、縦軸が回帰係数の推定結果を示している、らしい。

それぞれの折れ線が説明変数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")

f:id:MikuHatsune:20121225163324p:image

 

交差検証で最適な¥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

 

二値判別

二値判別にはmdrrデータセットを使ってみる。

Svetnik et al (2003) describe these data: "Bakken and Jurs studied a set of compounds originally discussed by Klopman et al., who were interested in multidrug resistance reversal (MDRR) agents. The original response variable is a ratio measuring the ability of a compound to reverse a leukemia cell’s resistance to adriamycin. However, the problem was treated as a classification problem, and compounds with the ratio >4.2 were considered active, and those with the ratio <= 2.0 were considered inactive. Compounds with the ratio between these two cutoffs were called moderate and removed from the data for twoclass classification, leaving a set of 528 compounds (298 actives and 230 inactives). (Various other arrangements of these data were examined by Bakken and Jurs, but we will focus on this particular one.) We did not have access to the original descriptors, but we generated a set of 342 descriptors of three different types that should be similar to the original descriptors, using the DRAGON software."

The data and R code are in the Supplimental Data file for the article.

で、

mdrrDescr:the descriptors

mdrrClass:the categorical outcome ("Active" or "Inactive")

となっているので、xをmdrrDescr、yをmdrrClassとして使えばおk。

familyをbinomial指定して計算する。データ行列がちょっと大きいので計算に時間がかかる。

data(mdrr, package = "caret")
mdrr.glmnet <- glmnet(as.matrix(mdrrDescr), mdrrClass, family="binomial")
plot(mdrr.glmnet)

f:id:MikuHatsune:20121225161016p:image

回帰と同じように、最適な¥lambdaを求めて、判別分析をする。

mdrr.cv.glmnet <- cv.glmnet(as.matrix(mdrrDescr), mdrrClass, family="binomial")
table(as.numeric(mdrrClass), predict(mdrr.glmnet, as.matrix(mdrrDescr), s=mdrr.cv.glmnet$lambda.min, type="class"))
    Active Inactive
  1    268       30
  2     45      185

 

多値判別

最後に多値判別を行う。これにはおなじみのirisを使う。

irisは4つのパラメータに3種類の種別がついている。

iris.glmnet <- glmnet(as.matrix(iris[, -5]), iris[, 5], family="multinomial")
par(mfrow=c(1, 3))
plot(iris.glmnet)

f:id:MikuHatsune:20121225161017p:image

同じく予測する。

iris.cv.glmnet <- cv.glmnet(as.matrix(iris[, -5]), iris[, 5], family="multinomial")
table(iris[, 5], predict(iris.glmnet, as.matrix(iris[, -5]), s=iris.cv.glmnet$lambda.min, type="class"))
             setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         48         2
  virginica       0          1        49