ill-identifide氏がRの挙動についてブログにてRの挙動にてもとめたい結果と出力された結果が間違っていることに気づきにくい場合
久々に R の ... に騙されたので結果が正しいかどうか確認しないとダメなアンチパターン集を書く - ill-identified diary
について書いており、そこで
具体例を解説するのが面倒になってきたので他にも箇条書きで簡単に例を書いてみる. わかる人はわかるはず.
- ある変数の1日前の値を比較したくて
lag()を使用してラグ変数の列を作ったところ, 全期間で全く変化がないという結果になった. (ヒント:dplyrは使っていない)- 回帰モデルの一般化逆行列が特異 (= 完全多重共線性がある) である場合, 最小二乗法では係数を求められない. しかし最尤推定とみなしてニュートン法などを使うと, 結果自体は得られる. これは何を意味するのか? その結果を使用しても良いのか?
- R の
lm()の場合は特異でも結果が出てしまうことがある. この挙動を説明できるか?- 予測モデルを作って予測残差をプロットしたところ, 分布が偏っていたので, 目的変数を
log(y) ~ .のように対数変換 すればうまく当てはまりそうだと考えやり直してみた. しかしpredict()で得られた予測値の RMSE を計算したところ, 以前よりかなり大きな値が出たため対数変換に効果がないと考えた. (これは回答に複数のパターンがありうる)
とブログ中で3つ具体例を解説していない例があり、1つ目について考えてみた1。
{stats}のlag()の挙動に関する問題
多分、状況としては次のようなデータの処理を意図している気がする。
# 7日分程度の時系列データのサンプルデータを作る # ダミーデータでも良かったが、今回は東京都の2021年4月29日から5月5日の間のCovid-19の新規感染者発表数を使う covid.tokyo <- data.frame(Date = seq.Date(as.Date("2021-04-29"), as.Date("2021-05-05"), by = 1), New_Infections = c(1027, 698, 1050, 879, 708, 609, 621)) # こんなありがちな感じの時系列データができる covid.tokyo
## Date New_Infections
## 1 2021-04-29 1027
## 2 2021-04-30 698
## 3 2021-05-01 1050
## 4 2021-05-02 879
## 5 2021-05-03 708
## 6 2021-05-04 609
## 7 2021-05-05 621
ここで前の日の値をデータフレーム中に{stats}のlag()で作りたいと思ったりしたとする。lag()で-1個データをずらしたデータを作りたいと思ったら次のようなコードをイメージしそう。2
covid.tokyo.dame <- transform(covid.tokyo, Previous_Day_Infections = stats::lag(covid.tokyo$New_Infections, k = -1))
エラーが出ないで処理が終了する。しかし中身は、
covid.tokyo.dame
## Date New_Infections Previous_Day_Infections
## 1 2021-04-29 1027 1027
## 2 2021-04-30 698 698
## 3 2021-05-01 1050 1050
## 4 2021-05-02 879 879
## 5 2021-05-03 708 708
## 6 2021-05-04 609 609
## 7 2021-05-05 621 621
となり、新規感染者数を示すNew_Infectionsの列と、前の日の感染者数を意図したPrevious_Day_Infectionsの列が同じ値と意図した挙動をしていない。
ill-identifide氏がブログで言いたかったのではたぶん、ここで意図した値じゃない列ができているのにちゃんと気づこうよという話だと思う。誤ってここでできたPrevious_Day_InfectionsとNew_Infectionsで差をとって値が0になったのだけを見て、「前の日と差がなかったんだな。」などと解釈することがないようにということだと思う。
これくらい小さいデータや簡単な処理だとすぐに気付けるが、大きいデータを処理していて中身の確認を怠ったり、続けてなにか処理を続けているとエラーメッセージが出ていないので見落とすかもしれない。
ちなみに、{dplyr}のlag()を使うと似たようなコードでもちゃんとずれる。
library(dplyr)
##
## 次のパッケージを付け加えます: 'dplyr'
## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter, lag
## 以下のオブジェクトは 'package:base' からマスクされています:
##
## intersect, setdiff, setequal, union
covid.tokyo.dplyr <- covid.tokyo %>% mutate(Previous_Day_incfections = dplyr::lag(New_Infections, n = 1)) covid.tokyo.dplyr
## Date New_Infections Previous_Day_incfections
## 1 2021-04-29 1027 NA
## 2 2021-04-30 698 1027
## 3 2021-05-01 1050 698
## 4 2021-05-02 879 1050
## 5 2021-05-03 708 879
## 6 2021-05-04 609 708
## 7 2021-05-05 621 609
なんで{stats}のlag()では{dplyr}のlag()のような結果を得られないのかと言うことについて考えると、{stats}のlag()のドキュメントに、
Description
Compute a lagged version of a time series, shifting the time base back by a given number of observations.
Arguments
x
A vector or matrix or univariate or multivariate time series
と書いてあるように、取り扱うデータは時系列データを取り扱うことを想定している。時系列データとさらっと書いてあるが、これはRのtsオブジェクトという特殊な時系列データを取り扱うことを想定している。3
ちゃんと
A time series object with the same class as x.
とドキュメントにも書いている。
ざっくり{stats}のlag()の挙動を追いかける
ここで{stats}のlag()のソースコードをprint.function()で見てみると次のようになっている。
# {stats}のlag()はジェネリック関数なので、{stats}のlag.defalt()を直接呼び出す print.function(stats:::lag.default)
## function (x, k = 1, ...)
## {
## if (k != round(k)) {
## k <- round(k)
## warning("'k' is not an integer")
## }
## x <- hasTsp(x)
## p <- tsp(x)
## tsp(x) <- p - (k/p[3L]) * c(1, 1, 0)
## x
## }
## <bytecode: 0x0000000015e16c60>
## <environment: namespace:stats>
処理を追ってみる。
xがシフトしたいオブジェクトで、kがどの程度シフトさせたいのかという引数になっている。
最初のif文は、与えられたkが整数ではない場合はround()で四捨五入（五捨五入）し整数化しつつ、kが整数ではなかったことを注意喚起してくれている。
次のx <- hasTsp(x)はxがtsp属性を有しているか確認するもので、tsp属性を持っていない場合はc(1, NROW(x), 1)（NROW(X)はxの数）のtsp属性を付与する。
処理は至ってシンプルで次のように、if文でtsp属性が付与されているかいないか確認し、なければ上記のように属性を付与し、あればそのままの値を返す。
print.function(hasTsp)
## function (x)
## {
## if (is.null(attr(x, "tsp")))
## attr(x, "tsp") <- c(1, NROW(x), 1)
## x
## }
## <bytecode: 0x0000000015e192b8>
## <environment: namespace:stats>
このtsp属性はtsオブジェクトでは、時系列データの最初の値がStart(最初の観測の時)、次がEnd(最後の観測の時)、最後がFrequency（単位時間あたりの観測数）と解釈される。
例えば{datasets}パッケージにある、tsオブジェクトであるldeathsを見ると次のようになっている。
# イギリスの月別肺疾患死亡者数データ
ldeaths
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1974 3035 2552 2704 2554 2014 1655 1721 1524 1596 2074 2199 2512
## 1975 2933 2889 2938 2497 1870 1726 1607 1545 1396 1787 2076 2837
## 1976 2787 3891 3179 2011 1636 1580 1489 1300 1356 1653 2013 2823
## 1977 3102 2294 2385 2444 1748 1554 1498 1361 1346 1564 1640 2293
## 1978 2815 3137 2679 1969 1870 1633 1529 1366 1357 1570 1535 2491
## 1979 3084 2605 2573 2143 1693 1504 1461 1354 1333 1492 1781 1915
# tsp属性を見てみる tsp(ldeaths)
## [1] 1974.000 1979.917 12.000
# カレンダー形式じゃなく、StartとEndとFrequencyと値で表示する print(ldeaths, calendar = FALSE)
## Time Series:
## Start = c(1974, 1)
## End = c(1979, 12)
## Frequency = 12
## [1] 3035 2552 2704 2554 2014 1655 1721 1524 1596 2074 2199 2512 2933 2889 2938
## [16] 2497 1870 1726 1607 1545 1396 1787 2076 2837 2787 3891 3179 2011 1636 1580
## [31] 1489 1300 1356 1653 2013 2823 3102 2294 2385 2444 1748 1554 1498 1361 1346
## [46] 1564 1640 2293 2815 3137 2679 1969 1870 1633 1529 1366 1357 1570 1535 2491
## [61] 3084 2605 2573 2143 1693 1504 1461 1354 1333 1492 1781 1915
となっている。逆に、tsp属性が与えられていればclassをtsにするとtsオブジェクトができる。
lag()の話にもどして、次のp <- tsp(x)はそんなtsp属性の3つの値をpに代入しておいているだけ。
そしてtsp(x) <- p - (k/p[3L]) * c(1, 1, 0)はxのtsp属性の値からkをFrequencyで割った値の値分だけStartとEndの値を小さくする、つまり始まりと終わりの「観測時間」をシフトしたい分だけ観測頻度も考慮して前にずらすという処理をしている。
そして、最後のxで始まりと終わりの「観測時間」がシフトされた値が出力されているという流れになっている。
ここまでの一連の流れを見てわかるように、最初に与えられたxの値をずらしているのではなく、あくまで観測時間を決めていいる属性の値をずらして、ラグをだしているという処理になっている。
このlag()の処理はtsオブジェクトであるかどうかを確認する箇所がないので、ただの値だけのベクトルも受け付けてしまうが、値そのものをずらすのではなくtsオブジェクトの値の観測時間をいじっているだけなので、出力される値はそのままになってしまっている。しかも、tsでなければあまり意味のない属性を付与している割に、クラスはいじらないので、属性がただ付与されているだけなので、ただのベクトルを入力するとただのベクトルになぞの属性が付与されて出力されてくるだけになっている。
# 一見すると謎な属性が付与されるだけ stats::lag(covid.tokyo$New_Infections, k = -1)
## [1] 1027 698 1050 879 708 609 621
## attr(,"tsp")
## [1] 2 8 1
クラスをtsにするとこれも意味をなす。
as.ts(stats::lag(covid.tokyo$New_Infections, k = -1))
## Time Series:
## Start = 2
## End = 8
## Frequency = 1
## [1] 1027 698 1050 879 708 609 621
# pitnt.tsでcalendar = TRUEとして強制的に観測時間を付与する print(as.ts(stats::lag(covid.tokyo$New_Infections, k = -1)), calendar = TRUE)
## 2 3 4 5 6 7 8
## 1027 698 1050 879 708 609 621
# ちなみにもとのcovid.tokyo$New_Infectionsをただtsオブジェクトにしたものをみてみると print(as.ts(covid.tokyo$New_Infections), calendar = TRUE)
## 1 2 3 4 5 6 7
## 1027 698 1050 879 708 609 621
と観測時間が2からにシフトしている。ただし観測値自体はいじられていないのでこの値だけを取り出すともとの値のままになっている。
{stats}のlag()がなぜtsオブジェクト以外も受け付けるのかはまではよくわからないし、受け付けたとしてなんらかのメッセージを出してもいいんじゃないかとも思うが、そもそも論ちゃんとドキュメントよめばわかるし出力された結果をちゃんと確認しろといわれたらぐうの音も出ない。
たぶん、例題の答え？としてはこんな感じかと思う。
蛇足
じゃあ、{dplyr}とか使わないでずらすにはどうしたらいいかというとよくわからないが、適当に頭にNAつけてずらしてやればいいんじゃないかと思う。
transform(covid.tokyo, Previous_Day_Infections = c(NA, covid.tokyo$New_Infection[1:nrow(covid.tokyo)-1]))
## Date New_Infections Previous_Day_Infections
## 1 2021-04-29 1027 NA
## 2 2021-04-30 698 1027
## 3 2021-05-01 1050 698
## 4 2021-05-02 879 1050
## 5 2021-05-03 708 879
## 6 2021-05-04 609 708
## 7 2021-05-05 621 609
とりあえず1個ずれた。日付が欠損していたら予め日付を埋めるように処理しておけばよいかもしれないし、そもそもどういう処理をしたいのかによると思う。
さらなる蛇足
{zoo}パッケージのlag()ならちょっと工夫すればいける。zooのlag.zoo()はstats::lag()で呼べる（ジェネリック関数でzooオブジェクトが入れられるとそちらの処理に回る）。
require(zoo)
## 要求されたパッケージ zoo をロード中です
##
## 次のパッケージを付け加えます: 'zoo'
## 以下のオブジェクトは 'package:base' からマスクされています:
##
## as.Date, as.Date.numeric
transform(covid.tokyo, Previous_Day_Infections = stats::lag(zoo(covid.tokyo$New_Infections), k = -1, na.pad = TRUE))
## Date New_Infections Previous_Day_Infections
## 1 2021-04-29 1027 NA
## 2 2021-04-30 698 1027
## 3 2021-05-01 1050 698
## 4 2021-05-02 879 1050
## 5 2021-05-03 708 879
## 6 2021-05-04 609 708
## 7 2021-05-05 621 609
