mikutaifukuの雑記帳

個人的な雑記帳。データ分析とか読んだ本の感想とか。

Fitbitデータでスパース推定に入門してみる①　〜lassoで睡眠効率に影響を与える要因を探る〜

Fitbit R 機械学習 スパース推定

はじめに

2018年2月22日～2月23日(2日間)、統計数理研究所の2017年度リーディングDAT講座「【L-B2】 機械学習とデータサイエンスの現代的手法」（2017年度リーディングDAT）に参加してきました。（【L-B1】 実践ベイズモデリング」は抽選に外れたので来年あればリベンジしたい。） 従来の私ならこのような講座に参加するだけで満足しがちなのですが、しっかりと知識が定着するようにブログにまとめたいと思います。

講座内容は「カーネル法」「因果推論」「ニューラルネットワーク」「教師なし学習と行列分解」「スパース推定」と盛りだくさんでしたが、 今月末に阪大MMDS「機械学習・統計学 スプリングキャンプ」に参加するので、予習も兼ねてここでは「スパース推定（lassoを中心）」について取り上げます。

スパースとは

「スパース」とは「疎」という意味で、スカスカというイメージです。

私もこのような知的なツイートができるようになりたいです。

スパース推定とは

スパース推定とは統計学における推定方法の１つ（最小二乗法や最尤推定みたいな）で、特徴は

  • パラメータ推定
  • 変数選択（特徴量選択）

を同時に行うことにあります。データはものすごい大規模だけれど、意味ある情報はごく一部である時に使われます。具体例としては、以下のようなものがよく挙げられます。

  • 沢山ある遺伝子のうち、特定の病気に関係するものは少ない
  • 画像データでは、隣同士のピクセルの色は同じような色であることがほとんど

仕組みとしては、パラメータを0とすることで変数選択を可能とします。 f:id:mikutaifuku:20180304180916p:plain パラメータ推定値に0が含まれ、非0となる推定値がまばらにある様が、疎（スパース）のように見えるため、この推定方法をスパース推定とよびます。

一般的なパラメータ推定

f:id:mikutaifuku:20180302165636p:plain ただ、上記のような線形回帰モデルの場合、サンプルサイズより特徴量の方が大きい時（n < p）や、特徴量間の相関が高い時（多重共線性）は、最小二乗推定量は存在しないことが知られています。そこで、スパース推定の出番です。

スパース推定

スパース推定の考え方は非常にシンプルで上記に書いた最小二乗推定量を算出する式の後ろに、正則化項（罰則項）をつけて、パラメータが大きく推定されるのを防ぎます。特に、正則化項にL1ノルムを考えた場合をLassoと言います。Lassoはパラメータを0と推定することで、変数選択をすることが可能です。例えば相関が高い２つの説明変数があったとき、片方の説明変数のパラメータが0と推定されて、その変数はモデルから除外されます。「なぜLassoはパラメータを0に推定することができるのか？」はこの辺りの記事（https://research.preferred.jp/2014/05/group-lasso/）を参考にしていただければと思います。

f:id:mikutaifuku:20180302165849p:plain 医学領域では、DNAマイクロアレイデータ（遺伝子の発言量）やレセプトデータといった高次元小標本データは多く存在しており、スパース推定の応用が期待されています。

やってみた

タイトルにあるように、Fitbitデータに対してスパース推定を試してみました。fitbitデータ及び自作の生活行動データ（毎日Excelで手入力しています。）を利用して、睡眠効率に影響を与える要因を探ってみます。

データ

自作の生活行動データは欠損が非常に多く、データ数も多くありませんが気にせずやってみます。 項目は

  • totalTimeInBed：ベッドにいた時間（≠睡眠時間）
  • efficiency：睡眠効率（睡眠時間÷ベッドにいた時間）
  • Work ：仕事の日は1
  • Drink：お酒を飲んだ日は1
  • Glass：眼鏡をかけていた日は1、コンタクトの日は0
  • Bath：湯船につかった日は1
  • Stretch：お風呂上がりにストレッチをした日は1
  • weekday：平日は1
  • step8000：8000歩以上歩いた日は1

となります。

> glimpse(df)
Observations: 43
Variables: 9
$ totalTimeInBed <int> 423, 357, 590, 426, 463, 426, 440, 417, 381, 588...
$ efficiency     <dbl> 0.9290780, 0.9607843, 0.9067797, 0.9436620, 0.95...
$ Work           <int> 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, ...
$ Drink          <int> 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, ...
$ Glass          <int> 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, ...
$ Bath           <int> 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, ...
$ Stretch        <int> 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, ...
$ weekday        <int> 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, ...
$ step8000       <int> 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, ...

本来であれば、しっかりと可視化をしてデータの分布などを確認すべきですが割愛します。

lasso

glmnetパッケージを使用すれば、簡単にlassoを実施する事ができます。

y <- unlist(select(df, efficiency))
x <- scale(as.matrix(select(df, -efficiency)))

library("glmnet")
fit <- glmnet(x, y, standardize = FALSE)
plot(fit, label = TRUE)

f:id:mikutaifuku:20180304132526p:plain

簡単にsolution path（解パス）を求める事ができました。solution pathとは、正則化パラメータλを変化させた時の、回帰係数の変化を示したものになります。パラメータが大きい時（罰則が大きい時）は、一部の変数のパラメータが0に推定されている事がわかります。

ただ、各線がどの変数を示しているのかが分からず、グラフもイマイチです。「ggplot2形式でsolution pathを描けないのか」と思い調べた結果、やっぱりありました。しかも綺麗にまとめてくださっていて素晴らしいですね。これがRの良いところだと思います。

qiita.com

library("ggfortify")
autoplot(fit)

f:id:mikutaifuku:20180304151133p:plain

裏側でggplotが動いているのだと思います。個人的には「背景白色派」なので、もしかして…と思い下記を実行しました。

autoplot(fit) + theme_pander() + scale_color_pander()

f:id:mikutaifuku:20180304151424p:plain

いい感じになりました。この辺りは雰囲気でggplot2をやっている人間なので、オブジェクトの考え方含め一度しっかりと勉強しようと思います。

ハイパーパラメータの選択

上記の結果を見てわかるように回帰係数は、ハイパーパラメータの選び方によって大きく変化します。従って、最適なハイパーパラメータを交差検証法（cross validation）等によって選択する必要があります。これも簡単にRで実行する事ができます。

set.seed(0)
cvfit <- cv.glmnet(x, y, standardize = FALSE, nfolds = 10)
autoplot(cvfit) + theme_pander() + scale_color_pander()

f:id:mikutaifuku:20180304152045p:plain

左側の点線が交差検証法で算出した２乗誤差を最小にするλ (= λmin) を表しています。
λminの時の回帰係数は以下のように取得にします。 weekday（平日）とstep8000（歩数が8000歩以上）が変数として選択されました。

> coef(cvfit, s = "lambda.min")
9 x 1 sparse Matrix of class "dgCMatrix"
                          1
(Intercept)    0.9431179145
totalTimeInBed .           
Work           .           
Drink          .           
Glass          .           
Bath           .           
Stretch        .           
weekday        0.0025008443
step8000       0.0006026092

結果の解釈

推定結果やsolution pathから自分なりに結果を解釈してみます。

  • weekdayは睡眠効率にプラス
    • 平日はベッドの中でダラダラする事が少ないからだと思われる
  • step8000は睡眠効率にプラス
    • よく運動している（歩いている日）はぐっすり眠りにつく事ができるので睡眠効率にプラスの影響を与えそう
  • Drink（飲酒）は睡眠効率にプラス
    • 家に帰ってきた瞬間に着替えてベッドに入り爆睡するから睡眠効率が良いのかも…。これは意外な結果でした。

λが小さい時は各変数間に相関関係が見られるので（多重共線性）、少し回帰係数の推定結果が微妙な気がします。Workとweekdayはほぼ同じ日フラグが立っているのですが、おそらく相関が強いのでプラスとマイナスに両極端な推定結果が見られているのだと思います（たぶん）。ただ、cross validationで得られた結果ではweekdayしか選択されなかったのでとりあえずは大丈夫なのかなと思います（私見）。

さいごに

今回のエントリでは、スパース推定入門という事で、Fitbitデータの睡眠効率に影響を与える要因をlassoで探索しました。この記事を作成するにあたって昨年参加した(C)スパース推定(平成29年度)の資料を参考に書きました。Rのコードの部分は岩波データサイエンス Vol.5の前半部分を参考にしました。今後この辺りの本をしっかり読んで（スパース推定法による統計モデリング (統計学One Point 6)）理論面もしっかりと学ぼうと思います。

最後に個人的に気になった点、学びたい点をMEMOしておきます。

  • 説明変数が0か1のデータであった時にも、問題なくスパース推定出来るのか？説明変数の標準化について、ある行為を受けたか受けなったかみたいな01変数を正規化する事の意味・解釈の仕方がよくわからない。
  • ggplot2のオブジェクトについて諸々。雰囲気でRをやる人間なので、そろそろ真面目に勉強したい

分からないことだらけですね。まだまだINPUTが足らないようなので、引き続き勉強頑張ります。

なお、本エントリは雰囲気で書いた部分も多いので、間違いがあれば遠慮なくご指摘くださいませ。

mikutaifuku