SNAGeek

Sociologically Technological, and Technologically Sociological

ggplot2で多重クロス表を作る

  • 最近は大学院の研究でひたすら3重クロス表の分析を行っています。
  • Rでは xtable パッケージの xtable() を使うと2重クロス表をTeX形式で出力してくれるのですが、3重以上の多重クロス表には対応していません。
    • gridExtra::grid.table() を使う方法もあるかと思いますが、余白の扱いや多重の場合が少し面倒な印象があります。
  • そこで今回は、ggplot2を使った多重クロス表の出力方法を紹介します

セッション情報

  • 毎度のことながらrstudio cloud上で動かしています。
> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 LTS

Matrix products: default
BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C           LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

下準備

{MASS} パッケージ内の Cars93 データを用います。

library(tidyverse)
data <- MASS::Cars93

まずは普通に出力

  • 車種(Type)、駆動系(DriveTrain)、生産国(Origin)の3重クロス表を作成します。
crosstab <- data %>% select(Type,DriveTrain,Origin) %>% table
> crosstab
, , Origin = USA

         DriveTrain
Type      4WD Front Rear
  Compact   0     7    0
  Large     0     7    4
  Midsize   0     9    1
  Small     0     7    0
  Sporty    2     2    4
  Van       3     2    0

, , Origin = non-USA

         DriveTrain
Type      4WD Front Rear
  Compact   1     6    2
  Large     0     0    0
  Midsize   0     8    4
  Small     2    12    0
  Sporty    0     5    1
  Van       2     2    0

geom_tileを使ってplotする

  • 何はともあれtidyに。
tidy_crosstab <- as.data.frame.table(crosstab)
> tidy_crosstab
      Type DriveTrain  Origin Freq
1  Compact        4WD     USA    0
2    Large        4WD     USA    0
3  Midsize        4WD     USA    0
4    Small        4WD     USA    0
5   Sporty        4WD     USA    2
6      Van        4WD     USA    3
7  Compact      Front     USA    7
8    Large      Front     USA    7
9  Midsize      Front     USA    9
10   Small      Front     USA    7
11  Sporty      Front     USA    2
12     Van      Front     USA    2
13 Compact       Rear     USA    0
14   Large       Rear     USA    4
15 Midsize       Rear     USA    1
16   Small       Rear     USA    0
17  Sporty       Rear     USA    4
18     Van       Rear     USA    0
19 Compact        4WD non-USA    1
20   Large        4WD non-USA    0
21 Midsize        4WD non-USA    0
22   Small        4WD non-USA    2
23  Sporty        4WD non-USA    0
24     Van        4WD non-USA    2
25 Compact      Front non-USA    6
26   Large      Front non-USA    0
27 Midsize      Front non-USA    8
28   Small      Front non-USA   12
29  Sporty      Front non-USA    5
30     Van      Front non-USA    2
31 Compact       Rear non-USA    2
32   Large       Rear non-USA    0
33 Midsize       Rear non-USA    4
34   Small       Rear non-USA    0
35  Sporty       Rear non-USA    1
36     Van       Rear non-USA    0
  • ggplotにつっこむ。
g <- ggplot(data = tidy_crosstab, aes(x = Type, y = DriveTrain, label = Freq))
g <- g + geom_tile(fill = "white",col = "black")
g <- g + geom_text()
g <- g + facet_wrap(~Origin, nrow = 1)
  • 出力されるのはこんな感じの画像になります。

f:id:meana0:20200310220153p:plain

4重クロス表

  • AirBagという変数を増やした4重クロス表を出力します
  • 今度は、facet_grid() を使います。
crosstab2 <- data %>% select(Type,DriveTrain,AirBags,Origin) %>% table
tidy_crosstab2 <- as.data.frame.table(crosstab2)
g <- ggplot(data = tidy_crosstab2, aes(x = Type, y = DriveTrain, label = Freq))
g <- g + geom_tile(fill = "white",col = "black")
g <- g + geom_text()
g <- g + facet_grid(Origin ~ AirBags)

f:id:meana0:20200310222149p:plain

{gnm}を使ったクロス表の分析

  • 最近は対数線形・対数乗法モデルの復習をしています。
  • 色々と読み返してはいるのですが、以下の本は簡潔にまとまっており、非常に勉強になりました。
  • オッズ比から始まって、RC(M)-L modelまで取り扱っています。

Association Models (Quantitative Applications in the Social Sciences) by Raymond Sin-Kwok Wong(2010-02-09)

Association Models (Quantitative Applications in the Social Sciences) by Raymond Sin-Kwok Wong(2010-02-09)

  • 作者:Raymond Sin-Kwok Wong
  • 出版社/メーカー: SAGE Publications, Inc
  • 発売日: 1787
  • メディア: ペーパーバック

  • 今回は備忘録も兼ねて、同書のExample2.2をRの{gnm}パッケージを使って一部再現しようと思います。
    • Rでの対数線形・対数乗法モデルの実装方法に関しては、あまりまとまった日本語の資料がないようでした。
    • モデル自体の解説は日本語でもすぐに資料が見つかると思います。

github.com

データの用意

  • 同書に記載ある通り、学歴(EDUC)と現職(OCC)のクロス表 ex.2.2.txt を用意します。
  • 元ネタはGSS(1985-1990 cumulative, 女性のみ)のようです。
,上層ノンマニュアル,下層ノンマニュアル,上層マニュアル,下層マニュアル,農業
大卒以上,518,95,6,35,5
短大,81,67,4,49,2
高校,452,1003,67,630,5
高校未満,71,157,37,562,12

実行環境

> sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8    
 [5] LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8    LC_PAPER=C.UTF-8       LC_NAME=C             
 [9] LC_ADDRESS=C           LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] purrr_0.3.3 dplyr_0.8.3 gnm_1.1-0  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1       lattice_0.20-35  crayon_1.3.4     relimp_1.0-5     assertthat_0.2.1 MASS_7.3-51.5   
 [7] grid_3.4.4       R6_2.4.0         magrittr_1.5     pillar_1.4.2     rlang_0.4.2      qvcalc_1.0.1    
[13] rstudioapi_0.10  Matrix_1.2-18    tools_3.4.4      glue_1.3.1       compiler_3.4.4   pkgconfig_2.0.2 
[19] tidyselect_0.2.5 nnet_7.3-12      tibble_2.1.3    

ライブラリ読み込み

library(gnm)
library(dplyr)
library(purrr)

前処理

  • 何はともあれ、行ラベル・列ラベル・頻度からなるtidyなデータフレームにします。
tab <- read.table("ex.2.2.txt",header = TRUE,sep = ",",row.names = 1)
mat <- as.matrix(tab)
df <- as.data.frame.table(mat)
colnames(df) <- c("EDUC","OCC","Freq")
> df
       EDUC                OCC Freq
1  大卒以上 上層ノンマニュアル  518
2      短大 上層ノンマニュアル   81
3      高校 上層ノンマニュアル  452
4  高校未満 上層ノンマニュアル   71
5  大卒以上 下層ノンマニュアル   95
6      短大 下層ノンマニュアル   67
7      高校 下層ノンマニュアル 1003
8  高校未満 下層ノンマニュアル  157
9  大卒以上     上層マニュアル    6
10     短大     上層マニュアル    4
11     高校     上層マニュアル   67
12 高校未満     上層マニュアル   37
13 大卒以上     下層マニュアル   35
14     短大     下層マニュアル   49
15     高校     下層マニュアル  630
16 高校未満     下層マニュアル  562
17 大卒以上               農業    5
18     短大               農業    2
19     高校               農業    5
20 高校未満               農業   12

分析

  • まずはクロス表全体としてオッズ比が1になる independent model を作ります。
  • 一般化非線形モデルを推定するためのパッケージ{gnm}を使います。
> gnm(formula = Freq ~ EDUC + OCC, data = df, family = poisson)

Call:
gnm(formula = Freq ~ EDUC + OCC, family = poisson, data = df)

Coefficients:
          (Intercept)               EDUC短大               EDUC高校           EDUC高校未満  
               5.2557                -1.1775                 1.1858                 0.2415  
OCC下層ノンマニュアル      OCC上層マニュアル      OCC下層マニュアル                OCC農業  
               0.1640                -2.2867                 0.1286                -3.8448  

Deviance:            1373.176 
Pearson chi-squared: 1455.658 
Residual df:         12 
  • ここから階層的にモデルを作っていきます。
  • RモデルやCモデルを作るためには、推定の問題上、パラメータを標準化し、足して0になるように制約をかける必要があります。
Row <- scale(as.numeric(row(mat)), scale = TRUE)
Col <- scale(as.numeric(col(mat)), scale = TRUE)
  • 本来であればR+C+RCモデルなども記載があるのですが、パラメータ制約がややこしいので省略しています。
  • あとで map() するため、結果をlistに放り込んでいきます。
set.seed(1024)
gnm_list <- list()
gnm_list$O <- gnm(formula = Freq ~ EDUC + OCC, data = df, family = poisson)
gnm_list$U <- update(gnm_list$O, ~ . + Row:Col)
gnm_list$R <- update(gnm_list$O, ~ . + EDUC:Col)
gnm_list$C <- update(gnm_list$O, ~ . + OCC:Row)
gnm_list$RpC <- update(gnm_list$O, ~ . + EDUC:Col + OCC:Row) #R+Cモデル
gnm_list$RC <- update(gnm_list$O, ~ . + Mult(EDUC,OCC))
gnm_list$RC2 <- update(gnm_list$RC, ~ . + Mult(EDUC,OCC,inst = 2))
  • さて、無事に各モデルが推定できたはずです。

モデルの適合度

  • 次は各モデルの適合度を見ていきます。どのモデルが最もクロス表をよく説明するでしょうか。
  • 自由度df、逸脱度G2(-2log尤度)、p値(Χ2検定)、BICの4つを出力します。
fit_df <- function(gnm_object){
  df <- gnm_object$df.residual
  g_sq <- gnm_object$deviance
  p_value <- chi2.p.value(gnm_object)
  bic <- g_sq - gnm_object$df.residual * log(sum(gnm_object$data$Freq))
  fit_df <- data.frame(df,g_sq,p_value,bic)
  return(fit_df)
}

gof <- bind_rows(map(gnm_list,fit_df),.id = "model")
  • 結果を見てみます。
> gof
  model df         g_sq  p_value        bic
1     O 12 1373.1757657 0.000000 1274.08092
2     U 11  244.0177150 0.000000  153.18077
3     R  9  205.9735836 0.000000  131.65245
4     C  8  155.3722685 0.000000   89.30903
5   RpC  6   91.6077617 0.000000   42.06034
6    RC  6  125.0597486 0.000000   75.51232
7   RC2  2    0.6001022 0.750266  -15.91571
  • 最後のRC(2)モデルが最もBICが低く、かつp値も大きいため、このモデルが選択されます。

    • p値が大きいというのは、Χ2値が小さい、つまり、RC2モデルによって求められる各セルの期待度数が、元データからそれほど乖離していないことを示しています。
    • 本来は、自由度をいくつ使ってG2が改善したかなども併せて見ていく必要がありますが…。
  • RC(2)モデルのパラメータはこのような感じになります。

> gnm_list$RC2

Call:
gnm(formula = Freq ~ EDUC + OCC + Mult(EDUC, OCC) + Mult(EDUC, 
    OCC, inst = 2), family = poisson, data = df)

Coefficients:
                                  (Intercept)                                       EDUC短大  
                                     5.066410                                      -0.916769  
                                     EDUC高校                                   EDUC高校未満  
                                     1.067655                                       0.217877  
                        OCC下層ノンマニュアル                              OCC上層マニュアル  
                                     0.353842                                      -2.143860  
                            OCC下層マニュアル                                        OCC農業  
                                     0.129995                                      -3.973502  
                    Mult(., OCC).EDUC大卒以上                          Mult(., OCC).EDUC短大  
                                    -5.581051                                      -1.462538  
                        Mult(., OCC).EDUC高校                      Mult(., OCC).EDUC高校未満  
                                     0.975465                                       1.973569  
          Mult(EDUC, .).OCC上層ノンマニュアル            Mult(EDUC, .).OCC下層ノンマニュアル  
                                    -0.249405                                       0.099685  
              Mult(EDUC, .).OCC上層マニュアル                Mult(EDUC, .).OCC下層マニュアル  
                                     0.222067                                       0.314239  
                        Mult(EDUC, .).OCC農業            Mult(., OCC, inst = 2).EDUC大卒以上  
                                    -0.013403                                       1.240632  
              Mult(., OCC, inst = 2).EDUC短大                Mult(., OCC, inst = 2).EDUC高校  
                                     0.683454                                      -1.323081  
          Mult(., OCC, inst = 2).EDUC高校未満  Mult(EDUC, ., inst = 2).OCC上層ノンマニュアル  
                                     3.154406                                      -0.168216  
Mult(EDUC, ., inst = 2).OCC下層ノンマニュアル      Mult(EDUC, ., inst = 2).OCC上層マニュアル  
                                    -0.246388                                       0.007613  
    Mult(EDUC, ., inst = 2).OCC下層マニュアル                Mult(EDUC, ., inst = 2).OCC農業  
                                     0.094214                                       0.382994  

Deviance:            0.6001022 
Pearson chi-squared: 0.5746549 
Residual df:         2 
  • RC effectの1次元目は学歴・職業階層の一次元的な”序列”を表していそうです。
    • つまり、高い(低い)学歴が高い(低い)職業的地位と結びついている
  • 一方で2次元目は解釈が難しく、筆者によればこれはメリトクラシーでは説明できない、労働市場の参入障壁を表しているといいます。
    • RC(M)モデルの Mを増やせば説明力は上がりますが、その解釈は分析者の腕の見せどころですね…。

xtableのエスケープを抑制する

xtableを使うとmatrixやdataframeをTeX形式の表に変換してくれて便利です。 しかし、表の中に数式などが含まれている場合、xtableの出力結果の中で _{ の前にエスケープ文字列 \ がつく場合があります。 当然、出力結果を貼り付けてそのままコンパイルすると、適切に数式が表示されません。

> library(tidyverse)
> library(xtable)
> 
> d <- tibble(Variable = c("$\\hat{X}$"),
+             Value = "0.1")
> xtable(d)
% latex table generated in R 3.5.3 by xtable 1.8-4 package
% Tue Jan 28 14:28:26 2020
\begin{table}[ht]
\centering
\begin{tabular}{rll}
  \hline
 & Variable & Value \\ 
  \hline
1 & \$$\backslash$hat\{X\}\$ & 0.1 \\ 
   \hline
\end{tabular}
\end{table}

xtableオブジェクトをprint関数に渡して以下のようにオプションを指定すると、エスケープしません。これは複数指定できます。

#列ラベルのエスケープを抑制したい場合
xtable(d) %>% print(sanitize.colnames.function = identity)
#行ラベルのエスケープを抑制したい場合
xtable(d) %>% print(sanitize.rownames.function = identity)
#各要素のエスケープを抑制したい場合
xtable(d) %>% print(sanitize.text.function = identity)

上の例で実際にやってみると、確かにエスケープされなくなります。

xtable(d) %>% print(sanitize.text.function = identity)
% latex table generated in R 3.5.3 by xtable 1.8-4 package
% Tue Jan 28 14:28:59 2020
\begin{table}[ht]
\centering
\begin{tabular}{rll}
  \hline
 & Variable & Value \\ 
  \hline
1 & $\hat{X}$ & 0.1 \\ 
   \hline
\end{tabular}
\end{table}

参考:

print.xtable function | R Documentation

英文校正から自分の癖を知るシリーズ1

先日、あるオープンアクセス誌に英語論文を投稿しました。 後学のために、英文校正会社から返ってきた修正点のうち、自分がやらかしがちな部分を洗い出してみました。 同じような英語ライティングの癖を持っている人もいるかもしれないので、一般に公開してみることにします。基礎的なミスが多くて恥ずかしいですが...。 ちなみにイギリス英語の準拠して直してもらうようにオーダーしています。
 
基本的に過去形でいい
 
「(変数を)統制する」の意味でcontrolを使う時はforが必要
 
「日本の地方部で」はin rural area in Japanではなくin rural Japanで表せる
 
~was found to be...は単純にwasでよい
 
「に関する非対称性」のような時に使う「関する」はaboutではなくregarding
 
「観測できなかった」はcannot be observedではなくwas not observedでいい
 
APAの場合キーワードは頭文字が大文字
 
APAの場合ショートタイトルはすべて大文字。
 
It has been widely discussed that~は冗長
 
i.e.で言い換える時はカッコ内に書く
 
「~に対する偏見」はprejudice toではなくprejudice toward
 
「大きく依存する」はdeeply dependではなくdepends greatly
 
「~の黎明期」はat the early stage of~ではなくin the early stages of~
 
「~への調査」はresearch toではなくresearch on (surveyも同様)
 
「学部生」はungraduated studentsではなくundergraduate students
 
タイミング的に同時の場合はat the same timeではなくsimultaneously
 
「~を検討している研究」はstudies which examine...ではなくstudies examining...
 
もう一つ例を出す時はFor another exampleよりはIn another exampleのほうがベタ
 
examineとかevaluateだけじゃなくてassessも使えるようになりたい
 
長い文でも時制を保ちたい
 
electoronicではなくelectronic
 
「○○(学術用語)は~を表す」は~expresses...ではなくis an expression of...を使う
 
やたらとmean使わなくてもbe動詞でいい
 
settingsを使えるようになりたい
 
This studyの類義語としてThe current researchも使える
 
いろんな要素の間の違いを言う時はdifference amongではなくdifference across
 
いちいちtry to revealとか書かない。examineでいい。
 
色々列挙して他にも要素が考えられる場合はand moreで締める
 
In this perspectiveよりもFrom this perspectiveがベタ
 
「だけで引き起こされるわけではない」はbe not caused by only~ではなくbe not caused merely by~
 
「AとBは区別されなければならない」はmust be distinguished from one another
 
APAの場合、著者名にアクサンなどが入っていても本文中ではアルファベットで表記する
 
most famousよりもwell-knownのほうがこなれている
 
informationは不可算なのでthese informationとか言わない
 
「~を質問する」はaskじゃなくてask for
 
「後述の分析」はthe aforementhioned analysis
 
「~の問題に対応する」という意味ならto handle withじゃなくてto handleでいい
 
giveは口語なのでprovide等を使う
 
「~に同意する」なのにagree withを使わない。agree to。
 
「丸形の~」はround shapeではなくcircle shaped
 
読者の興味を引きたい時はEspeciallyではなくWorth noting,から始める
 
「言うなれば」の意味でnamely,を使いがちだが、要らない
 
対比を強調したい時はイタリック
 
アルゴリズムを採用する」はemploy algorithm
 
「この研究にはいくつかの限界がある」はThere are some limitations to this study. on this studyではない。
 
「たった一例」ということを強調したければ単にaではなくa singleで形容する
 
APAの場合、編著本の中の章はページ数を(pp. 1-14)のようにppをつける
 
APAの場合、論文名は最初だけ大文字
 

RでERGMを実装したかった(が、失敗した)

このエントリは、Sansan Advent Calendar 2019 21日目の記事です。

adventar.org

  • 基本的にタイトルの通りですが、この記事では統計的ネットワーク分析のデファクトスタンダードとなっているERGMをRで実装していきます。
    • {igraph}以外のパッケージは使わずにできるだけスクラッチで開発します。
  • 前もって断っておきますが、正しくパラメータ推定できるような実装には至ることができませんでした。

ERGMとは

  • ERGM = Exponential Random Graph Model
    • 日本語では「指数ランダムグラフモデル」と訳されることが多いです。
    • 略称は「あーぐむ」と読むようです。
  • 概要については以下の記事で紹介されています。

qiita.com

ERGMを使うと何が嬉しいのか

  • 観測されたネットワークが、どのような構造的なメカニズムによって生成されたのかを知ることができます。
  • ERGMの最大の強みは、ダイアド間の相互依存性を考慮できることにあります。
    • ダイアドを1つのケースとして、エッジの有無を従属変数としたロジスティック回帰も可能ですが、これは残差の独立性の仮定に反してしまいます。
  • 社会ネットワークにおいては、下の図のように、相互的なつながりが発生する傾向(相互性)や、友達同士は友達になりやすい傾向(トライアド閉鎖)がありますが、ERGMはダイアド同士の相互依存性を仮定できるので、このようなネットワーク内生的な効果を推定することができます。
    • 例えば、トライアド閉鎖は、トライアングルの数などのネットワーク統計量をモデルに投入して、その効果を確認することができます。

f:id:meana0:20191220140649p:plain

  • ERGMは一般に以下のような指数形式で定式化されます。
 \displaystyle
\Pr(\bf{Y}=\bf{y}) = \frac{1}{\kappa(\bf{y})} \exp\{\sum_A \theta_k g_k (\bf{y})\}
 \displaystyle
\kappa(\bf{y}) = \sum_{y' \in {Y}} \exp\{\theta^T g(y')\}.
  •  \bf {Y} : ネットワークの確率変数。
  •  \bf {y} : ネットワークの実現値(実際に観測されたネットワーク)。
  •  g_k (\bf {y}) : ネットワーク形象  k のネットワーク統計量。
    • 例えば、 y に4つのトライアングルが含まれている場合、 g_ {triangle}(\bf{k})= 4
  •   \theta_k : ネットワーク形象  k に対応するパラメータ。
  •  A :モデル内で検討するトライアングルなどのネットワーク形象の集合。
  •  \kappa(\bf{y}) : 正規化定数。可能なすべてのネットワークの  \theta_k z_k(y) の合計値。

    • N=20ほどの小規模ネットワークでも現実的に計算できない数の膨大なネットワークを検討しなければならないため、 \kappa(\bf{y}) の計算を回避するためのテクニックが色々と考案されています。
  • 例えば、ERGMで以下のようにパラメータが推定されたとします。

Pr=exp(-0.5*エッジ数+0.5*閉鎖トライアドの数+0.2*同じ性別間のつながりの数) 
  • この時、観測されたネットワークにおいて「友達同士が友達になりやすい」(トライアド閉鎖)という傾向と「同じ性別同士でつながりやすい」(ホモフィリー)という傾向が同時に作用していると解釈することができます。

  • ERGMを用いた研究事例は色々あるのですが、Goodreau et al. (2009)は、アメリカの学校内でのホモフィリーを、トライアド閉鎖の効果を統制しながら推定しています。

www.ncbi.nlm.nih.gov

  • ここからは、{ergm}パッケージを用いたERGMの推定と、RによるERGMパラメータ推定の実装の一例について紹介します。

セッション情報

  • RStudio Cloud上で動かしています。
> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 LTS

Matrix products: default
BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8   
 [6] LC_MESSAGES=C.UTF-8    LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C           LC_TELEPHONE=C        
[11] LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] coda_0.19-2           forcats_0.3.0         stringr_1.3.1         dplyr_0.8.3           purrr_0.3.3          
 [6] readr_1.3.1           tidyr_0.8.2           tibble_2.1.3          ggplot2_3.1.0         tidyverse_1.2.1      
[11] igraph_1.2.4.1        statnet_2019.6        tsna_0.3.0            sna_2.4               ergm.count_3.3.0     
[16] statnet.common_4.3.0  tergm_3.6.1           networkDynamic_0.10.0 ergm_3.10.4           network_1.15        

{ergm}の使い方

github.com

  • アメリカ西海岸のネットワーク分析の研究者らが中心となって開発されているパッケージです。
  • 社会ネットワーク分析用の{sna}やネットワークオブジェクトを管理する{network}などが全部入りの{statnet}というパッケージスイートがあり、こちらを導入するのが便利です。
> library(statnet)
statnet: version 2019.6, created on 2019-06-13
Copyright (c) 2019, Mark S. Handcock, University of California -- Los Angeles
                    David R. Hunter, Penn State University
                    Carter T. Butts, University of California -- Irvine
                    Steven M. Goodreau, University of Washington
                    Pavel N. Krivitsky, University of Wollongong
                    Skye Bender-deMoll
                    Martina Morris, University of Washington
Based on "statnet" project software (statnet.org).
For license and citation information see statnet.org/attribution
or type citation("statnet").
  • 今回は、メディチ家の婚姻ネットワークデータ flomarriage を用います。
  • とりあえず plot してみます。
data(florentine)
plot(flomarriage)

f:id:meana0:20191219230711p:plain

  • ネットワーク統計量は summary() 関数を使うと計算できます。
    • 今回は、エッジの本数(edges)とトライアングル(triangle)に着目します。
> summary(flomarriage ~ edges + triangle)
   edges triangle 
      20        3 
  • これらのネットワーク統計量からなるモデルを構築し、ERGMのパラメータを推定してみます。
> ergm_model <- ergm(flomarriage ~ edges + triangle)
Starting maximum pseudolikelihood estimation (MPLE):
Evaluating the predictor and response matrix.
Maximizing the pseudolikelihood.
Finished MPLE.
Starting Monte Carlo maximum likelihood estimation (MCMLE):
Iteration 1 of at most 20:
Optimizing with step length 1.
The log-likelihood improved by 0.003736.
Step length converged once. Increasing MCMC sample size.
Iteration 2 of at most 20:
Optimizing with step length 1.
The log-likelihood improved by 0.0002205.
Step length converged twice. Stopping.
Finished MCMLE.
Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
This model was fit using MCMC.  To examine model diagnostics and check
for degeneracy, use the mcmc.diagnostics() function.
  • 結果を確認してみます。通常の lm と同じく、summary() 関数で推定結果の詳細を確認できます。
>summary(ergm_model)
==========================
Summary of model fit
==========================

Formula:   flomarriage ~ edges + triangle

Iterations:  2 out of 20 

Monte Carlo MLE Results:
         Estimate Std. Error MCMC % z value Pr(>|z|)    
edges     -1.6743     0.3488      0  -4.801   <1e-04 ***
triangle   0.1626     0.5894      0   0.276    0.783    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

     Null Deviance: 166.4  on 120  degrees of freedom
 Residual Deviance: 108.1  on 118  degrees of freedom
 
AIC: 112.1    BIC: 117.6    (Smaller is better.)

実装編

  • さて、続いてはRでできるだけゼロからERGMを実装することを試みます。

下準備

  • flomarriageは network オブジェクトなので、 igraph オブジェクトに変換します。
    • {intergraph} を用いると、この作業が非常に簡単になります。
library(igraph)
library(intergraph)
g <- intergraph::asIgraph(flomarriage)
plot(g)

f:id:meana0:20191219233118p:plain

ネットワーク統計量を計算する関数

  • 本当は行列計算でも書けるのですが…
#エッジの数
count_edges <- function(g){
  edges <- length(E(g))
  return(edges)
}

#トライアングルの数
count_triangle <- function(g){
  triangle <- motifs(g)[4]
  if(is.na(triangle)){
    triangle <- 0
  }
  return(triangle)
}
  • 先ほどの結果と一致します。
> count_edges(g)
[1] 20
> count_triangle(g)
[1] 3

推定のステップ

  • ERGMのパラメータ推定は以下のステップで行われます。

(1) 疑似尤度最大化(MPLE)による  \theta_0 の決定

(2) MCMCによるランダムグラフ生成

(3) (2)で得られたサンプルを用いた対数尤度比の最大化 ->  \hat{\theta}の推定

  • (2),(3)のステップはモンテカルロ尤度最大化(MC-MLE)と呼ばれています。
  • 少しトリッキーになりますが、MC-MLEの説明から始めます。

変化統計量

  • MC-MLE自体の説明に入る前に、まず、変化統計量について説明します。
  • 他のすべてのダイアドの状態  y_{ij}^c を一定として、ダイアド  ij の状態を変化させた時の、ネットワーク統計量の変化を 変化統計量(change statistics)と呼びます。
 \displaystyle
\Delta g(y)_{ij} = g(y^+)_{ij} - g(y^-)_{ij}

MCMCによるランダムグラフ生成

  • ある固定された  \theta_0 のもとで、変化統計量の比 [\pi] を計算します。
  •  \theta_0 を上手く決めれば、観測されたネットワークに近いようなランダムグラフが生成されていきます。
 \displaystyle
\pi = \frac{ \Pr(Y_{ij}\,changed \ |\ Y_{ij}^c = y_{ij}^c ) }{ \Pr(Y_{ij}\,not\,changed \ |\ Y_{ij}^c = y_{ij}^c ) } = \theta_0^T\ [\pm \Delta g(y)_{ij} ]
  •  \min(1,\pi) を遷移確率として、次のサンプルに進みます。こうしてランダムグラフのサンプルを生成していきます。
##ダイアドの状態をスイッチさせる
change_edge <- function(graph,i,j){
  if(graph[i,j] == 1){
    graph[i,j] <- 0
  }else{
    graph[i,j] <- 1
  }
  return(graph)
}

##ランダムにダイアドの状態をスイッチさせたグラフを返す
random_propose <- function(graph){
  sampled_pairs <- sample(1:length(V(graph)),2,replace = FALSE)
  i <- sampled_pairs[1]
  j <- sampled_pairs[2]
  graph <- change_edge(graph,i,j)
  return(graph)
}

MC-MLE

metropolis <- function(graph, theta0_edges, theta0_triangle){
  
  ##現在のネットワーク
  current_graph <- graph
  current_edges <- count_edges(current_graph)
  current_triangle <- count_triangle(current_graph)
  
  ##提案されたネットワーク
  proposed_graph <- random_propose(graph)
  proposed_edges <- count_edges(proposed_graph)
  proposed_triangle <- count_triangle(proposed_graph)
  
  ##変化統計量の計算
  changed_edges <- proposed_edges - current_edges
  changed_triangle <- proposed_triangle - current_triangle
  
  ##piを計算し、提案されたグラフを採択するかどうか決める
  change_ratio <- exp(theta0_edges*changed_edges + theta0_triangle*changed_triangle)
  move_probability <- min(1,change_ratio)
  if(move_probability >= runif(1)){
    return(proposed_graph)
  }else{
    return(current_graph)
  }
}
MCMC <- function(graph, burn_in = 10000, samplesize = 50000, theta0_edges = 0.1, theta0_triangle = 0.1){
  graph <- graph.empty(length(V(graph)),directed=FALSE)
  edges_c <- count_edges(graph)
  triangle_c <- count_triangle(graph)
  par(mfrow = c(2,2))
  for(i in 1:samplesize){
    if(i %% 1000 == 0){
      progress_text <- paste0(i,"(",round(100*(i/samplesize),3),"%) done")
      print(progress_text)
      plot(edges_c,type="l")
      hist(edges_c)
      plot(triangle_c,type="l")
      hist(triangle_c)
    }
    after_graph <- metropolis(graph,theta0_edges,theta0_triangle)
    edges_c <- c(edges_c,count_edges(after_graph))
    triangle_c <- c(triangle_c,count_triangle(after_graph))
    graph <- after_graph
  }
  return(list(edges = edges_c[(burn_in+1):samplesize], triangle = triangle_c[(burn_in+1):samplesize],obs_graph = graph))
}

MPLE

  • さて、ここまでMCMCを用いたランダムグラフの生成は  \theta _ 0 を前提としていましたが、この値自体はどのように決めればいいのでしょうか。
  • 実は、 \theta _ 0 \theta の真値に近くないと、正しい分布からランダムグラフを生成してくれないという問題があるため、できるだけ近い値を  \theta _ 0 としたいです。
  • そこで、疑似尤度の最大化(Maximum Pseudo-Likelihood Estimation, MPLE)で推定されたパラメータを  \theta _ 0 とします。
    • これは、ダイアド同士の統計的連関を考慮しないモデル (dyad-independence) と同値であり、構造的特徴のないモデルではMCMCを行わず、MPLEのみが使われます。
  • これは非常にシンプルで、観測されたネットワークのダイアドにおけるエッジの有無を、ネットワーク統計量の変化統計量で回帰するロジスティック回帰を行うことで達成されます。
#変化統計量とエッジの存在からなるデータフレームを生成
generate_changed_data <- function(graph){
  changed_data <- data.frame()
  for(i in 1:length(V(graph))){
    for(j in 1:length(V(graph))){
      if(i >= j){
        next
      }
      existing_if <- graph[i,j]
      changed_graph <- change_edge(graph,i,j)
      changed_triangle <- abs(count_triangle(changed_graph) - count_triangle(graph))
      row <- data.frame(i = i, j = j, existing_if = existing_if, changed_triangle = changed_triangle)
      changed_data <- rbind(changed_data,row)
    }
  }
  return(changed_data)
}

#MPLEの実装
MPLE <- function(graph){
  data_for_mple <- generate_changed_data(graph)
  reg <- glm(data = data_for_mple, existing_if ~ changed_triangle, family = "binomial")
  return(list(edges = reg$coefficients["(Intercept)"], triangle = reg$coefficients["changed_triangle"]))
}
  • MPLE -> MC-MLEまでのプロセス
## SAMPLING
triangle_model <- function(graph, burn_in = 10000, samplesize = 50000){
  message("Starting MPLE fitting...")
  theta0_list <- MPLE(graph)
  print(theta0_list)
  message("MPLE fitting done. Now starting MC-MLE fitting...")
  fit <- MCMC(graph,
              burn_in = burn_in,
              samplesize = samplesize,
              theta0_edges = theta0_list$edges,
              theta0_triangle = theta0_list$triangle)
  fit$theta0_edges <- theta0_list$edges
  fit$theta0_triangle <- theta0_list$triangle
  return(fit)
}
  • 実際にMCMCからネットワーク統計量のサンプルを得ます。
fit <- triangle_model(g, burn_in = 10000, samplesize = 50000)

対数尤度比の最大化

  • さて、MCMCからネットワーク統計量のサンプルが得られました。次はこれを用いて、実際にパラメータを推定していきます。
  • 尤度それ自体は正規化定数を計算しなければ最大化できないので、代わりに、ある固定された \theta_0 を用いた  \ell(\theta)-\ell(\theta_0) =対数尤度比の最大化を目指します。
  •  \ell(\theta)-\ell(\theta_0) は、以下のように近似できます。
    •  m : MCMCで得られたランダムグラフで、近似に用いる数。
    •  y_i : MCMCで得られたランダムグラフで、近似に用いる  i 番目のサンプル。
 \displaystyle
\ell(\theta)-\ell(\theta_0) \approx (\theta - \theta_0)^T g(y) - \log[ \frac{1}{m}\sum_{i=1}^m \exp\{ (\theta-\theta_0)^T g(y_i) \} ]
#対数尤度比を計算する関数
loglikelihood_ratio <- function(theta,fit,n){
  sample_length <- length(fit$edges)
  sample_index <- (sample_length-n):sample_length
  sampled_edge_count <- fit$edges[sample_index]
  sampled_triangle_count <- fit$triangle[sample_index]
  theta0 <- c(fit$theta0_edges,fit$theta0_triangle)
  g_obs <- c(count_edges(fit$obs_graph),count_triangle(fit$obs_graph))
  first = (theta - theta0) %*% g_obs
  exp_vec <- c()
  for(i in 1:n){
    sam <- exp((theta - theta0) %*% c(sampled_edge_count[i],sampled_triangle_count[i]))
    exp_vec <- c(exp_vec,sam)
  }
  second <- log(sum(exp_vec)/length(exp_vec))
  loglikelihood_ratio <- first - second
  return(loglikelihood_ratio)
}
  • optim() 関数を使って、MCMCで得られたグラフのネットワーク特徴を与えて、対数尤度比を最大化してみます。
##対数尤度比の最大化
my_estimates <- optim(par = c(0,0), #初期値
                      fn = loglikelihood_ratio, #目的関数
                      control = list(fnscale=-1), #最大化するため
                      fit = fit,
                      n = 1000) #MCMCサンプルサイズ
  • さて、いよいよオレオレERGMの推定結果を確認します。
  • {ergm}で求めたパラメータと比べてみます。
> ergm_model$coef #推定結果 ({ergm})
    edges  triangle 
-1.674332  0.162630 
> my_estimates$par #推定結果 (オレオレergm)
[1] -1.8267176  0.2983575
  • それなりに近い推定値にはなっているものの、ややズレてしまっています。

問題:MCMCが上手くランダムグラフを生成していない

  • 当初は最適化関数の扱いに問題があるのかと思い、初期値などを色々試してみるも無駄に終わってしまいました。
  • 調べてみると、MCMCが上手くランダムグラフを生成していまいことがわかりました。
  • まず、{ergm}を用いて生成したmcmcsampleは以下のようなtraceplotを呈します。
    • わかりやすくするために、観測されたグラフのネットワーク統計量が0になるようにセンタリングしています。
> mcmc.diagnostics(ergm_model, esteq = FALSE, center = TRUE)

f:id:meana0:20191219235557p:plain

  • 一方、自分で実装したmcmcsampleを見てみます。
  • ergm::mcmc.diagnostics() と似たプロットを行う関数を作ります。
##オレオレMCMCのサンプル生成結果をplot
diagnosis <- function(fit){
  edges_c <- fit$edges
  triangle_c <- fit$triangle
  obs_edges <- count_edges(fit$obs_graph)
  obs_triangle <- count_triangle(fit$obs_graph)
  edges_centered <- edges_c - obs_edges
  triangle_centered <- triangle_c - obs_triangle
  par(mfrow = c(2,2))
  plot(edges_centered,type="l",)
  hist(edges_centered)
  plot(triangle_centered,type="l")
  hist(triangle_centered)
}
  • 結果がこちらです。
> diagnosis(fit)

f:id:meana0:20191220000006p:plain

  • 観測値に比べるとマイナス側に偏った分布になってしまっています。

さいごに

  • 今回は時間が足りず、これ以上のことはできませんでした。
  • 尻切れトンボにようになってしまいますが、なぜ上手く推定できないのかも含めて、この課題に引き続き取り組んでいこうと思います。

参考文献

{parallelDist}が鬼速い

R Advent Calendar 2019 15日目の記事です。

qiita.com

タイトルだけで内容の99%は伝わっていると思うのですが、{parallelDist} ライブラリが並列でベクトル間の距離計算をしてくれて、鬼のように速いという話をします。

github.com

締め切り前なのに距離計算が間に合わない!という時に重宝します。

Dynamic Time Warpingとは

  • 今回はDynamic Time Warpingという距離を例にparallelDistの威力を紹介します。
  • 時系列シーケンス間の距離。時系列クラスタリングによく用いられる。
    • 平行移動や長さ・周期の違いに対してロバストな類似度を計算できる。
    • 先人による解説:

sinhrks.hatenablog.com

  • 一方で計算量がとても大きい。
    • 2つのシーケンスを総当りで計算するので、長さnの系列と長さmの系列のDTW距離を計算しようとするとO(nm)かかる。

セッション情報

  • WSL上で動かしています。
sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

パッケージインストール

install.packages("TSclust")
install.packages("dtw")
install.packages("parallelDist")

ライブラリ読み込み

library(magrittr)
library(dplyr)
library(tidyr)
library(purrr)

ランダムウォークによる系列の生成

random_walk <- function(x,t){
  return(cumsum(rnorm(t)))
}

set.seed(1024)
series_data <- map_dfc(1:100,random_walk,t=100) #この記法は重宝しますね
plot_data <- series_data %>% 
  mutate(period = row_number()) %>%
  gather(key = series, value = value,-period)
g <- ggplot(plot_data,aes(x = period, y = value, group = series))
g <- g + geom_line(stat="identity")
g

f:id:meana0:20191215014654p:plain

{TSclust}を使う

  • DTWを計算する時に一番多く使われているのはこのパッケージかと思います。
    • 先人による解説:

sinhrks.hatenablog.com

tsclust_time <- system.time({
  tsclust_dist <- TSclust::diss(series_data,"DTWARP")
})
tsclust_time["user.self"]
user.self 
   19.609 

{dtw}を使う

dtw_time <- system.time({
  dtw_dist <- dtw::dtwDist(t(series_data)) #dtwの場合は時系列を行として配置
})
dtw_time["user.self"]
user.self 
   22.141 

parallelDistを使う

  • とにかく鬼速く距離計算してくれるライブラリ
  • 何故速いのか?
    • RcppParallel によるマルチスレッド計算
    • RcppArmadillo による行列計算の最適化
    • CRTPによる仮想関数の回避(よくわからない…)
  • GitHub : https://github.com/alexeckert/parallelDist
pardist_time <- system.time({
  pardist_dist <- parallelDist::parDist(t(series_data),
                                        method = "dtw",
                                        step.pattern="symmetric2") #dtwのデフォルトと合わせる
})
pardist_time["user.self"]
user.self 
    2.922 
  • お、鬼速い…((((;゚Д゚))))
  • 約7倍の速度向上。

結果は同じ?

  • どのライブラリでも結果は同じなのか恐る恐る確認してみる。
as.matrix(tsclust_dist)[1:10]
as.matrix(dtw_dist)[1:10]
as.matrix(pardist_dist)[1:10]
 [1]    0.0000  399.4580  207.6915  150.0798  733.7981 2201.4722  345.4040
 [8] 1569.0761 1201.3074  415.9536
 [1]    0.0000  399.4580  207.6915  150.0798  733.7981 2201.4722  345.4040
 [8] 1569.0761 1201.3074  415.9536
 [1]    0.0000  399.4580  207.6915  150.0798  733.7981 2201.4722  345.4040
 [8] 1569.0761 1201.3074  415.9536

階層的クラスタリング

  • せっかくなので階層的クラスタリングも高速な {fastcluster}パッケージを使う。
hc <- fastcluster::hclust(pardist_dist)
plot(hc)
trees <- cutree(hc,3)

f:id:meana0:20191215015024p:plain

クラスタリング結果の可視化

cluster_data <- tibble(series = names(trees),
                       cluster = as.factor(trees))
plot_data <- series_data %>% 
  mutate(period = row_number()) %>%
  gather(key = series, value = value,-period) %>%
  left_join(cluster_data,by="series")
g <- ggplot(plot_data,aes(x = period, y = value, group = series, col = cluster))
g <- g + geom_line(stat="identity")
g

f:id:meana0:20191215015054p:plain

  • 時点数Tと系列数Nと、{dtw}比の計算時間の関係を検証したかったのですが、時間とマシンスペックが足りず…。
    • 暇な時にやろうかと思います。

補足

  • {TSclust}が依存している{rgl}を入れる時にunable to open X11 displayと怒られたが、Sys.setenv("RGL_USE_NULL"=TRUE)で解決した。

朝型人間に擬態するためにやっていること

親譲りでない夜型人間で小供の時から損ばかりしている。
―― 夏目漱石『夜型人間ちゃん』

夜型人間 Advent Calendar 2019 1日目の記事です。

2年ほど前にフルタイムワーカーを始めてからというもの、自分なりに色々と試行錯誤しながら、朝型人間に「擬態」する術を身につけてきました。

自分は根っからの「夜型人間」だと自認しており、修士課程の時に概日リズムをぶっ壊してからも、いまだに朝は自分の時間ではないと思っています。作業が捗り、良いアイデアが降ってくるのは大抵22時以降です。

wired.jp

最近、「夜型か朝型かはゲノムレベルで規定されており、この違いを努力でカバーすることは困難である」という趣旨の記事が話題になりました。この知見に関する妥当性は、専門外のため保留させていただきますが、「朝型人間と夜型人間の共存」についての議論のきっかけになることを願っています。

しかし、規範や制度が急速に変わることは稀です。生活リズムというテーマに関しても、明日いきなり社会の寛容度が向上することは、おそらく見込めないでしょう。さしあたり個人の幸福にとって重要なのは、「どうすれば朝型の社会をやり過ごすことができるか」です。この記事では、「環境構築」を中心に、自分が朝型人間に「擬態」するために導入したガジェットや、普段実践している生活習慣を紹介します。

といっても、人類はみな朝型人間か夜型人間のどちらかというわけではなく、両極の間のスペクトラムだとは思いますし、僕以上に深刻な夜型の方もいらっしゃることでしょう。場合によっては、医療機関への相談のほうが近道かもしれません。

それに、僕自身がこの記事で紹介することを実践しながらも少なくとも月に1,2回「やらかして」いるので、話半分に読んでいただけると幸いです。専門的な見地からみて間違っていることがあればご指摘ください。

この記事を通して、「いつも遅刻してくるあいつ」も実は頑張っているかもしれないですよ、ということが伝えられたらよいと思っています。

環境構築

朝型人間への擬態をする上での僕の基本戦略は「環境の改変」です。夜型人間なら誰しも、朝に起きたいと思っても起きられない朝をいくつも経験したことでしょう。眠る前の自分はあまりに楽天家で、寝覚めの瞬間の自分は、私たちが想像するよりも遥かに狡猾です。自分の意思ほど信頼できないものはありません。

環境の中には、寝具や耳栓のように入眠確率を上げる類のものと、目覚ましのように起床確率を上げる類のものがあると思っています。

寝具

博士課程に入りたての頃、祈るようにしてマニフレックスマットレスと枕を購入しました。当時としてはかなり思い切った判断でした。東急ハンズの寝具売り場を足繁く通っていた記憶があります。長い間すのこベッドに薄い敷ふとんを敷いて寝ていたのですが、睡眠環境は見違えるように変化しました。マットレスも枕も肉厚なため、体重85kg超えの自分でも底に身体が当たることはなく、非常に快適です。他のメーカーも試しましたが、最終的には厚みで選んだ形です。この投資は非常に効果があり、「なんか固くて寝にくいな…」という入眠時の思考のノイズを減らせた気がします。

最近は、パラマウントベッドが販売しているActive Sleepというベッド+マットレスが気になっていますが、いかんせん高価なので手が出ませんね…。

activesleep.jp

耳栓

耳栓も、入眠確率を上げるための重要なツールです。僕は人の声や環境音が気になって眠れなくなるタイプなので、耳栓はほぼ必須です。AmazonでMOLDEXの耳栓が200ペアまとめ買いできるので、これを定期的に購入しています。

耳栓をすると目覚ましの音が聞こえないのではないかと思われるかもしれませんが、自分の場合は寝返りのせいか、朝起きたときには外れていることが多く、そこまで起床の妨げにはなっていません。もちろん個人差あり、起床確率が下がる方向に作用する方もいらっしゃると思うので、注意してください。

目覚まし

続いては起床確率を上げるための目覚ましシステムです。普段使いしているのは5つです。音・光・振動の感覚複合的なアプローチを採用しています。

  • デジタル目覚まし時計

    • 普通のデジタル時計です。もう10年くらい使っていますね。
  • スマートフォンのアラーム

    • LGV30+デフォルトのアラームアプリです。
  • スマホのSleepCycle

    • マイクで睡眠中の寝返りの音などを拾って、睡眠の浅い時に起こしてくれるアプリです。
    • これを導入してからは本当に高確率で起きられるようになりました。
    • 過去の睡眠の傾向を分析してくれるアナリティクス機能(有料)は使用していません。

Sleep Cycle: スマートアラーム目覚まし時計

Sleep Cycle: スマートアラーム目覚まし時計

  • Sleep Cycle AB
  • ヘルスケア/フィットネス
  • 無料
apps.apple.com

  • 簡易的な光目覚まし
    • スイッチを押してくれるIoTデバイスSwitchBotを利用して、朝8時に部屋の照明がつくように設定しています。
    • 室内照明だけで十分効果的を感じています。
      • intiのようなとても強力な光目覚ましもありますが、試したことはないです。
    • 部屋の照明を使うという観点では、Philips HueやNature Remoでも思います。

Nature スマートリモコン Nature Remo mini Remo-2W1

Nature スマートリモコン Nature Remo mini Remo-2W1

  • HuaweiBand
    • スマートウォッチの一種です、アラーム機能が優秀です。
    • 寝返りを記録して、睡眠の浅い時にバイブレーションで起こしてくれます。
    • Apple WatchやFitBitでも同様の機能があるかと思います。

  • その他
    • 早朝にフライトがある時など、どうしても起きなければいけない場合は、ここにもう一つタブレットのアラームアプリを足しています。
    • 一時期はBose Sleepbudsを愛用していたのですが、バッテリーの不具合が出て以降は使うのをやめてしまいました。個人的には最高のガジェットだったのですが…。

生活習慣

環境の改変ではないですが、普段の生活習慣については、該日リズムを壊した時に医者から勧められたことを実践し続けています。つくづく「入眠・起床の確率を上げるための戦いは、日中から既に始まっている」ということに尽きます。

  • 19時以降にカフェインを摂取しない。お茶も飲まない。
    • カフェインの半減期は4-6時間のようです。24時過ぎに就寝することを目標としているので、そこから逆算しています。

e-heartclinic.com

  • 23時を過ぎたらPCを見ない。
    • 「アイデアを思いつくのは決まって深夜問題」がありますが、最近では楽しいアイデアを思いたら、どこかにメモしておいて、とりあえず眠るようになりました。
    • それでも夜にPCを触らなければいけない時のために、f.luxなどを常にオンにしておき、モニターの輝度を抑えています。

justgetflux.com

  • 運動する

    • 週2回にジム、毎日1万歩歩くことを目標としています。
      • 普通に通勤するだけだと6000歩ほどですが、帰宅時にたまに1駅前で降りたりしています。
        • HuaweiBandは1万歩を超えると知らせてくれます。
    • リングフィットアドベンチャーやフィットボクシングでもよいですね。神ゲーの登場以降ほとんどジムに行っていません。
    • 一方で、遅い時間の運動はむしろ起床確率を下げる気もします。
  • 休日に寝溜めをしない。

    • こればっかりは環境でどうこうできる問題ではないような気もしますが、休日の朝に楽しい予定を入れておくと、リズムが保ちやすいのではないかと思っています。何より、寝すぎると偏頭痛がつらいです。
      • 自分の場合は最近マイブームの水族館や野鳥撮影が大きなインセンティブになっています。
    • 大型連休は、睡眠リズムが壊れやすいので天敵です…。

以上、根っからの夜型人間の僕が朝型人間に擬態するためにやっていることを紹介してきました。 前述したとおり、これだけやってもうまく眠れなかったり朝起きられないこともありますが、「なんとかやれる」水準まで持って行けている気がしています。毎日、夜に眠気が来るとホッとします。

皆さまの朝型人間への擬態も、うまくいくようにお祈りしています。

"Thanks God I'm Falling asleep!"