臨床予測モデルの開発に必要なサンプルサイズの算出(Calculating the sample size required for developing a clinical prediction model)
概要
ヘルスケア領域のデータ分析に携わっていると、よく質問されるのが
「どれくらいサンプルサイズが必要なの?」
という質問です。統計的仮説検定であれば、
「有意水準、検出力…等を指定して、必要サンプル数を求めることが可能だと思います」
とそれらしく回答しているのですが、予測モデルを構築する場合は、
「なるべく多くのサンプル集める。特徴量の次元数の10倍のデータが必要みたいなことは、経験則として提示されていますね…」
というような感じで、もごもごと回答しています。
そこで、google検索してみたところ「Calculating the sample size required for developing a clinical prediction model」という論文が存在していたようなので、流し読みしてみました。
この論文では、binary outcome(二値の目的変数)に対しての、logistic regression(ロジスティック回帰モデル)を用いた場合に限り、確率的な性質を利用してのサンプルサイズの見積もり手法が提案されていました。しかも、Rのpmsampsizeパッケージも紹介されていました。
そこで、今回の記事では、背景や細かい理論は気にせず、個人的なメモとしてRのpmsampsizeパッケージの使い方を書いています。
pmsampsizeの使い方
まずは必要なlibraryを読み込みます。最後にちょっとしたデータ加工と可視化をするので、tidyverseも読みます。
library(pmsampsize) library(tidyverse)
pmsampsizeのpmsampsize関数を用いることで、サンプルサイズを見積もることが出来ます。
例えば、目的変数が0/1の2値変数、ロジスティック回帰モデルで、説明変数が5個、1の出現割合が0.3、AUC=0.70とすると、以下のような結果が出力され、サンプルサイズ412と見積もることができました。
> pmsampsize(type = "b", cstatistic = 0.70, parameters = 5, prevalence = 0.3) Given input C-statistic = 0.7 & prevalence = 0.3 Cox-Snell R-sq = 0.1029 NB: Assuming 0.05 acceptable difference in apparent & adjusted R-squared NB: Assuming 0.05 margin of error in estimation of intercept NB: Events per Predictor Parameter (EPP) assumes prevalence = 0.3 Samp_size Shrinkage Parameter Rsq Max_Rsq EPP Criteria 1 412 0.900 5 0.1029 0.71 24.72 Criteria 2 131 0.743 5 0.1029 0.71 7.86 Criteria 3 323 0.900 5 0.1029 0.71 19.38 Final 412 0.900 5 0.1029 0.71 24.72 Minimum sample size required for new model development based on user inputs = 412, with 124 events (assuming an outcome prevalence = 0.3) and an EPP = 24.72
また、1の出現割合を、0.1, 0.2, ..., 0.4と動かしたり、説明変数の個数を1,2,...,20と動かしたりすることで、様々なパターンのサンプルサイズを見積もることが出来ます。
df_pm <- tibble() p_s <- c(0.1, 0.2, 0.3, 0.4) for (p in p_s) { for (para in seq(1, 20, 1)) { pm_tmp <- pmsampsize(type = "b", cstatistic = 0.70, parameters = para, prevalence = p) df_tmp <- tibble(1の出現割合=p, パラメータ数=para, サンプル数=pm_tmp$sample_size) df_pm <- bind_rows(df_pm, df_tmp) } } df_pm$1の出現割合 <- factor(as.character(df_pm$1の出現割合)) g <- ggplot(df_pm ,aes(パラメータ数, サンプル数, colour=1の出現割合, group=1の出現割合)) + geom_line(size=1) + geom_point() + ggtitle("必要なサンプル") plot(g)
さいごに
- clinical prediction model(臨床予測モデル)、しかも0/1の2値のロジスティック回帰モデルに限定しての話で、おそらく一般的な機械学習だとこの限りではないと思います。
- また、論文内では、train/testのようなデータの分割(ホールドアウト検証)は推奨されておらず、ブートストラップ法による検証が推奨されており、このあたりの感覚は少し違和感があります。
- 雰囲気でブログを書いただけなので、おかしな事をやっていたり、間違っていたらこっそりご連絡ください。
参考
Calculating the sample size required for developing a clinical prediction model | The BMJ
MENTA登録してみました
知り合いが登録しているとのことで、チャレンジしてみることにしました。
【初心者〜中級者向け】R/Pythonでのデータ分析スキル習得に向けてサポートします! mikutaifuku|メンターに教えてもらおう
こういうのをきっかけに、ちょっとした弱い繋がりができればと期待しつつ、がんばってみます。
ggplot2でfacetごとのヒストグラムに色んな情報を付与する② 〜Tidy evalによる関数化と繰り返し処理〜
概要
前回の続きで、自分用メモです。
タイトルの通り、Tidy evalにより関数を作り*1、全ての変数に対して処理を繰り返し行います。 要は、以下のようなグラフを Sepal.Width, Petal.Length, Petal.Width 全てに対して作ります。
なお、前提知識として、ブログ下部の「参考」にあるdplyrの基礎知識(使い方)は必要です。「やってみた」系なので、特に有益な情報もないですが、予めご容赦ください。
関数化
データ加工処理
group_var, x_varは、quo()されたもので、関数の中で!!を用いてアンクオートします。他は前回のブログと何も変わりません。
Make_Data <- function(d, group_var, x_var){ df <- d %>% select(x = !!x_var, group = !!group_var) df_summary <- df %>% group_by(group) %>% summarise(n = n(), mean_x = mean(x), min_x = min(x), max_x = max(x), sd_x =sd(x)) %>% mutate(group_label = paste0(group, " (n = ", n, ")"), mean_label = paste0("Mean : ", round(mean_x, 2)), sd_label = paste0("SD : ", round(sd_x, 2)), min_label = paste0("Min : ", round(min_x, 2)), max_label = paste0("Max : ", round(max_x, 2))) %>% mutate(text_label = paste0(mean_label,"\n",sd_label,"\n",min_label,"\n",max_label)) df_plot <- df %>% left_join(df_summary %>% select(group, group_label), by = "group") return(list(df_plot, df_summary)) }
可視化処理
x_varはquo()されたもので、関数の中で列名だけを使用したいので、quo_name(x_var)により文字列に変換します。他は前回のブログと何も変わりません。
Make_Plot <- function(d_summary, d_plot, x_var){ x_pos <- max(d_summary$max_x) y_pos <- Inf x_var_name <- quo_name(x_var) g <- ggplot(d_plot, aes(x)) + theme_minimal() + geom_histogram(fill="darkgrey") + geom_vline(data=d_summary, aes(xintercept=mean_x), linetype="dashed", color="darkred") + geom_vline(data=d_summary, aes(xintercept=min_x), linetype="dashed", color="darkred") + geom_vline(data=d_summary, aes(xintercept=max_x), linetype="dashed", color="darkred") + geom_text(data=d_summary, aes(x=x_pos, y=y_pos, label=text_label), colour="darkred", inherit.aes=FALSE, vjust="inward", hjust="inward") + facet_wrap(~group_label, ncol=1) + labs(x=x_var_name) }
繰り返し処理
以下のようにループ処理すると、冒頭に述べたようなグラフが全ての変数に対して作られます。*2
x_vars <- quos(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) for (i in 1:length(x_vars)) { data_list <- Make_Data(d=iris, x_var=x_vars[[i]], group_var=quo(Species)) df_plot <- data_list[[1]] df_summary <- data_list[[2]] g_x <- Make_Plot(d_summary=df_summary, d_plot=df_plot, x_var=x_vars[[i]]) ggsave(plot = g_x, filename = paste0(quo_name(x_vars[[i]]), ".png")) }
参考
さいごに
上記のような処理をする際に、Tidy evalなるものは役立つのではないかと勝手に思っています。 雰囲気で作成したので、間違っていたら、遠慮なくご指摘ください。 また、他に良いやり方があれば、教えてください。
ggplot2でfacetごとのヒストグラムに色んな情報を付与する① 〜基本作図編〜
概要
大した内容ではないですが、仕事で使ったので自分用メモです。なので、やや適当です。
以下のようなグラフを作ります。
具体的には、以下のような特徴を備えたグラフです。
- facetごとのヒストグラムに、要約統計量(平均値、最小値、最大値)の線を引く
- 各グラフの右上に、テキストで要約統計量を出力する
- グループ名の後ろに(n = xx)というような、n数を付与する
データ
データはirisデータを使用します。
グループごとの要約統計量
SpeciesごとにSepal.Lengthの要約統計用を計算します。また、出力用のテキストを作成します。
library(tidyverse) df <- iris %>% select(x = Sepal.Length, group = Species) df_summary <- df %>% group_by(group) %>% summarise(n = n(), mean_x = mean(x), min_x = min(x), max_x = max(x), sd_x =sd(x)) %>% mutate(group_label = paste0(group, " (n = ", n, ")"), mean_label = paste0("Mean : ", round(mean_x, 2)), sd_label = paste0("SD : ", round(sd_x, 2)), min_label = paste0("Min : ", round(min_x, 2)), max_label = paste0("Max : ", round(max_x, 2))) %>% mutate(text_label = paste0(mean_label,"\n",sd_label,"\n",min_label,"\n",max_label))
作成したデータは以下のような感じです。
> df_summary # A tibble: 3 x 12 group n mean_x min_x max_x sd_x group_label mean_label sd_label min_label max_label text_label <fct> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> 1 setosa 50 5.01 4.3 5.8 0.352 setosa (n = 50) Mean : 5.01 SD : 0.35 Min : 4.3 Max : 5.8 "Mean : 5.01\nSD : 0.35\nMin : 4.3\nMax : 5.8" 2 versicolor 50 5.94 4.9 7 0.516 versicolor (n = 50) Mean : 5.94 SD : 0.52 Min : 4.9 Max : 7 "Mean : 5.94\nSD : 0.52\nMin : 4.9\nMax : 7" 3 virginica 50 6.59 4.9 7.9 0.636 virginica (n = 50) Mean : 6.59 SD : 0.64 Min : 4.9 Max : 7.9 "Mean : 6.59\nSD : 0.64\nMin : 4.9\nMax : 7.9"
元データの加工
元データに各グループのn数を付与します。
df_plot <- df %>% left_join(df_summary %>% select(group, group_label), by = "group")
作成したデータは以下のような感じです。
> head(df_plot) x group group_label 1 5.1 setosa setosa (n = 50) 2 4.9 setosa setosa (n = 50) 3 4.7 setosa setosa (n = 50) 4 4.6 setosa setosa (n = 50) 5 5.0 setosa setosa (n = 50) 6 5.4 setosa setosa (n = 50)
プロット(完成)
いい感じのグラフできました。
x_pos <- max(df_summary$max_x) y_pos <- Inf g <- ggplot(df_plot, aes(x)) + theme_minimal() + geom_histogram(fill="darkgrey") + geom_vline(data=df_summary, aes(xintercept=mean_x), linetype="dashed", color="darkred") + geom_vline(data=df_summary, aes(xintercept=min_x), linetype="dashed", color="darkred") + geom_vline(data=df_summary, aes(xintercept=max_x), linetype="dashed", color="darkred") + geom_text(data=df_summary, aes(x=x_pos, y=y_pos, label=text_label), colour="darkred", inherit.aes=FALSE, vjust="inward", hjust="inward") + facet_wrap(~group_label, ncol=1) + labs(x="Sepal.Length") plot(g)
参考
いつも有益な情報を提供してくださる、株式会社ホクソエムのブログを参考にしました。 blog.hoxo-m.com
テキストがグラフ内におさまる方法を参考にしました。
ggplot(df, aes(x, y)) + geom_text(aes(label = text), vjust = "inward", hjust = "inward")
さいごに
パート②では、加工処理とプロット処理を関数化し、プロットを全ての変数に対して繰り返し行います。
雰囲気で作成したので、間違っていたら、遠慮なくご指摘ください。
Fitbitデータでスパース推定に入門してみる④ 〜Graphical lassoで変数間の関係性について調べる〜
はじめに
これまで述べてきたのは、主にモデルの回帰係数の推定に基づく話でした。
今回のエントリでは、これまでの回帰係数の推定ではなく、共分散構造に注目した手法を取り上げます。多数の変数間の依存関係を推定するために、ガウシアングラフィカルモデルにL1正則化の考え方を応用する、Graphical lassoとよばれる手法を使います。 なお、今回の記事も今までと同様に以下の書籍(5.1 グラフィカルモデルにおけるスパース推定)を参考に書いたものですので、詳細は書籍を参考にしていただければと思います。
スパース推定法による統計モデリング (統計学One Point)
- 作者: 川野秀一,松井秀俊,廣瀬慧
- 出版社/メーカー: 共立出版
- 発売日: 2018/03/08
- メディア: 単行本
- この商品を含むブログを見る
ガウシアングラフィカルモデル
が多変量正規分布 に従うとき、変数 との関係性を推定し、それをグラフ化する手法です。分散共分散行列の逆行列 とし、の成分をとします。詳細は書籍に譲りますが、の時に、とは他の変数を与えた元で条件付き独立と言います。
データ分析でよくみられる相関係数は、他の変数からの影響も含めて測っているのに対して、行列 は他の変数からの影響を除いて、(すなわちと以外の変数の条件のもとで)関連の強さを測るものになります。
例えば、以下のような行列 (左図)が得られている場合は、右のようなグラフで条件付き独立性の有無を表現することができます*1。
このグラフが表現するのは「条件付き独立性」があるか否かということだけであって、「変数間の関係の向き」を推定することができません。 もし、グラフの辺の向きを推定したい場合は、以前紹介したようにLingamモデルを仮定することで推定することができます。
パラメータ推定
グラフの辺の数は最大で頂点の2乗のオーダー程度になりますが、実際に辺が存在する数は少ないという意味ではスパースと言えます。データからスパースな分散共分散行列の逆行列(精度行列) を推定する手法をgraphical lassoと言います。 データ が観測され、これらが正規分布に従っていると仮定します。この時、標本から計算される分散共分散行列をSとすると、対数尤度関数に基づく損失関数は、
となります(定数省略)*2。これに正則化項を加えた式を最小化にすることで、精度行列 を推定します。
Rでやってみた
データ
Fitbitデータから得られる
- efficiency:睡眠効率
- timeInBed:ベッドにいた時間
- sedentaryMinutes:座っている時間
- restingHeartRate:安静時心拍数
- steps:歩数
- caloriesOut:消費カロリー
の項目を使用しました。
> glimpse(d) Observations: 210 Variables: 6 $ efficiency <int> 94, 95, 96, 97, 98, 94, 96, 98, 97, 99, 99, 95, 98, 97, 97, 97, 97, ... $ timeInBed <int> 413, 444, 568, 531, 308, 718, 332, 329, 303, 309, 309, 283, 295, 359... $ sedentaryMinutes <int> 884, 660, 669, 717, 658, 797, 653, 905, 941, 925, 943, 530, 1033, 94... $ restingHeartRate <int> 67, 64, 64, 65, 65, 65, 66, 67, 68, 67, 68, 67, 67, 66, 67, 68, 69, ... $ steps <int> 6488, 17042, 18910, 7970, 9075, 9175, 2161, 10136, 8588, 10373, 8997... $ caloriesOut <int> 2585, 3583, 3573, 2737, 3060, 3409, 2188, 2986, 2790, 2950, 2842, 25...
散布図行列を作成し、データの傾向をみてみます*3。コードはこのブログの最後に載せています。
- 歩数と消費カロリー(5行6列)は相関係数が0.92と強い正の相関があります。
- 睡眠効率とベッドにいた時間(1行2列)は相関係数が-0.16とやや負の相関があります。ベッドにだらだらと長い間居ると、睡眠効率が悪くなりそうなことから容易に想像がつきます。
- ベッドにいた時間と座って居る時間が(2行3列)は相関係数が-0.11とやや負の相関があります。普通に考えてベッドに長くいれば(≠座っている状態)、座る時間も短くなりそうです。
- 睡眠効率と消費カロリー(1行6列)は相関係数が0.18とやや正の相関があります。消費カロリーが大きい日(よく歩いたり、運動したりした日)は、効率の良い睡眠が取れてそうです。
というように、この結果からも色々なことが考えられます。(もっと色々相関があってほしかった…)
graphical lasso
glassoパッケージを使用します。プログラムはものすごく単純です*4。
library(glasso) library(igraph) # グラフを描くのに必要 cordat <- cor(d) fit.glasso <- glasso(cordat, rho = 0.4) fit.glasso$wi # 共分散逆行列 makeadjmat <- function(mat, cordat) { ans <- mat diag(ans) <- 0 ans <- ans != 0 ans <- matrix(as.numeric(ans), nrow(ans), nrow(ans)) colnames(ans) <- rownames(ans) <- colnames(cordat) ans } adjmat <- makeadjmat(fit.glasso$wi, cordat) # 隣接行列 g <- graph.adjacency(adjmat, mode = "undirected") plot(g,vertex.size = 70, vertex.shape = "rectangle", vertex.color = "white")
λ=0.1の時のグラフィカルモデルです
λ=0.4の時のグラフィカルモデルです。
- 安静時心拍数(RestingHeartRate)はどの変数とも関係がなさそうです。
- 歩数(steps)と消費カロリー(caloriesOut)は関係がありそうです。
推定された から、偏相関係数という量を算出して、可視化することもできますが、次回の機会にとっておきます。
相関行列の可視化Rコード
library("ggplot2") library("ggthemes") library("scales") library("GGally") N_col <- ncol(d) ggp <- ggpairs(d, upper='blank', diag='blank', lower='blank') for(i in 1:N_col) { x <- d[,i] colnames(x) <- "tmp_name" p <- ggplot(x, aes(tmp_name)) p <- p + theme(text=element_text(size=14), axis.text.x=element_text(angle=40, vjust=1, hjust=1)) bw <- (max(x)-min(x))/10 p <- p + geom_histogram(binwidth=bw, color='darkgrey', fill='darkgrey', alpha=0.5) p <- p + geom_line(eval(bquote(aes(y=..count..*.(bw)))), stat='density', alpha=0.5) p <- p + geom_label(data=data.frame(x=-Inf, y=Inf, label=colnames(d)[i]), aes(x=x, y=y, label=label), hjust=0, vjust=1) p <- p + theme_bw() + scale_color_pander() ggp <- putPlot(ggp, p, i, i) } zcolat <- seq(-1, 1, length=81) zcolre <- c(zcolat[1:40]+1, rev(zcolat[41:81])) for(i in 1:(N_col-1)) { for(j in (i+1):N_col) { x <- as.numeric(unlist(d[,i])) y <- as.numeric(unlist(d[,j])) r <- cor(x, y, method='spearman', use='pairwise.complete.obs') zcol <- lattice::level.colors(r, at=zcolat, col.regions=colorRampPalette(c(scales::muted('red'), 'white', scales::muted('blue')), space='rgb')(81)) textcol <- ifelse(abs(r) < 0.4, 'grey20', 'white') ell <- ellipse::ellipse(r, level=0.95, type='l', npoints=50, scale=c(.2, .2), centre=c(.5, .5)) p <- ggplot(data.frame(ell), aes(x=x, y=y)) p <- p + theme_bw() + theme( plot.background=element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.border=element_blank(), axis.ticks=element_blank() ) p <- p + geom_polygon(fill=zcol, color=zcol) p <- p + geom_text(data=NULL, x=.5, y=.5, label=100*round(r, 2), size=6, col=textcol) ggp <- putPlot(ggp, p, i, j) } } for(j in 1:(N_col-1)) { for(i in (j+1):N_col) { x <- d[,j] y <- d[,i] colnames(x) <- "tmp_x_name" colnames(y) <- "tmp_y_name" p <- ggplot(data.frame(x, y), aes(x=tmp_x_name, y=tmp_y_name)) p <- p + theme(text=element_text(size=14), axis.text.x=element_text(angle=40, vjust=1, hjust=1)) p <- p + geom_point(size=1, color='darkgrey', fill='darkgrey', alpha=0.5) p <- p + theme_bw() + scale_color_pander() p <- p + theme(text=element_text(size=14), axis.text.x=element_text(angle=40, vjust=1, hjust=1)) ggp <- putPlot(ggp, p, i, j) } } print(ggp)
さいごに
今回のエントリでは、graphical lassoという手法を用いてFitbitデータの変数間の関係性をみました。またgraphical lassoによる異常検知の手法というのも存在しているらしく、相当変なデータを使用してもおかしな結果を出しにくいという意味で実用性が高いらしいので*5、機会があればチャレンジしてみようかと思います。
ということで、4回に渡ってスパース推定入門エントリを書いてきましたが、このシリーズは一旦これにて終わりにしようと思います。 今回はスパース推定法による統計モデリング (統計学One Point)の本を読んで自分なりに整理したわけですが、計算アルゴリズムを飛ばしたり、式を手計算できっちりと追いかけた訳ではないので、時間を見つけて再勉強しようと思います。
短時間で書いたので間違いがあるかもしれません。その際は優しくご指摘いただきますと幸いです。
mikutaifuku
*1:書籍p.107を参考に書きました
*2:詳細な計算手順等は書籍を参考してください
*3:散布図行列を描くには (corrplot, pairs, GGally) - StatModeling Memorandumの記事を参考にさせて頂きました。本当にすごいですね。
*4:書籍に書かれていたものを写しました
*5:岩波データサイエンス Vol.5に書いていました
Fitbitデータでスパース推定に入門してみる③ 〜心拍数データに対してfused lassoを使い、変化点を検出する〜
はじめに
前回の記事では、個々の変数を選択するための様々なスパース推定の手法をまとめ、Fitbitデータに適用することでそれぞれの手法を比較しました。
今回のエントリでは、少し趣向を変えて心拍数データに対して、fused lassoと呼ばれる隣合う変数の関係性を加味した手法使ってみます。なお、今回のエントリは、書籍の第3章のfused lassoの部分を自分なりに整理してまとめただけなので、丁寧に学びたい方は、書籍を読んでいただければと思います。
スパース推定法による統計モデリング (統計学One Point)
- 作者: 川野秀一,松井秀俊,廣瀬慧
- 出版社/メーカー: 共立出版
- 発売日: 2018/03/08
- メディア: 単行本
- この商品を含むブログを見る
fused lassoの概要
fused lassoは以下のような性質をもつデータに対して提案された手法です。
- p個の説明変数に順序関係がある
- 隣接する説明変数のいくつかは目的変数に同程度の寄与がある
最小二乗推定量を算出する式の後ろに、正則化項にL1ノルムと隣接する係数の差の絶対値を考えたものです。
この式を最小にすることでパラメータβを推定します。正則化項として係数の差の絶対値を考えているため、推定値が等しくなるように働くという仕組みです。
時系列データへの応用
Fitbitの睡眠時心拍数データ*1を使って、時間的に変化する値の推定を考えます。隣り合う時点の値はほぼ等しく、ときどきある時点でジャンプしています。
時刻tでの観測値をとし、
とモデリングします。εは観測ノイズです。隣り合う時点でβの値が変わりにくいという性質を、正則化項
で表現します。この項の効果で不連続に飛ぶ点の数を少なくしながら推定を行います。つまり、以下のような式を最小にする事で推定パラメータを求めます。
これは、上記の(1)式においてn = p, Xが単位行列、λ1=0を考えたのと同じです。*2
Rでやってみた
genlassoパッケージのfusedlasso1dを用いることでfused lassoによる推定を行えます。
データ
> glimpse(df) Observations: 414 Variables: 2 $ time <dttm> 2018-01-11 00:27:00, 2018-01-11 00:28:00, 2018-01-11 0... $ value <int> 102, 110, 106, 84, 79, 80, 79, 78, 77, 80, 77, 83, 83, ...
fused lasso
library(genlasso) y <- unlist(df$value) res <- fusedlasso1d(y) # fused lasso plot(res, col="grey") # plot
正則化パラメータλを変えた時のフィッティングを折れ線で示しています。L1正則化により変化点が検出できそうです。が、デフォルトのplotがわかりにくいので自作します。
resには様々なλ(全部で383個)において推定されたパラメータβが格納されています。適当なλを選択し(今回は50, 75, 100, 125, 150番目を使いました)、可視化しやすいようにデータを縦持ち(tidy化)にします。
library("dplyr") library("tidyr") step <- c(50,75,100,125,150) beta_df <- data.frame(res$beta)[, step] # 指定した列のβを抽出 beta_df$time <- df$time # 時点の列を追加 beta_df$HeartRate_ob <- df$value # 実測値の列を追加 colnames(beta_df) <- c(res$lambda[step], "time", "HeartRate_ob") # ヘッダー名をλの値に変更 # tidy化 plot_df <- beta_df %>% gather(key="lambda", value = "HeartRate_es", -time, -HeartRate_ob) %>% mutate(lambda = paste0("λ = ", round(as.double(lambda), 2)))
データが良い感じになりました。
> glimpse(plot_df) Observations: 2,070 Variables: 4 $ time <dttm> 2018-01-11 00:27:00, 2018-01-11 00:28:00, 2018-... $ HeartRate_ob <int> 102, 110, 106, 84, 79, 80, 79, 78, 77, 80, 77, 8... $ lambda <chr> "λ = 13", "λ = 13", "λ = 13", "λ = 13", "λ = 13"... $ HeartRate_es <dbl> 101.66667, 101.66667, 101.66667, 84.00000, 82.11...
可視化
library("ggplot2") library("ggthemes") ggplot(plot_df, aes(time, HeartRate_ob)) + geom_point(color = "darkgrey", alpha = 0.5) + geom_line(aes(time, HeartRate_es, color=lambda)) + facet_wrap(~lambda, scales = "free", ncol = 1) + theme_pander() + scale_color_pander() + theme(legend.position = 'none') + xlab("") + ylab("Heart Rate")
そこそこ綺麗なグラフが完成しました。λが大きいと直線に近づき、λが小さいと観測値に沿った折れ線が得られていることがわかります。λ=4のグラフを見ると4:00や6:00頃の変化点をしっかりと検出できています。
さいごに
今回のエントリでは、Fitbitの心拍数データに対してfused lassoを適用し、変化点の検出を試みました。ただ、実務においては、このように全てのデータが揃ってから変化点を検出しても遅く、要望として例えば「異常を5分前に予兆したい」みたいなことが挙げられます。が、この手法では残念ながら解決できません(おそらく)。という事で、これら課題をどうやったら解決できそうか、うまく問題設定を考えてチャレンジしてみようと思います。
このブログではパッケージを駆使して「やってみた」系の内容がほとんどです(てか全部)。自分でしっかり頭を使ってモデリングをしている訳でもなければ、自分で数式をしっかり理解した上でアルゴリズムを実装している訳でもありません。この状況は個人的にはなかなかヤバイと思いながらも、時間コストを考えるとこれが限界な気がしなくもありません。色々難しいですね、ちょっと色々考えてみます。
引き続き頑張ります。
mikutaifuku
Fitbitデータでスパース推定に入門してみる② 〜lasso, elastic net, SCAD, adaptive lassoの概要とRでの数値例〜
はじめに
前回の記事では、スパース推定の概要と代表的な手法であるlassoについてまとめ、fitbitデータに適用することで、睡眠効率に影響を与える要因を探索しました。
今回のエントリでは、lasso以外の様々なスパース正則化法を紹介し、前回と同様にFitbitデータに適用することで、それぞれの手法を比較したいと思います。なお、今回のエントリは、以下の書籍の第1・2章を自分なりに簡単にまとめた*1だけなので、丁寧に学びたい方は、書籍を読んでいただければと思います。
スパース推定法による統計モデリング (統計学One Point)
- 作者: 川野秀一,松井秀俊,廣瀬慧
- 出版社/メーカー: 共立出版
- 発売日: 2018/03/08
- メディア: 単行本
- この商品を含むブログを見る
lassoの概要と弱点
概要
lassoは以下のように、最小二乗推定量を算出する式の後ろに、正則化項にL1ノルムを考えたものです。 Lassoはパラメータを0と推定することで、変数選択をすることが可能です。例えば相関が高い2つの説明変数があったとき、片方の説明変数のパラメータが0と推定されて、その変数はモデルから除外されます。
弱点
しかし、lassoは2つの問題があることが知られています。
- 相関の高い2つの説明変数が目的変数に関係している場合、一方の説明変数しかモデルに含まれない
- サンプルサイズより説明変数の数の方が大きい場合(所謂 n < p の場合)、高々サンプルサイズ分の変数までしか選択されない
これらの問題を解決するという観点でelastic netが提案されました。さらにオラクル性質(oracle property)*2を持たせるという観点で、正則化項として非凸正則化項を使うSCADやadaptive lassoが提案されています。以降ではこれらの概要を説明してから、Rで数値計算してみます。
lasso正則化項の拡張
elastic net
elastic netは上記2つの問題を解決させるために提案された手法です。最小2乗推定量を算出する式の後ろに、下記のような正則化項を加えたものを最小化することで、パラメータを推定します。
これはL1ノルムとL2ノルムを合わせた形になっており、α=1の時はlasso, α=0の時はリッジとなります。このような正則化項を取ることで、上記の問題が解決されます。つまり、相関が高い変数を両方採用したい場合や、サンプルサイズ以上の変数を採用したい場合はelastic netを使用するのが良いということになります。
SCAD
以下のような正則化項をもつスパース推定法をSCADと言います。
λは正則化パラメータ、a(>2)は調整パラメータです*3。SCADはオラクル性質を満たしており、漸近正規性や変数選択の一致性を有します。
adaptive lasso
SCADはオラクル性質を満たしますが、複雑な正則化項の形をしており、lassoの形とは大きく異なります。そこで、lassoのL1ノルムの部分に適当な重みwを置くadaptive lassoが提案されています。
これはある正則条件の下でオラクル性質を満たします。重みwは既知でなければなりませんが、例えば最小二乗推定量βを使って、w=1/βと推定して用います。 wが大きい(βが小さい)とモデルに含む必要がないため回帰係数を0に推定する方向に持っていき、逆にwが小さい(βが大きい)と最小2乗推定量になるべく近い推定量を採用することを意味します。上記推定式は適切な式変形をすればlassoと同じアルゴリズムで推定値を得ることができますが、これについては以下のRでの数値計算例の時に簡単に説明します。
数値例
前回と同様にfitbitデータと自作の生活行動データを使用して、睡眠効率を目的変数としました。最小2乗推定、lasso、elastic net、SCAD、adaptive lassoそれぞれのアプローチで数値計算し結果を比較します。
データ*4
> 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, ...
標準化します。
y <- unlist(select(df, efficiency)) y <- y - mean(y) x <- scale(as.matrix(select(df, -efficiency)))
本来であれば、しっかりと可視化をしてデータの分布などを確認すべきですが割愛します。
lasso
lassoは前回の記事で紹介した通りです。
library("glmnet") fit <- glmnet(x, y, standardize = FALSE) autoplot(fit) + theme_pander() + scale_color_pander() # solution path set.seed(0) cvfit <- cv.glmnet(x, y, standardize = FALSE, nfolds = 10) coef(cvfit, s = "lambda.min") # 回帰係数
elastic net
glmnetパッケージのαパラメータに0~1の数値を代入すればよく(ここではα=0.5とした)、それ以外はlassoと全く同じです。
fit <- glmnet(x, y, alpha = 0.5, standardize = FALSE) autoplot(fit) + theme_pander() + scale_color_pander() # solution path set.seed(0) cvfit <- cv.glmnet(x, y, alpha = 0.5, standardize = FALSE, nfolds = 10) coef(cvfit, s = "lambda.min") # 回帰係数
lassoとほぼ同じ(と言うか全く同じ?)ですね。正直よくわかりません(やっぱりデータがよくないのかな…)。
SCAD
ncvregパッケージを使用します。
library("ncvreg") set.seed(0) fit <- ncvreg(x, y, family="gaussian", penalty = "SCAD") plot(fit) # solution path
影の部分はnonconvex region(非凸領域)です*5。 デフォルトのグラフでも綺麗ですが、各線がどの変数に紐づいているか不明なので、自作します。
lambda_df <- as.data.frame(fit$beta) # パラメータの抽出、data frame化 lambda_df$variable <- rownames(lambda_df) # tidy化 df <- lambda_df %>% gather(key="lambda", value="beta", -variable) %>% filter(variable != "(Intercept)") %>% # 切片変数を除外 mutate(lambda = as.double(lambda)) # 可視化 ggplot(df, aes(lambda, beta, color=variable)) + geom_line() + theme_pander() + scale_color_pander() + scale_x_reverse()
lassoやelastic netのsolution pathとは大きく異なりますが、出てくる変数(パラメータが非0となる変数)の順番は同じっぽいです。
cross validationによりパラメータλを選択し、回帰係数を算出します。
set.seed(0) cvfit <- cv.ncvreg(x, y ,family="gaussian", penalty = "SCAD") opt_lambda <- cvfit$lambda.min # 方程式を最小にするλ fit_lambda.min <- ncvreg(x, y, family="gaussian", penalty = "SCAD", lambda = opt_lambda) fit_lambda.min$beta # 回帰回数
adaptive lasso
adaptive lassoもlassoと同じglmnetパッケージを使用することができますが、少し工夫が必要です。以下のように式変形すると、
lassoのアルゴリズムが使用できます。以下、ざっとした計算の流れをまとめます。
- 重みwを最小2乗推定値により算出*6
- を算出・・・①
- glmnetによりパラメータ推定
によりパラメータを得ることができます。
W <- diag(1/abs(coef(lm(y~x))[-1])) xstar <- x %*% solve(W) # ①を算出 colnames(xstar) <- colnames(x) fit <- glmnet(xstar, y, standardize = FALSE) autoplot(fit) + theme_pander() + scale_color_pander() # solution path set.seed(0) cvfit <- cv.glmnet(xstar, y, standardize = FALSE, nfolds = 10) beta <- coef(cvfit, s = "lambda.min") solve(W) %*% beta[-1] # 回帰係数
※ を変換したに対するsolution pathです。参考までにやってみたものの、どのように解釈していいか不明です。
結果の解釈
推定結果をまとめた表が以下になります。
書籍には、
非凸正則化によるスパース推定法では、lassoに比べてより多くの回帰係数が0と推定されている。非0と推定されている部分に注目すると、これらの手法はlassoに比べて、最小2乗推定値に近い値が推定されていることがわかる。
と書いているが、ここではデータが微妙なせいか、よく分かりませんでした。
さいごに
今回のエントリでは、様々なスパース正則化法についての概要をまとめ数値例を示しました。*7
最後に個人的に気になった点、学びたい点をMEMOしておきます。
分からないことだらけですね。まだまだINPUTが足らないようなので、引き続き勉強頑張ります。
なお、本エントリは雰囲気で書いた部分も多いので、間違いがあれば遠慮なくご指摘くださいませ。
mikutaifuku
*1:時間の都合上(難しかったので)、計算アルゴリズムは飛ばしています。
*2:漸近正規性、変数選択の一致性を有する性質のことを言うらしいが、時間の都合上(難しかったので)詳細は書籍に譲ります
*3:aはリスク最小化の観点からa=3.7が用いられているらしいです。cross validationなどで選択することも可能です
*4:データ数少ないですね。最近サボり気味なのでしっかりデータ集めます。
*5:よくわからなかったので余力があればしっかりと調べます
*6:リッジでも問題ないが、リッジ推定を用いた場合はオラクル性質を満たさないらしい
*7:疲れて最後の方は少し適当になってしまった
*8:adaptive lassoでもある条件の下だけらしいが