時系列データ

[2017-10-01 追記] macOS High SierraにしてからタイムゾーンがJSTでなく警告が出てGMTになってしまう気がします。.bash_profileexport TZ=":/etc/localtime" と書いておけばよさそうです(How setting the TZ environment variable avoids thousands of system calls 参照)。Rに限っては Sys.setenv(TZ=":/etc/localtime") あるいは具体的に Sys.setenv(TZ="Asia/Tokyo") などと打ち込むか ~/.Rprofile に書いておくという手もありそうです。 ←このあたりはRのバージョンアップで直った?

Rでの日時の表し方

Rで現在の日時を調べるにはつぎのようにする。

Sys.time()
[1] "2011-03-20 12:34:56 JST"

この形式がISOで決められた時間の表し方である。最後のJSTは日本標準時を表すが,自明なら "2011-03-20 12:34:56" のように省略してもよい。さらに秒が不要なら "2011-03-20 12:34" でよいし,日付だけでよければ "2011-03-20" でよい(これは "2011-03-20 00:00" と同じことである)。

これらは文字列のように見えるが,内部表現は1970年元旦からの秒数で,POSIXct型という衣を着ている(「ct」はcalendar time)。衣(クラス)を剥がしてみれば,整数値が現れる。単に数値に変換しても同じ結果になる。

unclass(Sys.time())
[1] 1300592096
as.numeric(Sys.time())
[1] 1300592096

"2011-03-20 12:34:56" のような文字列からPOSIXct型の値に変換するには

t = as.POSIXct("2011-03-20 12:34:56")

のようにすればよい。"2011-03-20" の代わりに "2011-3-20" でも "2011/3/20" でもよい。別の方法として次のようにもできる。

t = ISOdatetime(2011, 3, 20, 12, 34, 56)

POSIXct型は内部では秒単位の整数であるため,秒を加算することができる。例えば時刻 t を起点として10分ごと1時間の数列を作るには次のようにする(最後の 600"10 min" でもよい)。

seq(t, t+3600, 600)

ちなみに,日本の環境では日時は日本標準時(JST)で扱われる(ただし下で述べるようにJDTという例外がある)。JSTは協定世界時(UTC)から9時間(32400秒)ずれている:

as.numeric(as.POSIXct("1970-01-01 00:00:00"))
[1] -32400
as.numeric(as.POSIXct("1970-01-01 00:00:00", tz="Asia/Tokyo"))
[1] -32400
as.numeric(as.POSIXct("1970-01-01 00:00:00", tz="Japan"))
[1] -32400
as.numeric(as.POSIXct("1970-01-01 00:00:00", tz="UTC"))
[1] 0

POSIXct型以外にPOSIXlt型がある(「lt」はlocal time)。これは年月日時分秒をリストとして並べたものである。

日時を表す文字列をPOSIXct型に変換する際に,まず strptime() によりPOSIXlt型に変換される。ここでシステムコールが発生し,遅い。より速く変換するには fasttime パッケージを使う(datatable 参照)。

日単位ならDate型を使う。これは内部的には1970-01-01からの日数である。as.numeric(as.Date("1970-01-01")) の値は0になる。

簡単な時系列グラフ

ここまでのまとめとして,簡単な時系列データのグラフを描いてみよう。文科省が都道府県別環境放射能水準調査結果[リンク切れ]で公開している福島第一原発の事故に伴う3月15日17-18時から1時間ごとの東京都新宿区で計測した放射線の量は次の通りである。

x = c(0.0941, 0.2, 0.361, 0.123, 0.0888, 0.0657, 0.0556)

このグラフは次のようにして描ける。

t0 = as.POSIXct("2011-03-15 17:00")
t = seq(t0, by=3600, along.with=x)
plot(t, x, type="o", pch=16, xlab="時刻", ylab="γ線(μSv/h)")

もっとデータを増やしてみよう。

x = c(0.0941, 0.2, 0.361, 0.123, 0.0888, 0.0657,
      0.0556, 0.0538, 0.0547, 0.0672, 0.101, 0.141,
      0.143, 0.142, 0.104, 0.0891, 0.069, 0.058,
      0.057, 0.056, 0.055, 0.054, 0.054, 0.054)

これくらいになると,時間軸の目盛りが思うようにならない。目盛りを省略(xaxt="n")して,別に目盛りだけを描くほうがよい。

plot(t, x, type="o", pch=16, xlab="時刻", ylab="γ線(μSv/h)", xaxt="n")
axis.POSIXct(1, at=t, format="%d日%H時")

目盛りを4時間ごとにしてみよう。

plot(t, x, type="o", pch=16, xlab="時刻", ylab="γ線(μSv/h)", xaxt="n")
r = as.POSIXct(round(range(t), "hours"))
axis.POSIXct(1, at=seq(r[1],r[2],by="4 hours"), format="%d日%H時")

このように,POSIXct型では sec(s)min(s)hour(s)day(s) という単位が使える。また,format の形式はRで strptime のヘルプを見れば出ている。よく使うものだけを挙げておく:

%Y年(4桁)
%y年(2桁)
%m月(01-12)
%b月( 1-12←日本語環境では)
%d日(01-31)
%e日( 1-31)
%H時(00-23)
%M分(00-59)
%S秒(00-61)

タイムゾーンについて

R 3.5.0/3.5.1 (Windows/Mac) で as.POSIXct("1948-05-02") がエラーになる:

as.POSIXct("1948-05-01")
[1] "1948-05-01 JST"
as.POSIXct("1948-05-02")
 as.POSIXlt.character(x, tz, ...) でエラー: 
   文字列は標準的な曖昧さのない書式にはなっていません 
as.POSIXct("1948-05-03")
[1] "1948-05-03 JDT"

原因はどうやら日本でも1948年〜1951年にあった夏時間(JDT)のためらしい。"1948-05-02""1948-05-02 00:00:00" と解釈されて,夏時間を考慮すれば存在しない時刻になるのだろう。ただ,同じコードでもLinux上ではエラーにならない:

as.POSIXct("1948-05-02")
[1] "1948-05-01 23:00:00 JST"

as.POSIXct() の第2引数は tz つまりタイムゾーンである。何も指定しなければシステムのタイムゾーンになる。日本語環境では Asia/Tokyo または同じことだが Japan になる。

as.POSIXct("1948-05-03", "Asia/Tokyo")
[1] "1948-05-03 JDT"
as.POSIXct("1948-05-03", "UTC")
[1] "1948-05-03 UTC"

"Japan" と書けば /usr/share/zoneinfo/Japan にある情報が使われ,Asia/Tokyo と書けば /usr/share/zoneinfo/Asia/Tokyo にある情報が使われるが,これらは同じものである。タイムゾーンを指定しなければ /etc/localtime の情報が使われるが,日本の環境ではこれは /usr/share/zoneinfo/Asia/Tokyo へのシンボリックリンクになっている(Macでは /var/db/timezone/zoneinfo/Asia/Tokyo へのシンボリックリンクだが,/usr/share/zoneinfo が /var/db/timezone/zoneinfo へのシンボリックリンクになっているので,結局同じことである)。これらのデータの内容は zdump コマンドで見ることができる:

zdump -v /usr/share/zoneinfo/Asia/Tokyo # または zdump -v Asia/Tokyo

この Asia/Tokyo にあたる部分にどのような選択肢があるかは,Rのコマンド OlsonNames() でも調べられる。

as.POSIXct() は,日本の環境では(あるいは明示的に Asia/Tokyo などと書けば)文字列を日本標準時(JST)として扱うが,1948年〜1951年の夏だけは,夏時間(JDT)として扱うので,1時間ずれる。春は進み(spring forward),秋は戻る(fall back)。春には存在しない時刻があり,エラーになるか,おかしな結果を返す。秋には同じ時刻が2回あり,結果はあいまいになる。このため,夏時間を採用している国では,時間の扱いに注意しなければならない。

気象庁の過去データの日時はJSTで統一されているそうである。この場合,タイムゾーンとして Etc/GMT-9 を指定すればよい:

as.POSIXct("1948-05-02 00:01", tz="Etc/GMT-9")
[1] "1948-05-02 00:01:00 +09"

日本以外の例を挙げると,米国東海岸の春(EST/EDT)は,Mac上のR 3.5.1では次のようになる:

as.POSIXct("2018-03-11 01:59", "US/Eastern")
[1] "2018-03-11 01:59:00 EST"
as.POSIXct("2018-03-11 02:00", "US/Eastern") // 存在しない
[1] "2018-03-11 EST"
as.POSIXct("2018-03-11 02:59", "US/Eastern") // 存在しない
[1] "2018-03-11 EST"
as.POSIXct("2018-03-11 03:00", "US/Eastern")
[1] "2018-03-11 03:00:00 EDT"

Linux上のR 3.5.1では次のようになる

as.POSIXct("2018-03-11 01:59", "US/Eastern")
[1] "2018-03-11 01:59:00 EST"
as.POSIXct("2018-03-11 02:00", "US/Eastern") // 存在しない
[1] "2018-03-11 01:00:00 EST"
as.POSIXct("2018-03-11 02:59", "US/Eastern") // 存在しない
[1] "2018-03-11 01:59:00 EST"
as.POSIXct("2018-03-11 03:00", "US/Eastern")
[1] "2018-03-11 03:00:00 EDT"

lubridate パッケージの parse_date_time("1948-05-02", "Ymd", "Asia/Tokyo") やその高速版 parse_date_time2("1948-05-02", "Ymd", "Asia/Tokyo") も同様である(parse_date_time2() はLinuxではエラーにならない)。fasttime パッケージの fastPOSIXct() は1970年以前は扱えないようである。

readrparse_datetime() は次のようにおかしくなる:

parse_datetime("1948-05-01", locale=locale(tz="Asia/Tokyo"))
[1] "1948-05-01 JST"
parse_datetime("1948-05-02", locale=locale(tz="Asia/Tokyo"))
[1] "1970-01-01 08:59:59 JST"
parse_datetime("1948-05-03", locale=locale(tz="Asia/Tokyo"))
[1] "1948-05-03 JDT"

[2018年8月追記] 2020年の東京オリンピックで+2時間の夏時間を採用して,7時開始を5時開始相当にしようという話があるらしい。夏の5時は米国西海岸では何時だろうか? as.POSIXlt() も援用して

as.POSIXlt(as.POSIXct("2020-07-24 05:00", tz="Asia/Tokyo"), tz="US/Pacific")
[1] "2020-07-23 13:00:00 PDT"

とすればよい(参考)。as.POSIXct() で統一したければ

as.POSIXct(as.numeric(as.POSIXct("2020-07-24 05:00", tz="Asia/Tokyo")), origin="1970-01-01 00:00", tz="US/Pacific")
[1] "2020-07-23 13:00:00 PDT"

とすればよいのだろう。

ちなみに,Excelにはこのような時間を扱う関数がないらしい(『Excelにタイムゾーンという概念がなく、サマータイムに対応しないということは』『これから未来永劫、時限サマータイムの2年間についての複雑な特別処理を仕込まないと、正しい集計ができなくなります』)。実際の処理については,例えば Automatically Converting to GMT 参照。

時系列オブジェクト

日時なら as.POSIXct(),日だけなら as.Date() で扱うのが基本であるが,月ごとのデータや曜日に意味がある日ごとのデータは,ts() で時系列オブジェクトに変換すると便利なことがある。例えばグラフの例:白血病による死亡率の推移のデータの場合,

leu = read.csv("https://oku.edu.mie-u.ac.jp/~okumura/stat/data/leukemia2.csv")
x = ts(leu$rate, start=c(2001,1), frequency=12)
x
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2001 5.4 5.2 5.4 5.8 5.6 5.8 5.5 5.7 5.7 5.3 5.7 5.3
2002 5.5 5.4 5.8 5.6 5.3 5.3 5.4 5.6 5.6 5.7 6.1 5.2
2003 5.8 5.8 5.5 5.2 5.4 5.6 5.4 5.0 6.0 6.1 5.5 5.5
……
2016 7.0 7.0 7.2 7.4 6.8 7.4 6.8 7.0 6.9 7.3 7.0 6.8
class(x)
[1] "ts"
head(time(x))
[1] 2001.000 2001.083 2001.167 2001.250 2001.333 2001.417

のように「ts」というクラス(時系列クラス)のオブジェクトになる。plot(x) とすれば期待通りの折れ線グラフを描いてくれる。plot(aggregate(x, FUN=mean)) とすれば年単位の平均値にまとめて描いてくれる。boxplot(x ~ cycle(x)) とすれば月ごとの箱ひげ図を描いてくれる。


Last modified: