mikutaifukuの雑記帳

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

臨床予測モデルの開発に必要なサンプルサイズの算出(Calculating the sample size required for developing a clinical prediction model)

概要

ヘルスケア領域のデータ分析に携わっていると、よく質問されるのが

どれくらいサンプルサイズが必要なの?

という質問です。統計的仮説検定であれば、

有意水準、検出力…等を指定して、必要サンプル数を求めることが可能だと思います

とそれらしく回答しているのですが、予測モデルを構築する場合は、

なるべく多くのサンプル集める。特徴量の次元数の10倍のデータが必要みたいなことは、経験則として提示されていますね…

というような感じで、もごもごと回答しています。

そこで、google検索してみたところ「Calculating the sample size required for developing a clinical prediction model」という論文が存在していたようなので、流し読みしてみました。

www.bmj.com

この論文では、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)

f:id:mikutaifuku:20211015232738p:plain

さいごに

  • clinical prediction model(臨床予測モデル)、しかも0/1の2値のロジスティック回帰モデルに限定しての話で、おそらく一般的な機械学習だとこの限りではないと思います。
  • また、論文内では、train/testのようなデータの分割(ホールドアウト検証)は推奨されておらず、ブートストラップ法による検証が推奨されており、このあたりの感覚は少し違和感があります。
  • 雰囲気でブログを書いただけなので、おかしな事をやっていたり、間違っていたらこっそりご連絡ください。

参考

Calculating the sample size required for developing a clinical prediction model | The BMJ

CRAN - Package pmsampsize

MENTA登録してみました

知り合いが登録しているとのことで、チャレンジしてみることにしました。

【初心者〜中級者向け】R/Pythonでのデータ分析スキル習得に向けてサポートします! mikutaifuku|メンターに教えてもらおう

こういうのをきっかけに、ちょっとした弱い繋がりができればと期待しつつ、がんばってみます。

ggplot2でfacetごとのヒストグラムに色んな情報を付与する② 〜Tidy evalによる関数化と繰り返し処理〜

概要

前回の続きで、自分用メモです。

mikutaifuku.hatenablog.com

タイトルの通り、Tidy evalにより関数を作り*1、全ての変数に対して処理を繰り返し行います。 要は、以下のようなグラフを Sepal.Width, Petal.Length, Petal.Width 全てに対して作ります。

f:id:mikutaifuku:20190511035424p:plain

なお、前提知識として、ブログ下部の「参考」にある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"))
    
}

参考

tidyeval.tidyverse.org

さいごに

上記のような処理をする際に、Tidy evalなるものは役立つのではないかと勝手に思っています。 雰囲気で作成したので、間違っていたら、遠慮なくご指摘ください。 また、他に良いやり方があれば、教えてください。

*1:この表現は適切でないかもしれませんが、ご容赦ください

*2:下記のようにループ処理をしなくても、うまく出来る方法があるはず…?

ggplot2でfacetごとのヒストグラムに色んな情報を付与する① 〜基本作図編〜

概要

大した内容ではないですが、仕事で使ったので自分用メモです。なので、やや適当です。

以下のようなグラフを作ります。

f:id:mikutaifuku:20190511035424p:plain

具体的には、以下のような特徴を備えたグラフです。

  • 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)

f:id:mikutaifuku:20190511035424p:plain

参考

いつも有益な情報を提供してくださる、株式会社ホクソエムのブログを参考にしました。 blog.hoxo-m.com

テキストがグラフ内におさまる方法を参考にしました。

ggplot(df, aes(x, y)) +
  geom_text(aes(label = text), vjust = "inward", hjust = "inward")

ggplot2.tidyverse.org

さいごに

パート②では、加工処理とプロット処理を関数化し、プロットを全ての変数に対して繰り返し行います。

雰囲気で作成したので、間違っていたら、遠慮なくご指摘ください。

Fitbitデータでスパース推定に入門してみる④ 〜Graphical lassoで変数間の関係性について調べる〜

はじめに

これまで述べてきたのは、主にモデルの回帰係数の推定に基づく話でした。

mikutaifuku.hatenablog.com

mikutaifuku.hatenablog.com

mikutaifuku.hatenablog.com

今回のエントリでは、これまでの回帰係数の推定ではなく、共分散構造に注目した手法を取り上げます。多数の変数間の依存関係を推定するために、ガウシアングラフィカルモデルにL1正則化の考え方を応用する、Graphical lassoとよばれる手法を使います。 なお、今回の記事も今までと同様に以下の書籍(5.1 グラフィカルモデルにおけるスパース推定)を参考に書いたものですので、詳細は書籍を参考にしていただければと思います。

スパース推定法による統計モデリング (統計学One Point)

スパース推定法による統計モデリング (統計学One Point)

ガウシアングラフィカルモデル

 X = (X_1, ..., X_p)^{T}が多変量正規分布 N( \boldsymbol{\mu} ,\Sigma) に従うとき、変数  X_j X_kの関係性を推定し、それをグラフ化する手法です。分散共分散行列の逆行列 \Omega = \Sigma^{-1}とし、\Omega (j,  k)成分を w_{jk}とします。詳細は書籍に譲りますが、 w_{jk} = w_{kj} = 0の時に、 X_j X_kは他の変数を与えた元で条件付き独立と言います。

データ分析でよくみられる相関係数は、他の変数からの影響も含めて測っているのに対して、行列 \Omegaは他の変数からの影響を除いて、(すなわち X_j X_k以外の変数の条件のもとで)関連の強さを測るものになります。

例えば、以下のような行列 \Omega(左図)が得られている場合は、右のようなグラフで条件付き独立性の有無を表現することができます*1

f:id:mikutaifuku:20180323152702p:plain

このグラフが表現するのは「条件付き独立性」があるか否かということだけであって、「変数間の関係の向き」を推定することができません。 もし、グラフの辺の向きを推定したい場合は、以前紹介したようにLingamモデルを仮定することで推定することができます。

mikutaifuku.hatenablog.com

パラメータ推定

グラフの辺の数は最大で頂点の2乗のオーダー程度になりますが、実際に辺が存在する数は少ないという意味ではスパースと言えます。データからスパースな分散共分散行列の逆行列(精度行列) \Omega を推定する手法をgraphical lassoと言います。 データ  \boldsymbol{x_1} , … ,\boldsymbol{x_n} が観測され、これらが正規分布に従っていると仮定します。この時、標本から計算される分散共分散行列をSとすると、対数尤度関数に基づく損失関数は、

f:id:mikutaifuku:20180323172632p:plain:w220

となります(定数省略)*2。これに L_1正則化 \Sigma_{j,k} |w_{jk}|を加えた式を最小化にすることで、精度行列\Omega を推定します。

f:id:mikutaifuku:20180323173357p:plain:w250

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。コードはこのブログの最後に載せています。

f:id:mikutaifuku:20180323212010p:plain

  • 歩数と消費カロリー(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の時のグラフィカルモデルです

f:id:mikutaifuku:20180323222606p:plain

λ=0.4の時のグラフィカルモデルです。 f:id:mikutaifuku:20180323222627p:plain

  • 安静時心拍数(RestingHeartRate)はどの変数とも関係がなさそうです。
  • 歩数(steps)と消費カロリー(caloriesOut)は関係がありそうです。

推定された\Omega から、偏相関係数という量を算出して、可視化することもできますが、次回の機会にとっておきます。

相関行列の可視化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データに適用することでそれぞれの手法を比較しました。

mikutaifuku.hatenablog.com

今回のエントリでは、少し趣向を変えて心拍数データに対して、fused lassoと呼ばれる隣合う変数の関係性を加味した手法使ってみます。なお、今回のエントリは、書籍の第3章のfused lassoの部分を自分なりに整理してまとめただけなので、丁寧に学びたい方は、書籍を読んでいただければと思います。

スパース推定法による統計モデリング (統計学One Point)

スパース推定法による統計モデリング (統計学One Point)

fused lassoの概要

fused lassoは以下のような性質をもつデータに対して提案された手法です。

  • p個の説明変数に順序関係がある
  • 隣接する説明変数のいくつかは目的変数に同程度の寄与がある

最小二乗推定量を算出する式の後ろに、正則化項にL1ノルムと隣接する係数の差の絶対値を考えたものです。

f:id:mikutaifuku:20180315135856p:plain:w300

この式を最小にすることでパラメータβを推定します。正則化項として係数の差の絶対値を考えているため、推定値が等しくなるように働くという仕組みです。

時系列データへの応用

Fitbitの睡眠時心拍数データ*1を使って、時間的に変化する値の推定を考えます。隣り合う時点の値はほぼ等しく、ときどきある時点でジャンプしています。

f:id:mikutaifuku:20180315133942p:plain

時刻tでの観測値を y_{t}とし、

f:id:mikutaifuku:20180315152913p:plain:w200

モデリングします。εは観測ノイズです。隣り合う時点でβの値が変わりにくいという性質を、正則化

f:id:mikutaifuku:20180315153127p:plain:w210

で表現します。この項の効果で不連続に飛ぶ点の数を少なくしながら推定を行います。つまり、以下のような式を最小にする事で推定パラメータを求めます。

f:id:mikutaifuku:20180315151614p:plain:w220

これは、上記の(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

f:id:mikutaifuku:20180315140924p:plain

正則化パラメータλを変えた時のフィッティングを折れ線で示しています。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")

f:id:mikutaifuku:20180315144910p:plain

そこそこ綺麗なグラフが完成しました。λが大きいと直線に近づき、λが小さいと観測値に沿った折れ線が得られていることがわかります。λ=4のグラフを見ると4:00や6:00頃の変化点をしっかりと検出できています。

さいごに

今回のエントリでは、Fitbitの心拍数データに対してfused lassoを適用し、変化点の検出を試みました。ただ、実務においては、このように全てのデータが揃ってから変化点を検出しても遅く、要望として例えば「異常を5分前に予兆したい」みたいなことが挙げられます。が、この手法では残念ながら解決できません(おそらく)。という事で、これら課題をどうやったら解決できそうか、うまく問題設定を考えてチャレンジしてみようと思います。

このブログではパッケージを駆使して「やってみた」系の内容がほとんどです(てか全部)。自分でしっかり頭を使ってモデリングをしている訳でもなければ、自分で数式をしっかり理解した上でアルゴリズムを実装している訳でもありません。この状況は個人的にはなかなかヤバイと思いながらも、時間コストを考えるとこれが限界な気がしなくもありません。色々難しいですね、ちょっと色々考えてみます。

引き続き頑張ります。

mikutaifuku

*1:データのばらつきが小さい睡眠時のデータを使います。

*2:書籍にはλ1 > 0と書いてたが誤植かな?

Fitbitデータでスパース推定に入門してみる② 〜lasso, elastic net, SCAD, adaptive lassoの概要とRでの数値例〜

はじめに

前回の記事では、スパース推定の概要と代表的な手法であるlassoについてまとめ、fitbitデータに適用することで、睡眠効率に影響を与える要因を探索しました。

mikutaifuku.hatenablog.com

今回のエントリでは、lasso以外の様々なスパース正則化法を紹介し、前回と同様にFitbitデータに適用することで、それぞれの手法を比較したいと思います。なお、今回のエントリは、以下の書籍の第1・2章を自分なりに簡単にまとめた*1だけなので、丁寧に学びたい方は、書籍を読んでいただければと思います。

スパース推定法による統計モデリング (統計学One Point)

スパース推定法による統計モデリング (統計学One Point)

lassoの概要と弱点

概要

lassoは以下のように、最小二乗推定量を算出する式の後ろに、正則化項にL1ノルムを考えたものです。 f:id:mikutaifuku:20180302165849p:plain 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乗推定量を算出する式の後ろに、下記のような正則化項を加えたものを最小化することで、パラメータを推定します。

f:id:mikutaifuku:20180313171548p:plain:w300

これはL1ノルムとL2ノルムを合わせた形になっており、α=1の時はlasso, α=0の時はリッジとなります。このような正則化項を取ることで、上記の問題が解決されます。つまり、相関が高い変数を両方採用したい場合や、サンプルサイズ以上の変数を採用したい場合はelastic netを使用するのが良いということになります。

SCAD

以下のような正則化項をもつスパース推定法をSCADと言います。

f:id:mikutaifuku:20180313171929p:plain:w300

λは正則化パラメータ、a(>2)は調整パラメータです*3。SCADはオラクル性質を満たしており、漸近正規性や変数選択の一致性を有します。

adaptive lasso

SCADはオラクル性質を満たしますが、複雑な正則化項の形をしており、lassoの形とは大きく異なります。そこで、lassoのL1ノルムの部分に適当な重みwを置くadaptive lassoが提案されています。

f:id:mikutaifuku:20180313172145p:plain:w200

これはある正則条件の下でオラクル性質を満たします。重み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") # 回帰係数

f:id:mikutaifuku:20180313181455p:plain

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") # 回帰係数

f:id:mikutaifuku:20180313181614p:plain

lassoとほぼ同じ(と言うか全く同じ?)ですね。正直よくわかりません(やっぱりデータがよくないのかな…)。

SCAD

ncvregパッケージを使用します。

library("ncvreg")
set.seed(0)
fit <- ncvreg(x, y, family="gaussian", penalty = "SCAD")
plot(fit) # solution path

f:id:mikutaifuku:20180313181818p:plain

影の部分は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()

f:id:mikutaifuku:20180313182238p:plain

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パッケージを使用することができますが、少し工夫が必要です。以下のように式変形すると、

f:id:mikutaifuku:20180313222128p:plain:w400

lassoのアルゴリズムが使用できます。以下、ざっとした計算の流れをまとめます。

  • 重みwを最小2乗推定値により算出*6
  • X^{*} = XW^{-1}を算出・・・①
  • glmnetによりパラメータ推定
  •  \hat{\boldsymbol{\beta}} = W^{-1} \hat{\boldsymbol{\eta}}

によりパラメータを得ることができます。

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] # 回帰係数

f:id:mikutaifuku:20180313185001p:plain

Xを変換したX^{*}に対するsolution pathです。参考までにやってみたものの、どのように解釈していいか不明です。

結果の解釈

推定結果をまとめた表が以下になります。

f:id:mikutaifuku:20180313184548p:plain

書籍には、

非凸正則化によるスパース推定法では、lassoに比べてより多くの回帰係数が0と推定されている。非0と推定されている部分に注目すると、これらの手法はlassoに比べて、最小2乗推定値に近い値が推定されていることがわかる。

と書いているが、ここではデータが微妙なせいか、よく分かりませんでした。

さいごに

今回のエントリでは、様々なスパース正則化法についての概要をまとめ数値例を示しました。*7

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

  • SCADの活躍場所。オラクル性質を満たす正則化項を採用したいのならadaptive lassoを使えば良いということでOKか。*8
  • 計算アルゴリズムも余力があればしっかりと理解したい

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

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

mikutaifuku

*1:時間の都合上(難しかったので)、計算アルゴリズムは飛ばしています。

*2:漸近正規性、変数選択の一致性を有する性質のことを言うらしいが、時間の都合上(難しかったので)詳細は書籍に譲ります

*3:aはリスク最小化の観点からa=3.7が用いられているらしいです。cross validationなどで選択することも可能です

*4:データ数少ないですね。最近サボり気味なのでしっかりデータ集めます。

*5:よくわからなかったので余力があればしっかりと調べます

*6:リッジでも問題ないが、リッジ推定を用いた場合はオラクル性質を満たさないらしい

*7:疲れて最後の方は少し適当になってしまった

*8:adaptive lassoでもある条件の下だけらしいが