統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

summaryの罠

年月を6桁の数値(YYYYMM)で表すために以下のように書いて何気なくsummaryを実行したところ、思わぬ挙動となった。

> Year  <- rep(2012:2016, each=12)
> Month <- rep(1:12, 5)
> YM    <- Year * 100 + Month
> summary(YM)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 201200  201300  201400  201400  201500  201600 

上記の通り、月は1〜12の範囲なのにいずれも下2桁が00となっている。
念のため他の関数を使ってみると、ちゃんとした値が返る。

> fivenum(YM)
[1] 201201.0 201303.5 201406.5 201509.5 201612.0
> min(YM)
[1] 201201
> max(YM)
[1] 201612

summaryはgeneric functionなので、今回の例のようなnumeric型に対するsummaryの挙動を確認するためにmethodを調べてみると:

> methods("summary")
 [1] summary.aov                    summary.aovlist*               summary.aspell*               
 [4] summary.check_packages_in_dir* summary.connection             summary.data.frame            
 [7] summary.Date                   summary.default                summary.ecdf*                 
[10] summary.factor                 summary.glm                    summary.infl*                 
[13] summary.lm                     summary.loess*                 summary.manova                
[16] summary.matrix                 summary.mlm*                   summary.nls*                  
[19] summary.packageStatus*         summary.PDF_Dictionary*        summary.PDF_Stream*           
[22] summary.POSIXct                summary.POSIXlt                summary.ppr*                  
[25] summary.prcomp*                summary.princomp*              summary.proc_time             
[28] summary.srcfile                summary.srcref                 summary.stepfun               
[31] summary.stl*                   summary.table                  summary.tukeysmooth* 

色々と出てきたが、どうやらsummary.defaultがそれっぽい。
summary.defaultの中身を調べてみる。

> summary.default
…
中略
…
      }) else if (is.numeric(object)) {
         nas <- is.na(object)
         object <- object[!nas]
         qq <- stats::quantile(object)
         qq <- signif(c(qq[1L:3L], mean(object), qq[4L:5L]), digits)
         names(qq) <- c("Min.", "1st Qu.", "Median", "Mean", "3rd Qu.", 
                        "Max.")
         if (any(nas)) 
            c(qq, `NA's` = sum(nas))
         else qq
      } 
…
中略
…

どうやらsummary.defaultの中身は{stats}のquantileで、間にmeanを挟んでいるだけの様子。そこでquantileを直接実行してみると:

> quantile(YM)
      0%      25%      50%      75%     100% 
201201.0 201303.8 201406.5 201509.2 201612.0 

それっぽい値が返ってきている。どうやら次の処理のsignifの引数のdigitsが怪しい。
そこで今度はsummaryのdigitsを指定して実行してみる。

> summary(YM, digits = 8)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
201201.0 201303.8 201406.5 201406.5 201509.2 201612.0 

狙った通りの挙動となった。どうやら6桁の整数では大きすぎて小数点以下が丸められてしまっていたようだ。
めでたしめでたし・・・・と思ったら、よく見るとsummary(つまりquantile)の返り値とfivenumの返り値が微妙に違っている。

> summary(YM, digits = 8)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
201201.0 201303.8 201406.5 201406.5 201509.2 201612.0 
> fivenum(YM)
[1] 201201.0 201303.5 201406.5 201509.5 201612.0

25パーセンタイルと75パーセンタイルの値が違う。調べてみると、fivenumは四分位数ではなくヒンジを返すのだそう。

http://cse.naro.affrc.go.jp/takezawa/r-tips/r/59.html

なのでデータが偶数個の場合に先ほどのような挙動となる様子。
こんな感じ。

> hoge <- 1:10
> summary(hoge)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    3.25    5.50    5.50    7.75   10.00 
> quantile(hoge)
   0%   25%   50%   75%  100% 
 1.00  3.25  5.50  7.75 10.00 
> fivenum(hoge)
[1]  1.0  3.0  5.5  8.0 10.0 # 合わない!
> foo <- c(0, hoge)
> summary(foo)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     2.5     5.0     5.0     7.5    10.0 
> quantile(foo)
  0%  25%  50%  75% 100% 
 0.0  2.5  5.0  7.5 10.0 
> fivenum(foo)
[1]  0.0  2.5  5.0  7.5 10.0 # 一致した!