この方、 twitter を見てたら 今日は 食べ歩きをしているようだが、その分ちゃんと運動しているのか 余計なお世話ながら気になったので データから確認したい。

何のモチベーションかこういったデータ が公開されているので、運動状況として歩数データをみればよさそうだ。が、直近数日は更新がないため、過去のデータをモデル化して直近のトレンド / 休日の歩数の傾向を推測したい。

必要パッケージのインストール

library(devtools)
install_github('hoxo_m/pforeach')
install_github('sinhrks/ggfortify')

パッケージのロード

library(rstan)
library(Nippon)
library(pforeach)
library(ggfortify)

前処理

天下り的ではあるが、以下のデータを追加する

# ore.csv は 以下 URL から入手
# https://github.com/dichika/mydata
ore <- read.csv('ore.csv')
ore$time <- as.Date(ore$time)

is.holiday <- function(date) {
  # 祝日以外に休日として扱う日
  calendar <- as.Date(c('2015-01-02'))
  (weekdays(date) %in% c("土曜日", "日曜日")) | Nippon::is.jholiday(date) | (date %in% calendar)
}
ore$holiday <- is.holiday(ore$time)

get_weeknum <- function(date) {
  as.numeric(format(date, "%U"))
}
ore$weeknum <- get_weeknum(ore$time)

データの確認

モデル作成のため、周期性などをざっと確認。

walk <- as.ts(ore$data)
# 時系列プロット
autoplot(walk, geom = 'bar', fill = 'blue4')

# 自己相関はなさげ
gglagplot(walk, lag = 6, ncol = 3)

autoplot(acf(walk, plot = FALSE))

# 周期性 7 ごとでプロットして曜日の周期性を確認。同じ曜日でも結構ばらつく
ggfreqplot(walk, freq = 7, geom = 'point', ncol = 7)

# 曜日よりは休日かどうかでまとめたほうが傾向ありそう
ggplot(ore) + geom_point(aes(x = time, y = data)) + facet_wrap(~holiday) 

# (日時だとばらつきすぎるので) 週ごとにトレンドとればよさそう
ggplot(ore) + geom_boxplot(aes(x = as.factor(weeknum), y = data))

モデル作成

各変数の歩数への影響が見たいので状態空間モデルを使う。上の結果から、以下の仮定でモデルつくってみる。

状態空間時系列分析入門的にいうと 干渉変数のある週次トレンド ore モデル。案ができたので Stan でコード書く。

standata <- within(list(), {
  y <- as.numeric(ore$data)
  w <- as.numeric(ore$weeknum)
  h <- as.numeric(ore$holiday)
  n <- length(y)
  wn <- length(unique(ore$weeknum))
})

code = "data {
  int<lower=1> n;
  # 週番号の要素数
  int<lower=1> wn;
  vector<lower=0>[n] y;

  # 週番号
  int w[n];
  # 休日かどうかのフラグ
  int h[n];
}
parameters {
  # 週番号ごとのトレンド
  vector<lower=0>[wn] mu;
  # 休日の場合の干渉
  real holiday;
  # 平日、休日の撹乱項
  vector<lower=0>[2] sigma_irreg;
}
transformed parameters {
  vector[n] yhat;
  for(t in 1:n) {
    yhat[t] <- mu[w[t] + 1] + h[t] * holiday;
  }
}
model {
  for(t in 1:n)
    y[t] ~ normal(yhat[t], sigma_irreg[h[t] + 1]);
}"

stan_fit <- stan(model_code = code, chains = 0)
## 
## TRANSLATING MODEL 'code' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'code' NOW.
sflist <- pforeach(i = 1:3)({
  stan(fit = stan_fit, data = standata, chains = 1, seed = i)
})
fit <- sflist2stanfit(sflist)

# 収束したかチェック
is.converged <- function(stanfit) {
  summarized <- summary(stanfit)  
  all(summarized$summary[, 'Rhat'] < 1.1)
}
stopifnot(is.converged(fit))

# 結果の取り出し
(mu <- get_posterior_mean(fit, par = 'mu')[, 'mean-all chains'])
##     mu[1]     mu[2]     mu[3]     mu[4]     mu[5]     mu[6]     mu[7] 
##  5173.077  6430.400  7059.242  6761.704  6778.334  7266.128 12298.069 
##     mu[8]     mu[9]    mu[10] 
##  7141.936  7015.875  6306.878
yhat <- get_posterior_mean(fit, par = 'yhat')[, 'mean-all chains']
(holiday <- get_posterior_mean(fit, par = 'holiday')[, 'mean-all chains'])
## [1] -2088.726

結果のプロット

モデルでは休日の歩数が 平日と比べ 2000 歩ほど少なく出ている。が、逆の動きをしている場合も結構あるので 上のモデルでは説明が厳しそうだ。

p <- autoplot(as.ts(ore$data), colour = 'blue4')
autoplot(as.ts(yhat), p = p, colour = 'blue', linetype = 'dashed')

今後の展望

改善のためにはその日の活動をより直接的に示す説明変数が必要そうだ。room.csv のデータを追加して ストーキング モデルを改善したい。